    -----------------------------------------------------------------------------
    -- |
    -- 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 :: forall a. Floating a => 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
    ^ :: forall a. Num 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' :: forall a. (Ord a, Num a) => 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 :: forall d. (Floating d, Ord d) => 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' :: forall d. (Floating d, Ord d) => 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 :: forall d. (Floating d, Ord d) => 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' :: forall d. (Floating d, Ord d) => 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 :: forall d. (Floating d, Ord d) => 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