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

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

import Diagrams.Util (tau)

------------------------------------------------------------
------------------------------------------------------------

-- | 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     = []

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