    {-# OPTIONS_GHC -fno-warn-unused-binds #-}
    -----------------------------------------------------------------------------
    -- |
    -- Module      :  Diagrams.Solve
    -- Copyright   :  (c) 2011 diagrams-lib team (see LICENSE)
    -- License     :  BSD-style (see LICENSE)
    -- Maintainer  :  diagrams-discuss@googlegroups.com
    --
    -- Exact solving of low-degree (n <= 3) polynomials.
    --
    -----------------------------------------------------------------------------
    module Diagrams.Solve
           ( quadForm
           , cubForm
           ) where

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

    import           Diagrams.Util (tau)

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

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

        -- There are infinitely many solutions in this case,
        -- so arbitrarily return 0
      | a == 0 && b == 0 && c == 0 = [0]

        -- c /= 0
      | a == 0 && b == 0 = []

        -- linear
      | a == 0    = [-c/b]

        -- no real solutions
      | d < 0     = []

        -- ax^2 + c = 0
      | b == 0    = [sqrt (-c/a), -sqrt (-c/a)]

        -- multiplicity 2 solution
      | d == 0    = [-b/(2*a)]

        -- 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
           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

    ------------------------------------------------------------
    -- 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.
    cubForm :: (Floating d, Ord d) => d -> d -> d -> d -> [d]
    cubForm a b c d
      | aboutZero a             = quadForm b c d

        -- three real roots, use trig method to avoid complex numbers
      | delta >  0              = map trig [0,1,2]

        -- one real root of multiplicity 3
      | delta == 0 && disc == 0 = [ -b/(3*a) ]

        -- 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)
                                  ]

        -- 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]
                  | 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)
           trig k = 2 * sqrt(-p/3) * cos(1/3*acos(3*q/(2*p)*sqrt(-3/p)) - 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

    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
               -- 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?