diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-08 08:48:12 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-08 08:48:12 +0200 |
commit | 1925c123d7d8184a1d2ddc0a413e0fd2776e1083 (patch) | |
tree | fad79f909d9c3be53d68e6ebd67202650536d387 /packages/hmatrix/src/Numeric/GSL/Polynomials.hs | |
parent | eb3f702d065a4a967bb754977233e6eec408fd1f (diff) |
empty hmatrix-base
Diffstat (limited to 'packages/hmatrix/src/Numeric/GSL/Polynomials.hs')
-rw-r--r-- | packages/hmatrix/src/Numeric/GSL/Polynomials.hs | 58 |
1 files changed, 58 insertions, 0 deletions
diff --git a/packages/hmatrix/src/Numeric/GSL/Polynomials.hs b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs new file mode 100644 index 0000000..290c615 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs | |||
@@ -0,0 +1,58 @@ | |||
1 | {-# LANGUAGE CPP, ForeignFunctionInterface #-} | ||
2 | ----------------------------------------------------------------------------- | ||
3 | {- | | ||
4 | Module : Numeric.GSL.Polynomials | ||
5 | Copyright : (c) Alberto Ruiz 2006 | ||
6 | License : GPL-style | ||
7 | |||
8 | Maintainer : Alberto Ruiz (aruiz at um dot es) | ||
9 | Stability : provisional | ||
10 | Portability : uses ffi | ||
11 | |||
12 | Polynomials. | ||
13 | |||
14 | <http://www.gnu.org/software/gsl/manual/html_node/General-Polynomial-Equations.html#General-Polynomial-Equations> | ||
15 | |||
16 | -} | ||
17 | ----------------------------------------------------------------------------- | ||
18 | module Numeric.GSL.Polynomials ( | ||
19 | polySolve | ||
20 | ) where | ||
21 | |||
22 | import Data.Packed.Internal | ||
23 | import Data.Complex | ||
24 | import System.IO.Unsafe (unsafePerformIO) | ||
25 | |||
26 | #if __GLASGOW_HASKELL__ >= 704 | ||
27 | import Foreign.C.Types (CInt(..)) | ||
28 | #endif | ||
29 | |||
30 | {- | Solution of general polynomial equations, using /gsl_poly_complex_solve/. | ||
31 | |||
32 | For example, the three solutions of x^3 + 8 = 0 | ||
33 | |||
34 | >>> polySolve [8,0,0,1] | ||
35 | [(-2.0) :+ 0.0,1.0 :+ 1.7320508075688776,1.0 :+ (-1.7320508075688776)] | ||
36 | |||
37 | |||
38 | The example in the GSL manual: To find the roots of x^5 -1 = 0: | ||
39 | |||
40 | >>> polySolve [-1, 0, 0, 0, 0, 1] | ||
41 | [(-0.8090169943749472) :+ 0.5877852522924731, | ||
42 | (-0.8090169943749472) :+ (-0.5877852522924731), | ||
43 | 0.30901699437494756 :+ 0.9510565162951535, | ||
44 | 0.30901699437494756 :+ (-0.9510565162951535), | ||
45 | 1.0000000000000002 :+ 0.0] | ||
46 | |||
47 | -} | ||
48 | polySolve :: [Double] -> [Complex Double] | ||
49 | polySolve = toList . polySolve' . fromList | ||
50 | |||
51 | polySolve' :: Vector Double -> Vector (Complex Double) | ||
52 | polySolve' v | dim v > 1 = unsafePerformIO $ do | ||
53 | r <- createVector (dim v-1) | ||
54 | app2 c_polySolve vec v vec r "polySolve" | ||
55 | return r | ||
56 | | otherwise = error "polySolve on a polynomial of degree zero" | ||
57 | |||
58 | foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TVCV | ||