    {-# LANGUAGE CPP #-}
    {- |
    Module      :  Numeric.GSL.Polynomials
    Copyright   :  (c) Alberto Ruiz 2006
    License     :  GPL
    Maintainer  :  Alberto Ruiz
    Stability   :  provisional

    Polynomials.

    <http://www.gnu.org/software/gsl/manual/html_node/General-Polynomial-Equations.html#General-Polynomial-Equations>

    -}

    module Numeric.GSL.Polynomials (
        polySolve
    ) where

    import Numeric.LinearAlgebra.HMatrix
    import Numeric.GSL.Internal
    import System.IO.Unsafe (unsafePerformIO)

    #if __GLASGOW_HASKELL__ >= 704
    import Foreign.C.Types (CInt(..))
    #endif

    {- | Solution of general polynomial equations, using /gsl_poly_complex_solve/.

    For example, the three solutions of x^3 + 8 = 0

    >>> polySolve [8,0,0,1]
    [(-2.0) :+ 0.0,1.0 :+ 1.7320508075688776,1.0 :+ (-1.7320508075688776)]

    The example in the GSL manual: To find the roots of x^5 -1 = 0:

    >>> polySolve [-1, 0, 0, 0, 0, 1]
    [(-0.8090169943749472) :+ 0.5877852522924731,
    (-0.8090169943749472) :+ (-0.5877852522924731),
    0.30901699437494756 :+ 0.9510565162951535,
    0.30901699437494756 :+ (-0.9510565162951535),
    1.0000000000000002 :+ 0.0]

    -}
    polySolve :: [Double] -> [Complex Double]
    polySolve :: [Double] -> [Complex Double]
    polySolve = Vector (Complex Double) -> [Complex Double]
    forall a. Storable a => Vector a -> [a]
    toList (Vector (Complex Double) -> [Complex Double])
    -> ([Double] -> Vector (Complex Double))
    -> [Double]
    -> [Complex Double]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Double -> Vector (Complex Double)
    polySolve' (Vector Double -> Vector (Complex Double))
    -> ([Double] -> Vector Double)
    -> [Double]
    -> Vector (Complex Double)
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList

    polySolve' :: Vector Double -> Vector (Complex Double)
    polySolve' :: Vector Double -> Vector (Complex Double)
    polySolve' Vector Double
    v | Vector Double -> IndexOf Vector
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size Vector Double
    v Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    > Int
    1 = IO (Vector (Complex Double)) -> Vector (Complex Double)
    forall a. IO a -> a
    unsafePerformIO (IO (Vector (Complex Double)) -> Vector (Complex Double))
    -> IO (Vector (Complex Double)) -> Vector (Complex Double)
    forall a b. (a -> b) -> a -> b
    $ do
        r <- Int -> IO (Vector (Complex Double))
    forall a. Storable a => Int -> IO (Vector a)
    createVector (Vector Double -> IndexOf Vector
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size Vector Double
    vInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1)
        (v `applyRaw` (r `applyRaw` id)) c_polySolve #| "polySolve"
        return r
                 | Bool
    otherwise = String -> Vector (Complex Double)
    forall a. HasCallStack => String -> a
    error String
    "polySolve on a polynomial of degree zero"

    foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TV (TCV Res)