{-# 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)