From 7f486ddf956b2728b1ff92bb69f95b532aea5526 Mon Sep 17 00:00:00 2001 From: Mathnerd314 Date: Sat, 17 May 2014 21:12:11 -0600 Subject: [PATCH 1/4] Quartic formula, no obvious bugs --- src/Diagrams/Solve.hs | 64 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 58 insertions(+), 6 deletions(-) diff --git a/src/Diagrams/Solve.hs b/src/Diagrams/Solve.hs index 09751c2a..72f9ee78 100644 --- a/src/Diagrams/Solve.hs +++ b/src/Diagrams/Solve.hs @@ -6,12 +6,13 @@ -- License : BSD-style (see LICENSE) -- Maintainer : diagrams-discuss@googlegroups.com -- --- Exact solving of low-degree (n <= 3) polynomials. +-- Exact solving of low-degree (n <= 4) polynomials. -- ----------------------------------------------------------------------------- module Diagrams.Solve ( quadForm , cubForm + , quartForm ) where import Data.List (maximumBy) @@ -92,12 +93,10 @@ cubForm a b c d xx = 2*b*b*b - 9*a*b*c + 27*a*a*d p = disc/(3*a*a) q = xx/(27*a*a*a) - trig k = 2 * sqrt(-p/3) * cos(1/3*acos(3*q/(2*p)*sqrt(-3/p)) - k*tau/3) - - b/(3*a) - + phi = 1/3*acos(3*q/(2*p)*sqrt(-3/p)) + trig k = 2 * sqrt(-p/3) * cos(phi - k*tau/3) - b/(3*a) cubert x | x < 0 = -((-x)**(1/3)) | otherwise = x**(1/3) - aboutZero x = abs x < toler toler = 1e-10 @@ -115,4 +114,57 @@ cubForm_prop a b c d = all (aboutZero . eval) (cubForm a b c d) -- with numerical stability. If this turns out to be an -- issue in practice we could, say, use the solutions -- generated here as very good guesses to a numerical - -- solver which can give us a more precise answer? \ No newline at end of file + -- solver which can give us a more precise answer? + +------------------------------------------------------------ +-- Quartic formula +------------------------------------------------------------ + +-- Based on http://tog.acm.org/resources/GraphicsGems/gems/Roots3b/and4.c +-- as of 5/12/14, with help from http://en.wikipedia.org/wiki/Quartic_function + +-- | Solve the quartic equation c4 x^4 + c3 x^3 + c2 x^2 + c1 x + c0 = 0, returning a +-- list of all real roots. +quartForm :: (Floating d, Ord d) => d -> d -> d -> d -> d -> [d] +quartForm c4 c3 c2 c1 c0 + -- obvious cubic + | aboutZero c4 = cubForm c3 c2 c1 c0 + -- x(ax^3+bx^2+cx+d) + | aboutZero c0 = 0 : cubForm c4 c3 c2 c1 + | otherwise = + let -- eliminate c4: x^4+ax^3+bx^2+cx+d + a=c3/c4; b=c2/c4; c=c1/c4; d=c0/c4 + -- substitute x = u - a/4 to eliminate cubic term: + -- u^4 + pu^2 + qu + r = 0 + -- + p = b - 3/8*a*a + q = 1/8*a*a*a-a*b/2+c + r = (-3/256)*a*a*a*a+a*a*b/16-a*c/4+d + in + map (\x->x-(a/4)) $ if aboutZero r then + -- no constant term: u(u^3 + pu + q) = 0 + 0 : cubForm 1 0 p q + else do + -- solve the resolvent cubic - should have exactly one solution, but cubic solver is inaccurate + let z:_ = cubForm 1 (-p/2) (-r) (p*r/2 - q*q/8) + -- quadratic equations + let u = z * z - r + v = 2 * z - p + if u < 0 || v < 0 then [] else + let u' = if aboutZero u then 0 else sqrt u + v' = if aboutZero v then 0 else sqrt v + s1 = quadForm 1 (if' (q<0) (-v') v') (z-u') + s2 = quadForm 1 (if' (q<0) v' (-v')) (z+u') + in + s1++s2 + where + aboutZero x = abs x < toler + toler = 1e-10 + if' a b c = if a then b else c + +quartForm_prop :: Double -> Double -> Double -> Double -> Double -> Bool +quartForm_prop a b c d e = all (aboutZero . eval) (quartForm a b c d e) + where eval x = a*x*x*x*x + b*x*x*x + c*x*x + d*x + e + aboutZero x = abs x < tolerance + tolerance = 1e-5 + -- Same note about tolerance as for cubic From 4291c4f0e31582ced5a658fbddc187cd033a0b68 Mon Sep 17 00:00:00 2001 From: Mathnerd314 Date: Sun, 18 May 2014 14:25:59 -0600 Subject: [PATCH 2/4] Various nits and cleanups --- src/Diagrams/Solve.hs | 141 +++++++++++++++++++++++------------------- 1 file changed, 78 insertions(+), 63 deletions(-) diff --git a/src/Diagrams/Solve.hs b/src/Diagrams/Solve.hs index 72f9ee78..189147ba 100644 --- a/src/Diagrams/Solve.hs +++ b/src/Diagrams/Solve.hs @@ -1,4 +1,3 @@ -{-# OPTIONS_GHC -fno-warn-unused-binds #-} ----------------------------------------------------------------------------- -- | -- Module : Diagrams.Solve @@ -20,6 +19,20 @@ import Data.Ord (comparing) import Diagrams.Util (tau) +import Prelude hiding ((^)) +import qualified Prelude as P ((^)) + +-- | A specialization of (^) to Integer +-- c.f. http://comments.gmane.org/gmane.comp.lang.haskell.libraries/21164 +-- for discussion. "The choice in (^) and (^^) to overload on the +-- power's Integral type... was a genuinely bad idea." - Edward Kmett +(^) :: (Num a) => a -> Integer -> a +(^) = (P.^) + +-- | Utility function used to avoid singularities +aboutZero' :: (Ord a, Num a) => a -> a -> Bool +aboutZero' toler x = abs x < toler + ------------------------------------------------------------ -- Quadratic formula ------------------------------------------------------------ @@ -49,14 +62,12 @@ quadForm a b c -- see http://www.mpi-hd.mpg.de/astrophysik/HEA/internal/Numerical_Recipes/f5-6.pdf | otherwise = [q/a, c/q] - where d = b*b - 4*a*c + where d = b^2 - 4*a*c q = -1/2*(b + signum b * sqrt d) quadForm_prop :: Double -> Double -> Double -> Bool -quadForm_prop a b c = all (aboutZero . eval) (quadForm a b c) - where eval x = a*x*x + b*x + c - aboutZero x = abs x < tolerance - tolerance = 1e-10 +quadForm_prop a b c = all (aboutZero' 1e-10 . eval) (quadForm a b c) + where eval x = a*x^2 + b*x + c ------------------------------------------------------------ -- Cubic formula @@ -65,10 +76,10 @@ quadForm_prop a b c = all (aboutZero . eval) (quadForm a b c) -- See http://en.wikipedia.org/wiki/Cubic_formula#General_formula_of_roots -- | Solve the cubic equation ax^3 + bx^2 + cx + d = 0, returning a --- list of all real roots. -cubForm :: (Floating d, Ord d) => d -> d -> d -> d -> [d] -cubForm a b c d - | aboutZero a = quadForm b c d +-- list of all real roots. First argument is tolerance. +cubForm' :: (Floating d, Ord d) => d -> d -> d -> d -> d -> [d] +cubForm' toler a b c d + | aboutZero' toler a = quadForm b c d -- three real roots, use trig method to avoid complex numbers | delta > 0 = map trig [0,1,2] @@ -78,33 +89,35 @@ cubForm a b c d -- two real roots, one of multiplicity 2 | delta == 0 && disc /= 0 = [ (b*c - 9*a*d)/(2*disc) - , (9*a*a*d - 4*a*b*c + b*b*b)/(a * disc) + , (9*a^2*d - 4*a*b*c + b^3)/(a * disc) ] -- one real root (and two complex) | otherwise = [-b/(3*a) - cc/(3*a) + disc/(3*a*cc)] - where delta = 18*a*b*c*d - 4*b*b*b*d + b*b*c*c - 4*a*c*c*c - 27*a*a*d*d - disc = 3*a*c - b*b - qq = sqrt(-27*a*a*delta) - qq' | aboutZero disc = maximumBy (comparing (abs . (+xx))) [qq, -qq] + where delta = 18*a*b*c*d - 4*b^3*d + b^2*c^2 - 4*a*c^3 - 27*a^2*d^2 + disc = 3*a*c - b^2 + qq = sqrt(-27*(a^2)*delta) + qq' | aboutZero' toler disc = maximumBy (comparing (abs . (+xx))) [qq, -qq] | otherwise = qq cc = cubert (1/2*(qq' + xx)) - xx = 2*b*b*b - 9*a*b*c + 27*a*a*d - p = disc/(3*a*a) - q = xx/(27*a*a*a) + xx = 2*b^3 - 9*a*b*c + 27*a^2*d + p = disc/(3*a^2) + q = xx/(27*a^3) phi = 1/3*acos(3*q/(2*p)*sqrt(-3/p)) trig k = 2 * sqrt(-p/3) * cos(phi - k*tau/3) - b/(3*a) cubert x | x < 0 = -((-x)**(1/3)) | otherwise = x**(1/3) - aboutZero x = abs x < toler - toler = 1e-10 + +-- | Solve the cubic equation ax^3 + bx^2 + cx + d = 0, returning a +-- list of all real roots within 1e-10 tolerance +-- (although currently it's closer to 1e-5) +cubForm :: (Floating d, Ord d) => d -> d -> d -> d -> [d] +cubForm = cubForm' 1e-10 cubForm_prop :: Double -> Double -> Double -> Double -> Bool -cubForm_prop a b c d = all (aboutZero . eval) (cubForm a b c d) - where eval x = a*x*x*x + b*x*x + c*x + d - aboutZero x = abs x < tolerance - tolerance = 1e-5 +cubForm_prop a b c d = all (aboutZero' 1e-5 . eval) (cubForm a b c d) + where eval x = a*x^3 + b*x^2 + c*x + d -- Basically, however large you set the tolerance it seems -- that quickcheck can always come up with examples where -- the returned solutions evaluate to something near zero @@ -124,47 +137,49 @@ cubForm_prop a b c d = all (aboutZero . eval) (cubForm a b c d) -- as of 5/12/14, with help from http://en.wikipedia.org/wiki/Quartic_function -- | Solve the quartic equation c4 x^4 + c3 x^3 + c2 x^2 + c1 x + c0 = 0, returning a --- list of all real roots. -quartForm :: (Floating d, Ord d) => d -> d -> d -> d -> d -> [d] -quartForm c4 c3 c2 c1 c0 +-- list of all real roots. First argument is tolerance. +quartForm' :: (Floating d, Ord d) => d -> d -> d -> d -> d -> d -> [d] +quartForm' toler c4 c3 c2 c1 c0 -- obvious cubic - | aboutZero c4 = cubForm c3 c2 c1 c0 + | aboutZero' toler c4 = cubForm c3 c2 c1 c0 -- x(ax^3+bx^2+cx+d) - | aboutZero c0 = 0 : cubForm c4 c3 c2 c1 - | otherwise = - let -- eliminate c4: x^4+ax^3+bx^2+cx+d - a=c3/c4; b=c2/c4; c=c1/c4; d=c0/c4 - -- substitute x = u - a/4 to eliminate cubic term: - -- u^4 + pu^2 + qu + r = 0 - -- - p = b - 3/8*a*a - q = 1/8*a*a*a-a*b/2+c - r = (-3/256)*a*a*a*a+a*a*b/16-a*c/4+d - in - map (\x->x-(a/4)) $ if aboutZero r then - -- no constant term: u(u^3 + pu + q) = 0 - 0 : cubForm 1 0 p q - else do - -- solve the resolvent cubic - should have exactly one solution, but cubic solver is inaccurate - let z:_ = cubForm 1 (-p/2) (-r) (p*r/2 - q*q/8) - -- quadratic equations - let u = z * z - r - v = 2 * z - p - if u < 0 || v < 0 then [] else - let u' = if aboutZero u then 0 else sqrt u - v' = if aboutZero v then 0 else sqrt v - s1 = quadForm 1 (if' (q<0) (-v') v') (z-u') - s2 = quadForm 1 (if' (q<0) v' (-v')) (z+u') - in - s1++s2 - where - aboutZero x = abs x < toler - toler = 1e-10 - if' a b c = if a then b else c + | aboutZero' toler c0 = 0 : cubForm c4 c3 c2 c1 + -- substitute solutions of y back to x + | otherwise = map (\x->x-(a/4)) roots + where + -- eliminate c4: x^4+ax^3+bx^2+cx+d + [a,b,c,d] = map (/c4) [c3,c2,c1,c0] + -- eliminate cubic term via x = y - a/4 + -- reduced quartic: y^4 + py^2 + qy + r = 0 + p = b - 3/8*a^2 + q = 1/8*a^3-a*b/2+c + r = (-3/256)*a^4+a^2*b/16-a*c/4+d + + -- | roots of the reduced quartic + roots | aboutZero' toler r = + 0 : cubForm 1 0 p q -- no constant term: y(y^3 + py + q) = 0 + | u < 0 || v < 0 = [] -- no real solutions due to square root + | otherwise = s1++s2 -- solutions of the quadratics + + -- solve the resolvent cubic - only one solution is needed + z:_ = cubForm 1 (-p/2) (-r) (p*r/2 - q^2/8) + + -- solve the two quadratic equations + -- y^2 ± v*y-(±u-z) + u = z^2 - r + v = 2*z - p + u' = if aboutZero' toler u then 0 else sqrt u + v' = if aboutZero' toler v then 0 else sqrt v + s1 = quadForm 1 (if q<0 then -v' else v') (z-u') + s2 = quadForm 1 (if q<0 then v' else -v') (z+u') + +-- | Solve the quartic equation c4 x^4 + c3 x^3 + c2 x^2 + c1 x + c0 = 0, returning a +-- list of all real roots within 1e-10 tolerance +-- (although currently it's closer to 1e-5) +quartForm :: (Floating d, Ord d) => d -> d -> d -> d -> d -> [d] +quartForm = quartForm' 1e-10 quartForm_prop :: Double -> Double -> Double -> Double -> Double -> Bool -quartForm_prop a b c d e = all (aboutZero . eval) (quartForm a b c d e) - where eval x = a*x*x*x*x + b*x*x*x + c*x*x + d*x + e - aboutZero x = abs x < tolerance - tolerance = 1e-5 +quartForm_prop a b c d e = all (aboutZero' 1e-5 . eval) (quartForm a b c d e) + where eval x = a*x^4 + b*x^3 + c*x^2 + d*x + e -- Same note about tolerance as for cubic From 07f41be6534b598e0f7e2fed6a6ae2d56e3d879b Mon Sep 17 00:00:00 2001 From: Mathnerd314 Date: Sun, 18 May 2014 17:27:41 -0600 Subject: [PATCH 3/4] Export cubForm' and quartForm' --- src/Diagrams/Solve.hs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Diagrams/Solve.hs b/src/Diagrams/Solve.hs index 189147ba..69b5f424 100644 --- a/src/Diagrams/Solve.hs +++ b/src/Diagrams/Solve.hs @@ -12,6 +12,8 @@ module Diagrams.Solve ( quadForm , cubForm , quartForm + , cubForm' + , quartForm' ) where import Data.List (maximumBy) From dd6f273073fa712a051e3f6932fbb87e2b1092ea Mon Sep 17 00:00:00 2001 From: Mathnerd314 Date: Sun, 18 May 2014 18:19:39 -0600 Subject: [PATCH 4/4] Add underscores to silence warnings, and add warnings to match with Travis CI --- diagrams-lib.cabal | 1 + src/Diagrams/Solve.hs | 12 ++++++------ 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/diagrams-lib.cabal b/diagrams-lib.cabal index a0a06be5..33362cfd 100644 --- a/diagrams-lib.cabal +++ b/diagrams-lib.cabal @@ -115,4 +115,5 @@ Library if impl(ghc < 7.6) Build-depends: ghc-prim Hs-source-dirs: src + ghc-options: -Wall default-language: Haskell2010 diff --git a/src/Diagrams/Solve.hs b/src/Diagrams/Solve.hs index 69b5f424..f027548d 100644 --- a/src/Diagrams/Solve.hs +++ b/src/Diagrams/Solve.hs @@ -67,8 +67,8 @@ quadForm a b c where d = b^2 - 4*a*c q = -1/2*(b + signum b * sqrt d) -quadForm_prop :: Double -> Double -> Double -> Bool -quadForm_prop a b c = all (aboutZero' 1e-10 . eval) (quadForm a b c) +_quadForm_prop :: Double -> Double -> Double -> Bool +_quadForm_prop a b c = all (aboutZero' 1e-10 . eval) (quadForm a b c) where eval x = a*x^2 + b*x + c ------------------------------------------------------------ @@ -117,8 +117,8 @@ cubForm' toler a b c d cubForm :: (Floating d, Ord d) => d -> d -> d -> d -> [d] cubForm = cubForm' 1e-10 -cubForm_prop :: Double -> Double -> Double -> Double -> Bool -cubForm_prop a b c d = all (aboutZero' 1e-5 . eval) (cubForm a b c d) +_cubForm_prop :: Double -> Double -> Double -> Double -> Bool +_cubForm_prop a b c d = all (aboutZero' 1e-5 . eval) (cubForm a b c d) where eval x = a*x^3 + b*x^2 + c*x + d -- Basically, however large you set the tolerance it seems -- that quickcheck can always come up with examples where @@ -181,7 +181,7 @@ quartForm' toler c4 c3 c2 c1 c0 quartForm :: (Floating d, Ord d) => d -> d -> d -> d -> d -> [d] quartForm = quartForm' 1e-10 -quartForm_prop :: Double -> Double -> Double -> Double -> Double -> Bool -quartForm_prop a b c d e = all (aboutZero' 1e-5 . eval) (quartForm a b c d e) +_quartForm_prop :: Double -> Double -> Double -> Double -> Double -> Bool +_quartForm_prop a b c d e = all (aboutZero' 1e-5 . eval) (quartForm a b c d e) where eval x = a*x^4 + b*x^3 + c*x^2 + d*x + e -- Same note about tolerance as for cubic