-----------------------------------------------------------------------------
-- |
-- Module      :  Diagrams.Solve.Polynomial
-- Copyright   :  (c) 2011-2015 diagrams-solve team (see LICENSE)
-- License     :  BSD-style (see LICENSE)
-- Maintainer  :  diagrams-discuss@googlegroups.com
--
-- Exact solving of low-degree (n <= 4) polynomials.
--
-----------------------------------------------------------------------------
module Diagrams.Solve.Polynomial
       ( quadForm
       , cubForm
       , quartForm
       , cubForm'
       , quartForm'
       ) where

import           Data.List (maximumBy)
import           Data.Ord  (comparing)

import           Prelude   hiding ((^))
import qualified Prelude   as P ((^))

-- | The fundamental circle constant, /i.e./ ratio between a circle's
--   circumference and radius.
tau :: Floating a => a
tau :: a
tau = a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
forall a. Floating a => a
pi

-- | 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
--
--   Note there are rewrite rules in GHC.Real to expand small exponents.
(^) :: (Num a) => a -> Integer -> a
^ :: a -> Integer -> a
(^) = a -> Integer -> a
forall a b. (Num a, Integral b) => a -> b -> a
(P.^)

-- | Utility function used to avoid singularities
aboutZero' :: (Ord a, Num a) => a -> a -> Bool
aboutZero' :: a -> a -> Bool
aboutZero' a
toler a
x = a -> a
forall a. Num a => a -> a
abs a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
toler
{-# INLINE aboutZero' #-}

------------------------------------------------------------
-- Quadratic formula
------------------------------------------------------------

-- | The quadratic formula.
quadForm :: (Floating d, Ord d) => d -> d -> d -> [d]
quadForm :: d -> d -> d -> [d]
quadForm d
a d
b d
c

    -- There are infinitely many solutions in this case,
    -- so arbitrarily return 0
  | d
a d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0 Bool -> Bool -> Bool
&& d
b d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0 Bool -> Bool -> Bool
&& d
c d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0 = [d
0]

    -- c /= 0
  | d
a d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0 Bool -> Bool -> Bool
&& d
b d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0 = []

    -- linear
  | d
a d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0    = [-d
cd -> d -> d
forall a. Fractional a => a -> a -> a
/d
b]

    -- no real solutions
  | d
d d -> d -> Bool
forall a. Ord a => a -> a -> Bool
< d
0     = []

    -- ax^2 + c = 0
  | d
b d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0    = [d -> d
forall a. Floating a => a -> a
sqrt (-d
cd -> d -> d
forall a. Fractional a => a -> a -> a
/d
a), -d -> d
forall a. Floating a => a -> a
sqrt (-d
cd -> d -> d
forall a. Fractional a => a -> a -> a
/d
a)]

    -- multiplicity 2 solution
  | d
d d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0    = [-d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/(d
2d -> d -> d
forall a. Num a => a -> a -> a
*d
a)]

    -- see http://www.mpi-hd.mpg.de/astrophysik/HEA/internal/Numerical_Recipes/f5-6.pdf
  | Bool
otherwise = [d
qd -> d -> d
forall a. Fractional a => a -> a -> a
/d
a, d
cd -> d -> d
forall a. Fractional a => a -> a -> a
/d
q]
 where d :: d
d = d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2 d -> d -> d
forall a. Num a => a -> a -> a
- d
4d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
c
       q :: d
q = -d
1d -> d -> d
forall a. Fractional a => a -> a -> a
/d
2d -> d -> d
forall a. Num a => a -> a -> a
*(d
b d -> d -> d
forall a. Num a => a -> a -> a
+ d -> d
forall a. Num a => a -> a
signum d
b d -> d -> d
forall a. Num a => a -> a -> a
* d -> d
forall a. Floating a => a -> a
sqrt d
d)
{-# INLINE quadForm #-}

_quadForm_prop :: Double -> Double -> Double -> Bool
_quadForm_prop :: Double -> Double -> Double -> Bool
_quadForm_prop Double
a Double
b Double
c = (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' Double
1e-10 (Double -> Bool) -> (Double -> Double) -> Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
eval) (Double -> Double -> Double -> [Double]
forall d. (Floating d, Ord d) => d -> d -> d -> [d]
quadForm Double
a Double
b Double
c)
  where eval :: Double -> Double
eval Double
x = Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^Integer
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c

------------------------------------------------------------
-- Cubic formula
------------------------------------------------------------

-- 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. First argument is tolerance.
cubForm' :: (Floating d, Ord d) => d -> d -> d -> d -> d -> [d]
cubForm' :: d -> d -> d -> d -> d -> [d]
cubForm' d
toler d
a d
b d
c d
d
  | d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
a      = d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> [d]
quadForm d
b d
c d
d

    -- three real roots, use trig method to avoid complex numbers
  | d
delta d -> d -> Bool
forall a. Ord a => a -> a -> Bool
>  d
0              = (d -> d) -> [d] -> [d]
forall a b. (a -> b) -> [a] -> [b]
map d -> d
trig [d
0,d
1,d
2]

    -- one real root of multiplicity 3
  | d
delta d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0 Bool -> Bool -> Bool
&& d
disc d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0 = [ -d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/(d
3d -> d -> d
forall a. Num a => a -> a -> a
*d
a) ]

    -- two real roots, one of multiplicity 2
  | d
delta d -> d -> Bool
forall a. Eq a => a -> a -> Bool
== d
0 Bool -> Bool -> Bool
&& d
disc d -> d -> Bool
forall a. Eq a => a -> a -> Bool
/= d
0 = [ (d
bd -> d -> d
forall a. Num a => a -> a -> a
*d
c d -> d -> d
forall a. Num a => a -> a -> a
- d
9d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
d)d -> d -> d
forall a. Fractional a => a -> a -> a
/(d
2d -> d -> d
forall a. Num a => a -> a -> a
*d
disc)
                              , (d
9d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2d -> d -> d
forall a. Num a => a -> a -> a
*d
d d -> d -> d
forall a. Num a => a -> a -> a
- d
4d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> d -> d
forall a. Num a => a -> a -> a
*d
c d -> d -> d
forall a. Num a => a -> a -> a
+ d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
3)d -> d -> d
forall a. Fractional a => a -> a -> a
/(d
a d -> d -> d
forall a. Num a => a -> a -> a
* d
disc)
                              ]

    -- one real root (and two complex)
  | Bool
otherwise               = [-d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/(d
3d -> d -> d
forall a. Num a => a -> a -> a
*d
a) d -> d -> d
forall a. Num a => a -> a -> a
- d
ccd -> d -> d
forall a. Fractional a => a -> a -> a
/(d
3d -> d -> d
forall a. Num a => a -> a -> a
*d
a) d -> d -> d
forall a. Num a => a -> a -> a
+ d
discd -> d -> d
forall a. Fractional a => a -> a -> a
/(d
3d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
cc)]

 where delta :: d
delta  = d
18d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> d -> d
forall a. Num a => a -> a -> a
*d
cd -> d -> d
forall a. Num a => a -> a -> a
*d
d d -> d -> d
forall a. Num a => a -> a -> a
- d
4d -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
3d -> d -> d
forall a. Num a => a -> a -> a
*d
d d -> d -> d
forall a. Num a => a -> a -> a
+ d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2d -> d -> d
forall a. Num a => a -> a -> a
*d
cd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2 d -> d -> d
forall a. Num a => a -> a -> a
- d
4d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
cd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
3 d -> d -> d
forall a. Num a => a -> a -> a
- d
27d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2d -> d -> d
forall a. Num a => a -> a -> a
*d
dd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2
       disc :: d
disc   = d
3d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
c d -> d -> d
forall a. Num a => a -> a -> a
- d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2
       qq :: d
qq     = d -> d
forall a. Floating a => a -> a
sqrt(-d
27d -> d -> d
forall a. Num a => a -> a -> a
*(d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2)d -> d -> d
forall a. Num a => a -> a -> a
*d
delta)
       qq' :: d
qq'    = if d -> d
forall a. Num a => a -> a
abs (d
xx d -> d -> d
forall a. Num a => a -> a -> a
+ d
qq) d -> d -> Bool
forall a. Ord a => a -> a -> Bool
> d -> d
forall a. Num a => a -> a
abs (d
xx d -> d -> d
forall a. Num a => a -> a -> a
- d
qq) then d
qq else -d
qq
       cc :: d
cc     = d -> d
forall a. (Ord a, Floating a) => a -> a
cubert (d
1d -> d -> d
forall a. Fractional a => a -> a -> a
/d
2d -> d -> d
forall a. Num a => a -> a -> a
*(d
qq' d -> d -> d
forall a. Num a => a -> a -> a
+ d
xx))
       xx :: d
xx     = d
2d -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
3 d -> d -> d
forall a. Num a => a -> a -> a
- d
9d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> d -> d
forall a. Num a => a -> a -> a
*d
c d -> d -> d
forall a. Num a => a -> a -> a
+ d
27d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2d -> d -> d
forall a. Num a => a -> a -> a
*d
d
       p :: d
p      = d
discd -> d -> d
forall a. Fractional a => a -> a -> a
/(d
3d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2)
       q :: d
q      = d
xxd -> d -> d
forall a. Fractional a => a -> a -> a
/(d
27d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
3)
       phi :: d
phi = d
1d -> d -> d
forall a. Fractional a => a -> a -> a
/d
3d -> d -> d
forall a. Num a => a -> a -> a
*d -> d
forall a. Floating a => a -> a
acos(d
3d -> d -> d
forall a. Num a => a -> a -> a
*d
qd -> d -> d
forall a. Fractional a => a -> a -> a
/(d
2d -> d -> d
forall a. Num a => a -> a -> a
*d
p)d -> d -> d
forall a. Num a => a -> a -> a
*d -> d
forall a. Floating a => a -> a
sqrt(-d
3d -> d -> d
forall a. Fractional a => a -> a -> a
/d
p))
       trig :: d -> d
trig d
k = d
2 d -> d -> d
forall a. Num a => a -> a -> a
* d -> d
forall a. Floating a => a -> a
sqrt(-d
pd -> d -> d
forall a. Fractional a => a -> a -> a
/d
3) d -> d -> d
forall a. Num a => a -> a -> a
* d -> d
forall a. Floating a => a -> a
cos(d
phi d -> d -> d
forall a. Num a => a -> a -> a
- d
kd -> d -> d
forall a. Num a => a -> a -> a
*d
forall a. Floating a => a
taud -> d -> d
forall a. Fractional a => a -> a -> a
/d
3) d -> d -> d
forall a. Num a => a -> a -> a
- d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/(d
3d -> d -> d
forall a. Num a => a -> a -> a
*d
a)
       cubert :: a -> a
cubert a
x | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0     = -((-a
x)a -> a -> a
forall a. Floating a => a -> a -> a
**(a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
3))
                | Bool
otherwise = a
xa -> a -> a
forall a. Floating a => a -> a -> a
**(a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
3)
{-# INLINE cubForm' #-}

-- | 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 :: d -> d -> d -> d -> [d]
cubForm = d -> d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> d -> [d]
cubForm' d
1e-10
{-# INLINE cubForm #-}

_cubForm_prop :: Double -> Double -> Double -> Double -> Bool
_cubForm_prop :: Double -> Double -> Double -> Double -> Bool
_cubForm_prop Double
a Double
b Double
c Double
d = (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' Double
1e-5 (Double -> Bool) -> (Double -> Double) -> Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
eval) (Double -> Double -> Double -> Double -> [Double]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm Double
a Double
b Double
c Double
d)
  where eval :: Double -> Double
eval Double
x = Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^Integer
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^Integer
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
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
           -- but larger than the tolerance (but it takes it more
           -- tries the larger you set the tolerance). Wonder if this
           -- is an inherent limitation or (more likely) a problem
           -- 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?

------------------------------------------------------------
-- 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. First argument is tolerance.
quartForm' :: (Floating d, Ord d) => d -> d -> d -> d -> d -> d -> [d]
quartForm' :: d -> d -> d -> d -> d -> d -> [d]
quartForm' d
toler d
c4 d
c3 d
c2 d
c1 d
c0
  -- obvious cubic
  | d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
c4 = d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm d
c3 d
c2 d
c1 d
c0
  -- x(ax^3+bx^2+cx+d)
  | d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
c0 = d
0 d -> [d] -> [d]
forall a. a -> [a] -> [a]
: d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm d
c4 d
c3 d
c2 d
c1
  -- substitute solutions of y back to x
  | Bool
otherwise = (d -> d) -> [d] -> [d]
forall a b. (a -> b) -> [a] -> [b]
map (\d
x->d
xd -> d -> d
forall a. Num a => a -> a -> a
-(d
ad -> d -> d
forall a. Fractional a => a -> a -> a
/d
4)) [d]
roots
    where
      -- eliminate c4: x^4+ax^3+bx^2+cx+d
      [d
a,d
b,d
c,d
d] = (d -> d) -> [d] -> [d]
forall a b. (a -> b) -> [a] -> [b]
map (d -> d -> d
forall a. Fractional a => a -> a -> a
/d
c4) [d
c3,d
c2,d
c1,d
c0]
      -- eliminate cubic term via x = y - a/4
      -- reduced quartic: y^4 + py^2 + qy + r = 0
      p :: d
p = d
b d -> d -> d
forall a. Num a => a -> a -> a
- d
3d -> d -> d
forall a. Fractional a => a -> a -> a
/d
8d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2
      q :: d
q = d
1d -> d -> d
forall a. Fractional a => a -> a -> a
/d
8d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
3d -> d -> d
forall a. Num a => a -> a -> a
-d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/d
2d -> d -> d
forall a. Num a => a -> a -> a
+d
c
      r :: d
r = (-d
3d -> d -> d
forall a. Fractional a => a -> a -> a
/d
256)d -> d -> d
forall a. Num a => a -> a -> a
*d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
4d -> d -> d
forall a. Num a => a -> a -> a
+d
ad -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2d -> d -> d
forall a. Num a => a -> a -> a
*d
bd -> d -> d
forall a. Fractional a => a -> a -> a
/d
16d -> d -> d
forall a. Num a => a -> a -> a
-d
ad -> d -> d
forall a. Num a => a -> a -> a
*d
cd -> d -> d
forall a. Fractional a => a -> a -> a
/d
4d -> d -> d
forall a. Num a => a -> a -> a
+d
d

      -- | roots of the reduced quartic
      roots :: [d]
roots | d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
r =
                d
0 d -> [d] -> [d]
forall a. a -> [a] -> [a]
: d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm d
1 d
0 d
p d
q   -- no constant term: y(y^3 + py + q) = 0
            | d
u d -> d -> Bool
forall a. Ord a => a -> a -> Bool
< -d
toler Bool -> Bool -> Bool
|| d
v d -> d -> Bool
forall a. Ord a => a -> a -> Bool
< -d
toler = []     -- no real solutions due to square root
            | Bool
otherwise      = [d]
s1[d] -> [d] -> [d]
forall a. [a] -> [a] -> [a]
++[d]
s2 -- solutions of the quadratics

      -- solve the resolvent cubic - only one solution is needed
      d
z:[d]
_ = d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> [d]
cubForm d
1 (-d
pd -> d -> d
forall a. Fractional a => a -> a -> a
/d
2) (-d
r) (d
pd -> d -> d
forall a. Num a => a -> a -> a
*d
rd -> d -> d
forall a. Fractional a => a -> a -> a
/d
2 d -> d -> d
forall a. Num a => a -> a -> a
- d
qd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2d -> d -> d
forall a. Fractional a => a -> a -> a
/d
8)

      -- solve the two quadratic equations
      -- y^2 ± v*y-(±u-z)
      u :: d
u = d
zd -> Integer -> d
forall a. Num a => a -> Integer -> a
^Integer
2 d -> d -> d
forall a. Num a => a -> a -> a
- d
r
      v :: d
v = d
2d -> d -> d
forall a. Num a => a -> a -> a
*d
z d -> d -> d
forall a. Num a => a -> a -> a
- d
p
      u' :: d
u' = if d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
u then d
0 else d -> d
forall a. Floating a => a -> a
sqrt d
u
      v' :: d
v' = if d -> d -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' d
toler d
v then d
0 else d -> d
forall a. Floating a => a -> a
sqrt d
v
      s1 :: [d]
s1 = d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> [d]
quadForm d
1 (if d
qd -> d -> Bool
forall a. Ord a => a -> a -> Bool
<d
0 then -d
v' else d
v') (d
zd -> d -> d
forall a. Num a => a -> a -> a
-d
u')
      s2 :: [d]
s2 = d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> [d]
quadForm d
1 (if d
qd -> d -> Bool
forall a. Ord a => a -> a -> Bool
<d
0 then d
v' else -d
v') (d
zd -> d -> d
forall a. Num a => a -> a -> a
+d
u')
{-# INLINE quartForm' #-}

-- | 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 :: d -> d -> d -> d -> d -> [d]
quartForm = d -> d -> d -> d -> d -> d -> [d]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> d -> d -> [d]
quartForm' d
1e-10
{-# INLINE quartForm #-}

_quartForm_prop :: Double -> Double -> Double -> Double -> Double -> Bool
_quartForm_prop :: Double -> Double -> Double -> Double -> Double -> Bool
_quartForm_prop Double
a Double
b Double
c Double
d Double
e = (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. (Ord a, Num a) => a -> a -> Bool
aboutZero' Double
1e-5 (Double -> Bool) -> (Double -> Double) -> Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
eval) (Double -> Double -> Double -> Double -> Double -> [Double]
forall d. (Floating d, Ord d) => d -> d -> d -> d -> d -> [d]
quartForm Double
a Double
b Double
c Double
d Double
e)
  where eval :: Double -> Double
eval Double
x = Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^Integer
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^Integer
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Integer -> Double
forall a. Num a => a -> Integer -> a
^Integer
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e
           -- Same note about tolerance as for cubic