Skip to content

Commit

Permalink
Various nits and cleanups
Browse files Browse the repository at this point in the history
  • Loading branch information
@Mathnerd314
Mathnerd314 committed May 18, 2014
1 parent 7f486dd commit 4291c4f
Showing 1 changed file with 78 additions and 63 deletions.
141 changes: 78 additions & 63 deletions src/Diagrams/Solve.hs
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
{-# OPTIONS_GHC -fno-warn-unused-binds #-}
-----------------------------------------------------------------------------
-- |
-- Module : Diagrams.Solve
Expand All @@ -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
------------------------------------------------------------
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand All @@ -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
Expand All @@ -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

1 comment on commit 4291c4f

@jeffreyrosenbluth
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good !

Please sign in to comment.