hmatrix 0.2.0.0 → 0.2.1.0
raw patch · 21 files changed
+1079/−928 lines, 21 filesdep +HUnitdep +QuickCheckdep +storable-complex
Dependencies added: HUnit, QuickCheck, storable-complex
Files
- README +16/−6
- examples/experiments/bigtests.hs +60/−0
- examples/pca1.hs +1/−1
- examples/tests.hs +0/−506
- hmatrix.cabal +12/−10
- lib/Data/Packed/Internal/Common.hs +1/−9
- lib/Numeric/GSL/Special/Airy.hs +2/−2
- lib/Numeric/GSL/Special/Debye.hs +16/−16
- lib/Numeric/GSL/Special/Ellint.hs +16/−16
- lib/Numeric/GSL/Special/Legendre.hs +91/−73
- lib/Numeric/GSL/Special/auto.hs +2/−0
- lib/Numeric/GSL/Special/autoall.sh +37/−38
- lib/Numeric/GSL/Special/gsl_sf_exp.h +0/−146
- lib/Numeric/GSL/Special/gsl_sf_legendre.h +319/−0
- lib/Numeric/GSL/Special/gsl_sf_log.h +0/−90
- lib/Numeric/GSL/Special/replace.hs +14/−0
- lib/Numeric/LinearAlgebra/Algorithms.hs +2/−2
- lib/Numeric/LinearAlgebra/Linear.hs +2/−13
- lib/Numeric/LinearAlgebra/Tests.hs +197/−0
- lib/Numeric/LinearAlgebra/Tests/Instances.hs +142/−0
- lib/Numeric/LinearAlgebra/Tests/Properties.hs +149/−0
README view
@@ -33,20 +33,30 @@ INSTALLATION -------------------------------------- -$ runhaskell Setup.lhs configure --prefix=$HOME --user-$ runhaskell Setup.lhs build-$ runhaskell Setup.lhs haddock-$ runhaskell Setup.lhs install+Automatic (using cabal-install and HackageDB): + $ cabal install hmatrix++Manual:++ Install storable-complex from HackageDB and then++ $ runhaskell Setup.lhs configure --prefix=$HOME --user+ $ runhaskell Setup.lhs build+ $ runhaskell Setup.lhs haddock+ $ runhaskell Setup.lhs install+ See below for installation on Windows. TESTS --------------------------------------------- -$ runhaskell examples/tests+$ ghci+GHCi, version 6.8.2: http://www.haskell.org/ghc/ :? for help+Prelude> Numeric.LinearAlgebra.Tests.runTests 20 Additional tests with big matrices (taking a few minutes): -$ runhaskell examples/tests --big+$ runhaskell examples/experiments bigtests EXAMPLES ------------------------------------------------------
+ examples/experiments/bigtests.hs view
@@ -0,0 +1,60 @@+module Main where++import Numeric.LinearAlgebra+import Numeric.LinearAlgebra.Tests+import System.Random(randomRs,mkStdGen)+import Test.HUnit hiding (test)+import System(getArgs)+++pseudorandomR seed (n,m) = reshape m $ fromList $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed++pseudorandomC seed (n,m) = toComplex (pseudorandomR seed (n,m), pseudorandomR (seed+1) (n,m))++bigmat = m + trans m+ where m = pseudorandomR 18 (1000,1000) :: Matrix Double+bigmatc = mc + ctrans mc+ where mc = pseudorandomC 19 (1000,1000) :: Matrix (Complex Double)++utest str b = TestCase $ assertBool str b++feye n = flipud (ident n) :: Matrix Double++infixl 4 |~|+a |~| b = dist a b < 10^^(-10)++dist a b = r+ where norm = pnorm Infinity+ na = norm a+ nb = norm b+ nab = norm (a-b)+ mx = max na nb+ mn = min na nb+ r = if mn < eps+ then mx+ else nab/mx++square m = rows m == cols m++unitary m = square m && m <> ctrans m |~| ident (rows m)++eigProp m = complex m <> v |~| v <> diag s+ where (s, v) = eig m++eigSHProp m = m <> v |~| v <> real (diag s)+ && unitary v+ && m |~| v <> real (diag s) <> ctrans v+ where (s, v) = eigSH m++bigtests = do+ putStrLn "--------- big matrices -----"+ runTestTT $ TestList+ [ utest "eigS" $ eigSHProp bigmat+ , utest "eigH" $ eigSHProp bigmatc+ , utest "eigR" $ eigProp bigmat+ , utest "eigC" $ eigProp bigmatc+ , utest "det" $ det (feye 1000) == 1 && det (feye 1002) == -1+ ]+ return ()++main = bigtests
examples/pca1.hs view
@@ -26,7 +26,7 @@ decode x = x <> vp + m m = mean dataSet c = cov dataSet- (_,v) = eigSH c+ (_,v) = eigSH' c vp = takeRows n (trans v) norm = pnorm PNorm2
− examples/tests.hs
@@ -1,506 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}--module Main where--import Numeric.GSL hiding (sin,cos,exp,choose)-import Numeric.LinearAlgebra-import Numeric.LinearAlgebra.LAPACK-import qualified Numeric.GSL.Matrix as GSL-import Test.QuickCheck hiding (test)-import Test.HUnit hiding ((~:),test)-import System.Random(randomRs,mkStdGen)-import System.Info-import Data.List(foldl1', transpose)-import System(getArgs)-import Debug.Trace(trace)--debug x = trace (show x) x--type RM = Matrix Double-type CM = Matrix (Complex Double)---- relative error-dist :: (Normed t, Num t) => t -> t -> Double-dist a b = r- where norm = pnorm Infinity- na = norm a- nb = norm b- nab = norm (a-b)- mx = max na nb- mn = min na nb- r = if mn < eps- then mx- else nab/mx--infixl 4 |~|-a |~| b = a :~10~: b--data Aprox a = (:~) a Int--(~:) :: (Normed a, Num a) => Aprox a -> a -> Bool-a :~n~: b = dist a b < 10^^(-n)---maxdim = 10--instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where- arbitrary = do- r <- arbitrary- i <- arbitrary- return (r:+i)- coarbitrary = undefined--instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where - arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)- m <- choose (1,maxdim)- n <- choose (1,maxdim)- l <- vector (m*n)- ctype <- arbitrary- let h = if ctype then (m><n) else (m>|<n)- trMode <- arbitrary- let tr = if trMode then trans else id- return $ tr (h l)- coarbitrary = undefined--data PairM a = PairM (Matrix a) (Matrix a) deriving Show-instance (Num a, Element a, Arbitrary a) => Arbitrary (PairM a) where- arbitrary = do- a <- choose (1,maxdim)- b <- choose (1,maxdim)- c <- choose (1,maxdim)- l1 <- vector (a*b)- l2 <- vector (b*c)- return $ PairM ((a><b) (map fromIntegral (l1::[Int]))) ((b><c) (map fromIntegral (l2::[Int])))- --return $ PairM ((a><b) l1) ((b><c) l2)- coarbitrary = undefined--data SqM a = SqM (Matrix a) deriving Show-sqm (SqM a) = a-instance (Element a, Arbitrary a) => Arbitrary (SqM a) where- arbitrary = do- n <- choose (1,maxdim)- l <- vector (n*n)- return $ SqM $ (n><n) l- coarbitrary = undefined--data Sym a = Sym (Matrix a) deriving Show-sym (Sym a) = a-instance (Linear Vector a, Arbitrary a) => Arbitrary (Sym a) where- arbitrary = do- SqM m <- arbitrary- return $ Sym (m + trans m)- coarbitrary = undefined--data Her = Her (Matrix (Complex Double)) deriving Show-her (Her a) = a-instance {-(Field a, Arbitrary a, Num a) =>-} Arbitrary Her where- arbitrary = do- SqM m <- arbitrary- return $ Her (m + ctrans m)- coarbitrary = undefined--data PairSM a = PairSM (Matrix a) (Matrix a) deriving Show-instance (Num a, Field a, Arbitrary a) => Arbitrary (PairSM a) where- arbitrary = do- a <- choose (1,maxdim)- c <- choose (1,maxdim)- l1 <- vector (a*a)- l2 <- vector (a*c)- return $ PairSM ((a><a) (map fromIntegral (l1::[Int]))) ((a><c) (map fromIntegral (l2::[Int])))- --return $ PairSM ((a><a) l1) ((a><c) l2)- coarbitrary = undefined--instance (Field a, Arbitrary a) => Arbitrary (Vector a) where - arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)- m <- choose (1,maxdim^2)- l <- vector m- return $ fromList l- coarbitrary = undefined--data PairV a = PairV (Vector a) (Vector a)-instance (Field a, Arbitrary a) => Arbitrary (PairV a) where - arbitrary = do --m <- sized $ \max -> choose (1,1+3*max)- m <- choose (1,maxdim^2)- l1 <- vector m- l2 <- vector m- return $ PairV (fromList l1) (fromList l2)- coarbitrary = undefined--------------------------------------------------------------------------test str b = TestCase $ assertBool str b--------------------------------------------------------------------------pseudorandomR seed (n,m) = reshape m $ fromList $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed --pseudorandomC seed (n,m) = toComplex (pseudorandomR seed (n,m), pseudorandomR (seed+1) (n,m))--bigmat = m + trans m :: RM- where m = pseudorandomR 18 (1000,1000)-bigmatc = mc + ctrans mc ::CM- where mc = pseudorandomC 19 (1000,1000)---------------------------------------------------------------------------m = (3><3)- [ 1, 2, 3- , 4, 5, 7- , 2, 8, 4 :: Double- ]--mc = (3><3)- [ 1, 2, 3- , 4, 5, 7- , 2, 8, i- ]---mr = (3><4)- [ 1, 2, 3, 4,- 2, 4, 6, 8,- 1, 1, 1, 2:: Double- ]--mrc = (3><4)- [ 1, 2, 3, 4,- 2, 4, 6, 8,- i, i, i, 2- ]--a = (3><4)- [ 1, 0, 0, 0- , 0, 2, 0, 0- , 0, 0, 0, 0 :: Double- ]--b = (3><4)- [ 1, 0, 0, 0- , 0, 2, 3, 0- , 0, 0, 4, 0 :: Double- ]--ac = (2><3) [1 .. 6::Double]-bc = (3><4) [7 .. 18::Double]--af = (2>|<3) [1,4,2,5,3,6::Double]-bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double]-----------------------------------------------------------feye n = flipud (ident n) :: Matrix Double---luTest1 m = m |~| p <> l <> u- where (l,u,p,_) = lu m--detTest1 = det m == 26- && det mc == 38 :+ (-3)- && det (feye 2) == -1--detTest2 m = s d1 |~| s d2- where d1 = det m- d2 = det' m * det q- det' m = product $ toList $ takeDiag r- (q,r) = qr m- s x = fromList [x]--invTest m = degenerate m || m <> inv m |~| ident (rows m)--pinvTest m = m <> p <> m |~| m- && p <> m <> p |~| p- && hermitian (m<>p)- && hermitian (p<>m)- where p = pinv m--square m = rows m == cols m--unitary m = square m && m <> ctrans m |~| ident (rows m)--hermitian m = m |~| ctrans m--upperTriang m = rows m == 1 || down == z- where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))- z = constant 0 (dim down)--upperHessenberg m = rows m < 3 || down == z- where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))- z = constant 0 (dim down)--svdTest svd m = u <> real d <> trans v |~| m- && unitary u && unitary v- where (u,d,v) = full svd m--svdTest' svd m = m |~| 0 || u <> real (diag s) <> trans v |~| m- where (u,s,v) = economy svd m--eigTest m = complex m <> v |~| v <> diag s- where (s, v) = eig m--eigTestSH m = m <> v |~| v <> real (diag s)- && unitary v- && m |~| v <> real (diag s) <> ctrans v- where (s, v) = eigSH m--zeros (r,c) = reshape c (constant 0 (r*c))--ones (r,c) = zeros (r,c) + 1--degenerate m = rank m < min (rows m) (cols m) --prec = 1E-15--singular m = s1 < prec || s2/s1 < prec- where (_,ss,_) = svd m- s = toList ss- s1 = maximum s- s2 = minimum s--nullspaceTest m = null nl || m <> n |~| zeros (r,c) -- 0- where nl = nullspacePrec 1 m- n = fromColumns nl- r = rows m- c = cols m - rank m------------------------------------------------------------------------polyEval cs x = foldr (\c ac->ac*x+c) 0 cs--polySolveTest' p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))---polySolveTest = test "polySolve" (polySolveTest' [1,2,3,4])-------------------------------------------------------------------------quad f a b = fst $ integrateQAGS 1E-9 100 f a b---- A multiple integral can be easily defined using partial application-quad2 f a b g1 g2 = quad h a b- where h x = quad (f x) (g1 x) (g2 x)--volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) - 0 r (const 0) (\x->sqrt (r*r-x*x))--epsTol = 1E-8::Double--integrateTest = test "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < epsTol)-------------------------------------------------------------------------besselTest = test "bessel_J0_e" ( abs (r-expected) < e )- where (r,e) = bessel_J0_e 5.0- expected = -0.17759677131433830434739701--exponentialTest = test "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 )- where (v,e,err) = exp_e10_e 30.0- expected = exp 30.0--gammaTest = test "gamma" (gamma 5 == 24.0)-------------------------------------------------------------------------cholRTest = chol ((2><2) [1,2,2,9::Double]) == (2><2) [1,2,0,2.23606797749979]-cholCTest = chol ((2><2) [1,2,2,9::Complex Double]) == (2><2) [1,2,0,2.23606797749979]-------------------------------------------------------------------------qrTest qr m = q <> r |~| m && unitary q && upperTriang r- where (q,r) = qr m-------------------------------------------------------------------------hessTest m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h- where (p,h) = hess m-------------------------------------------------------------------------schurTest1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s- where (u,s) = schur m--schurTest2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme- where (u,s) = schur m-------------------------------------------------------------------------nd1 = (3><3) [ 1/2, 1/4, 1/4- , 0/1, 1/2, 1/4- , 1/2, 1/4, 1/2 :: Double]--nd2 = (2><2) [1, 0, 1, 1:: Complex Double]--expmTest1 = expm nd1 :~14~: (3><3)- [ 1.762110887278176- , 0.478085470590435- , 0.478085470590435- , 0.104719410945666- , 1.709751181805343- , 0.425725765117601- , 0.851451530235203- , 0.530445176063267- , 1.814470592751009 ]--expmTest2 = expm nd2 :~15~: (2><2)- [ 2.718281828459045- , 0.000000000000000- , 2.718281828459045- , 2.718281828459045 ]--expmTestDiag m = expm (logm m) |~| complex m- where logm m = matFunc Prelude.log m---------------------------------------------------------------------------asFortran m = (rows m >|< cols m) $ toList (flatten $ trans m)-asC m = (rows m >< cols m) $ toList (flatten m)--mulC a b = a <> b-mulF a b = trans $ trans b <> trans a-----------------------------------------------------------------------------multiplyG a b = reshape (cols b) $ fromList $ concat $ multiplyL (toLists a) (toLists b)- where multiplyL a b = [[dotL x y | y <- transpose b] | x <- a]- dotL a b = sum (zipWith (*) a b)--r >|< c = f where- f l | dim v == r*c = reshapeF r v- | otherwise = error "(>|<)"- where v = fromList l- reshapeF r = trans . reshape r-------------------------------------------------------------------------rot :: Double -> Matrix Double-rot a = (3><3) [ c,0,s- , 0,1,0- ,-s,0,c ]- where c = cos a- s = sin a--fun n = foldl1' (<>) (map rot angles)- where angles = toList $ linspace n (0,1)--rotTest = fun (10^5) :~12~: rot 5E4-------------------------------------------------------------------------tests = do- setErrorHandlerOff- putStrLn "--------- internal -----"- quickCheck ((\m -> m == trans m).sym :: Sym Double -> Bool)- quickCheck ((\m -> m == trans m).sym :: Sym (Complex Double) -> Bool)- quickCheck $ \l -> null l || (toList . fromList) l == (l :: [Double])- quickCheck $ \l -> null l || (toList . fromList) l == (l :: [Complex Double])- quickCheck $ \m -> m == asC (m :: RM)- quickCheck $ \m -> m == asC (m :: CM)- quickCheck $ \m -> m == asFortran (m :: RM)- quickCheck $ \m -> m == asFortran (m :: CM)- quickCheck $ \m -> m == (asC . asFortran) (m :: RM)- quickCheck $ \m -> m == (asC . asFortran) (m :: CM)- runTestTT $ TestList- [ test "1E5 rots" rotTest- ]- putStrLn "--------- multiply ----"- quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: RM)- quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: CM)- quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: RM))- quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: CM))- quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: RM)- quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: CM)- putStrLn "--------- lu ---------"- quickCheck (luTest1 :: RM->Bool)- quickCheck (luTest1 :: CM->Bool)- quickCheck (detTest2 . sqm :: SqM Double -> Bool)- quickCheck (detTest2 . sqm :: SqM (Complex Double) -> Bool)- runTestTT $ TestList- [ test "det1" detTest1- ]- putStrLn "--------- svd ---------"- quickCheck (svdTest svdR)- quickCheck (svdTest svdRdd)- quickCheck (svdTest svdC)- quickCheck (svdTest' svdR)- quickCheck (svdTest' svdRdd)- quickCheck (svdTest' svdC)- quickCheck (svdTest' GSL.svdg)- putStrLn "--------- eig ---------"- quickCheck (eigTest . sqm :: SqM Double -> Bool)- quickCheck (eigTest . sqm :: SqM (Complex Double) -> Bool)- quickCheck (eigTestSH . sym :: Sym Double -> Bool)- quickCheck (eigTestSH . her :: Her -> Bool)- putStrLn "--------- inv ------"- quickCheck (invTest . sqm :: SqM Double -> Bool)- quickCheck (invTest . sqm :: SqM (Complex Double) -> Bool)- putStrLn "--------- pinv ------"- quickCheck (pinvTest ::RM->Bool)- if os == "mingw32"- then putStrLn "complex pinvTest skipped in this OS"- else quickCheck (pinvTest ::CM->Bool)- putStrLn "--------- chol ------"- runTestTT $ TestList- [ test "cholR" cholRTest- , test "cholC" cholRTest- ]- putStrLn "--------- qr ---------"- quickCheck (qrTest GSL.qr)- quickCheck (qrTest (GSL.unpackQR . GSL.qrPacked))- quickCheck (qrTest ( unpackQR . GSL.qrPacked))- quickCheck (qrTest qr ::RM->Bool)- quickCheck (qrTest qr ::CM->Bool)- putStrLn "--------- hess --------"- quickCheck (hessTest . sqm ::SqM Double->Bool)- quickCheck (hessTest . sqm ::SqM (Complex Double) -> Bool)- putStrLn "--------- schur --------"- quickCheck (schurTest2 . sqm ::SqM Double->Bool)- if os == "mingw32"- then putStrLn "complex schur skipped in this OS"- else quickCheck (schurTest1 . sqm ::SqM (Complex Double) -> Bool)- putStrLn "--------- expm --------"- runTestTT $ TestList- [ test "expmd" (expmTestDiag $ (2><2) [1,2,3,5 :: Double])- , test "expm1" (expmTest1)- , test "expm2" (expmTest2)- ]- putStrLn "--------- nullspace ------"- quickCheck (nullspaceTest :: RM -> Bool)- quickCheck (nullspaceTest :: CM -> Bool)- putStrLn "--------- vector operations ------"- quickCheck $ (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM))- quickCheck $ (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))- quickCheck $ (\u -> cos u * tan u |~| sin (u::RM))- quickCheck $ (\u -> (cos u * tan u) :~6~: sin (u::CM))- runTestTT $ TestList- [ test "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)- , test "arith2" $ (((1+i) .* ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( (140*i-51).*1 :: CM)- , test "arith3" $ exp (i.*ones(10,10)*pi) + 1 |~| 0- , test "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3]- ]- putStrLn "--------- GSL ------"- quickCheck $ \v -> ifft (fft v) |~| v- runTestTT $ TestList- [ gammaTest- , besselTest- , exponentialTest- , integrateTest- , polySolveTest- ]--bigtests = do- putStrLn "--------- big matrices -----"- runTestTT $ TestList- [ test "eigS" $ eigTestSH bigmat- , test "eigH" $ eigTestSH bigmatc- , test "eigR" $ eigTest bigmat- , test "eigC" $ eigTest bigmatc- , test "det" $ det (feye 1000) == 1 && det (feye 1002) == -1- ]--main = do- args <- getArgs- if "--big" `elem` args- then bigtests- else tests
hmatrix.cabal view
@@ -1,5 +1,5 @@ Name: hmatrix-Version: 0.2.0.0+Version: 0.2.1.0 License: GPL License-file: LICENSE Author: Alberto Ruiz@@ -13,26 +13,25 @@ . More information: <http://alberrto.googlepages.com/gslhaskell> Category: Numerical, Math-tested-with: GHC ==6.6.1, GHC ==6.8.1, GHC ==6.8.2+tested-with: GHC ==6.8.2 cabal-version: >=1.2+build-type: Simple flag splitBase description: Choose the new smaller, split-up base package. library if flag(splitBase)- build-depends: base >= 3, array+ build-depends: base >= 3, array, QuickCheck, HUnit, storable-complex else- build-depends: base < 3-- ghc-options: -O+ build-depends: base < 3, QuickCheck, HUnit, storable-complex if os(windows)- ghc-options: -fasm+ ghc-options: else if arch(x86_64)- ghc-options: -fasm+ ghc-options: else ghc-options: -fvia-C @@ -84,12 +83,15 @@ Numeric.LinearAlgebra.Instances, Numeric.LinearAlgebra.Interface, Numeric.LinearAlgebra.Algorithms,- Graphics.Plot+ Graphics.Plot,+ Numeric.LinearAlgebra.Tests other-modules: Data.Packed.Internal, Data.Packed.Internal.Common, Data.Packed.Internal.Vector, Data.Packed.Internal.Matrix,- Numeric.GSL.Special.Internal+ Numeric.GSL.Special.Internal,+ Numeric.LinearAlgebra.Tests.Instances+ Numeric.LinearAlgebra.Tests.Properties C-sources: lib/Data/Packed/Internal/auxi.c, lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c, lib/Numeric/GSL/gsl-aux.c
lib/Data/Packed/Internal/Common.hs view
@@ -22,16 +22,8 @@ import Debug.Trace import Foreign.C.String(peekCString) import Foreign.C.Types+import Foreign.Storable.Complex ------------------------------------------------------------------------instance (Storable a, RealFloat a) => Storable (Complex a) where --- alignment x = alignment (realPart x) --- sizeOf x = 2 * sizeOf (realPart x) --- peek p = do --- [re,im] <- peekArray 2 (castPtr p) --- return (re :+ im) --- poke p (a :+ b) = pokeArray (castPtr p) [a,b] ------------------------------------------------------------------------- -- | @debug x = trace (show x) x@ debug :: (Show a) => a -> a
lib/Numeric/GSL/Special/Airy.hs view
@@ -38,8 +38,8 @@ , airy_zero_Ai_deriv_e , airy_zero_Ai_deriv , airy_zero_Bi_deriv_e-, airy_zero_Bi_deriv,-Precision(..)+, airy_zero_Bi_deriv+, Precision(..) ) where import Foreign(Ptr)
lib/Numeric/GSL/Special/Debye.hs view
@@ -23,10 +23,10 @@ , debye_3 , debye_4_e , debye_4---, debye_5_e---, debye_5---, debye_6_e---, debye_6+, debye_5_e+, debye_5+, debye_6_e+, debye_6 ) where import Foreign(Ptr)@@ -92,27 +92,27 @@ -- | wrapper for int gsl_sf_debye_5_e(double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_debye_5_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>---debye_5_e :: Double -> (Double,Double)---debye_5_e x = createSFR "debye_5_e" $ gsl_sf_debye_5_e x---foreign import ccall "debye.h gsl_sf_debye_5_e" gsl_sf_debye_5_e :: Double -> Ptr () -> IO CInt+debye_5_e :: Double -> (Double,Double)+debye_5_e x = createSFR "debye_5_e" $ gsl_sf_debye_5_e x+foreign import ccall "debye.h gsl_sf_debye_5_e" gsl_sf_debye_5_e :: Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_debye_5(double x); -- -- <http://www.google.com/search?q=gsl_sf_debye_5&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>---debye_5 :: Double -> Double---debye_5 = gsl_sf_debye_5---foreign import ccall "debye.h gsl_sf_debye_5" gsl_sf_debye_5 :: Double -> Double+debye_5 :: Double -> Double+debye_5 = gsl_sf_debye_5+foreign import ccall "debye.h gsl_sf_debye_5" gsl_sf_debye_5 :: Double -> Double -- | wrapper for int gsl_sf_debye_6_e(double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_debye_6_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>---debye_6_e :: Double -> (Double,Double)---debye_6_e x = createSFR "debye_6_e" $ gsl_sf_debye_6_e x---foreign import ccall "debye.h gsl_sf_debye_6_e" gsl_sf_debye_6_e :: Double -> Ptr () -> IO CInt+debye_6_e :: Double -> (Double,Double)+debye_6_e x = createSFR "debye_6_e" $ gsl_sf_debye_6_e x+foreign import ccall "debye.h gsl_sf_debye_6_e" gsl_sf_debye_6_e :: Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_debye_6(double x); -- -- <http://www.google.com/search?q=gsl_sf_debye_6&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>---debye_6 :: Double -> Double---debye_6 = gsl_sf_debye_6---foreign import ccall "debye.h gsl_sf_debye_6" gsl_sf_debye_6 :: Double -> Double+debye_6 :: Double -> Double+debye_6 = gsl_sf_debye_6+foreign import ccall "debye.h gsl_sf_debye_6" gsl_sf_debye_6 :: Double -> Double
lib/Numeric/GSL/Special/Ellint.hs view
@@ -19,10 +19,10 @@ , ellint_Kcomp , ellint_Ecomp_e , ellint_Ecomp---, ellint_Pcomp_e---, ellint_Pcomp---, ellint_Dcomp_e---, ellint_Dcomp+, ellint_Pcomp_e+, ellint_Pcomp+, ellint_Dcomp_e+, ellint_Dcomp , ellint_F_e , ellint_F , ellint_E_e@@ -76,30 +76,30 @@ -- | wrapper for int gsl_sf_ellint_Pcomp_e(double k,double n,gsl_mode_t mode,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_ellint_Pcomp_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>---ellint_Pcomp_e :: Double -> Double -> Precision -> (Double,Double)---ellint_Pcomp_e k n mode = createSFR "ellint_Pcomp_e" $ gsl_sf_ellint_Pcomp_e k n (precCode mode)---foreign import ccall "ellint.h gsl_sf_ellint_Pcomp_e" gsl_sf_ellint_Pcomp_e :: Double -> Double -> Gsl_mode_t -> Ptr () -> IO CInt+ellint_Pcomp_e :: Double -> Double -> Precision -> (Double,Double)+ellint_Pcomp_e k n mode = createSFR "ellint_Pcomp_e" $ gsl_sf_ellint_Pcomp_e k n (precCode mode)+foreign import ccall "ellint.h gsl_sf_ellint_Pcomp_e" gsl_sf_ellint_Pcomp_e :: Double -> Double -> Gsl_mode_t -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_ellint_Pcomp(double k,double n,gsl_mode_t mode); -- -- <http://www.google.com/search?q=gsl_sf_ellint_Pcomp&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>---ellint_Pcomp :: Double -> Double -> Precision -> Double---ellint_Pcomp k n mode = gsl_sf_ellint_Pcomp k n (precCode mode)---foreign import ccall "ellint.h gsl_sf_ellint_Pcomp" gsl_sf_ellint_Pcomp :: Double -> Double -> Gsl_mode_t -> Double+ellint_Pcomp :: Double -> Double -> Precision -> Double+ellint_Pcomp k n mode = gsl_sf_ellint_Pcomp k n (precCode mode)+foreign import ccall "ellint.h gsl_sf_ellint_Pcomp" gsl_sf_ellint_Pcomp :: Double -> Double -> Gsl_mode_t -> Double -- | wrapper for int gsl_sf_ellint_Dcomp_e(double k,gsl_mode_t mode,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_ellint_Dcomp_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>---ellint_Dcomp_e :: Double -> Precision -> (Double,Double)---ellint_Dcomp_e k mode = createSFR "ellint_Dcomp_e" $ gsl_sf_ellint_Dcomp_e k (precCode mode)---foreign import ccall "ellint.h gsl_sf_ellint_Dcomp_e" gsl_sf_ellint_Dcomp_e :: Double -> Gsl_mode_t -> Ptr () -> IO CInt+ellint_Dcomp_e :: Double -> Precision -> (Double,Double)+ellint_Dcomp_e k mode = createSFR "ellint_Dcomp_e" $ gsl_sf_ellint_Dcomp_e k (precCode mode)+foreign import ccall "ellint.h gsl_sf_ellint_Dcomp_e" gsl_sf_ellint_Dcomp_e :: Double -> Gsl_mode_t -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_ellint_Dcomp(double k,gsl_mode_t mode); -- -- <http://www.google.com/search?q=gsl_sf_ellint_Dcomp&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>---ellint_Dcomp :: Double -> Precision -> Double---ellint_Dcomp k mode = gsl_sf_ellint_Dcomp k (precCode mode)---foreign import ccall "ellint.h gsl_sf_ellint_Dcomp" gsl_sf_ellint_Dcomp :: Double -> Gsl_mode_t -> Double+ellint_Dcomp :: Double -> Precision -> Double+ellint_Dcomp k mode = gsl_sf_ellint_Dcomp k (precCode mode)+foreign import ccall "ellint.h gsl_sf_ellint_Dcomp" gsl_sf_ellint_Dcomp :: Double -> Gsl_mode_t -> Double -- | wrapper for int gsl_sf_ellint_F_e(double phi,double k,gsl_mode_t mode,gsl_sf_result* result); --
lib/Numeric/GSL/Special/Legendre.hs view
@@ -15,24 +15,42 @@ ------------------------------------------------------------ module Numeric.GSL.Special.Legendre(- legendre_Pl+ legendre_Pl_e+, legendre_Pl+, legendre_P1_e+, legendre_P2_e+, legendre_P3_e , legendre_P1 , legendre_P2 , legendre_P3+, legendre_Q0_e , legendre_Q0+, legendre_Q1_e , legendre_Q1+, legendre_Ql_e , legendre_Ql+, legendre_Plm_e , legendre_Plm+, legendre_sphPlm_e , legendre_sphPlm , legendre_array_size+, conicalP_half_e , conicalP_half+, conicalP_mhalf_e , conicalP_mhalf+, conicalP_0_e , conicalP_0+, conicalP_1_e , conicalP_1+, conicalP_sph_reg_e , conicalP_sph_reg+, conicalP_cyl_reg_e , conicalP_cyl_reg+, legendre_H3d_0_e , legendre_H3d_0+, legendre_H3d_1_e , legendre_H3d_1+, legendre_H3d_e , legendre_H3d ) where @@ -40,12 +58,12 @@ import Foreign.C.Types(CInt) import Numeric.GSL.Special.Internal --- | wrapper for int gsl_sf_legendre_Pl_e(int l,double x,double* result);+-- | wrapper for int gsl_sf_legendre_Pl_e(int l,double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_Pl_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_Pl_e :: CInt -> Double -> Ptr Double -> CInt-legendre_Pl_e = gsl_sf_legendre_Pl_e-foreign import ccall "legendre.h gsl_sf_legendre_Pl_e" gsl_sf_legendre_Pl_e :: CInt -> Double -> Ptr Double -> CInt+legendre_Pl_e :: CInt -> Double -> (Double,Double)+legendre_Pl_e l x = createSFR "legendre_Pl_e" $ gsl_sf_legendre_Pl_e l x+foreign import ccall "legendre.h gsl_sf_legendre_Pl_e" gsl_sf_legendre_Pl_e :: CInt -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_legendre_Pl(int l,double x); --@@ -68,26 +86,26 @@ legendre_Pl_deriv_array = gsl_sf_legendre_Pl_deriv_array foreign import ccall "legendre.h gsl_sf_legendre_Pl_deriv_array" gsl_sf_legendre_Pl_deriv_array :: CInt -> Double -> Ptr Double -> Ptr Double -> CInt --- | wrapper for int gsl_sf_legendre_P1_e(double x,double* result);+-- | wrapper for int gsl_sf_legendre_P1_e(double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_P1_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_P1_e :: Double -> Ptr Double -> CInt-legendre_P1_e = gsl_sf_legendre_P1_e-foreign import ccall "legendre.h gsl_sf_legendre_P1_e" gsl_sf_legendre_P1_e :: Double -> Ptr Double -> CInt+legendre_P1_e :: Double -> (Double,Double)+legendre_P1_e x = createSFR "legendre_P1_e" $ gsl_sf_legendre_P1_e x+foreign import ccall "legendre.h gsl_sf_legendre_P1_e" gsl_sf_legendre_P1_e :: Double -> Ptr () -> IO CInt --- | wrapper for int gsl_sf_legendre_P2_e(double x,double* result);+-- | wrapper for int gsl_sf_legendre_P2_e(double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_P2_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_P2_e :: Double -> Ptr Double -> CInt-legendre_P2_e = gsl_sf_legendre_P2_e-foreign import ccall "legendre.h gsl_sf_legendre_P2_e" gsl_sf_legendre_P2_e :: Double -> Ptr Double -> CInt+legendre_P2_e :: Double -> (Double,Double)+legendre_P2_e x = createSFR "legendre_P2_e" $ gsl_sf_legendre_P2_e x+foreign import ccall "legendre.h gsl_sf_legendre_P2_e" gsl_sf_legendre_P2_e :: Double -> Ptr () -> IO CInt --- | wrapper for int gsl_sf_legendre_P3_e(double x,double* result);+-- | wrapper for int gsl_sf_legendre_P3_e(double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_P3_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_P3_e :: Double -> Ptr Double -> CInt-legendre_P3_e = gsl_sf_legendre_P3_e-foreign import ccall "legendre.h gsl_sf_legendre_P3_e" gsl_sf_legendre_P3_e :: Double -> Ptr Double -> CInt+legendre_P3_e :: Double -> (Double,Double)+legendre_P3_e x = createSFR "legendre_P3_e" $ gsl_sf_legendre_P3_e x+foreign import ccall "legendre.h gsl_sf_legendre_P3_e" gsl_sf_legendre_P3_e :: Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_legendre_P1(double x); --@@ -110,12 +128,12 @@ legendre_P3 = gsl_sf_legendre_P3 foreign import ccall "legendre.h gsl_sf_legendre_P3" gsl_sf_legendre_P3 :: Double -> Double --- | wrapper for int gsl_sf_legendre_Q0_e(double x,double* result);+-- | wrapper for int gsl_sf_legendre_Q0_e(double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_Q0_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_Q0_e :: Double -> Ptr Double -> CInt-legendre_Q0_e = gsl_sf_legendre_Q0_e-foreign import ccall "legendre.h gsl_sf_legendre_Q0_e" gsl_sf_legendre_Q0_e :: Double -> Ptr Double -> CInt+legendre_Q0_e :: Double -> (Double,Double)+legendre_Q0_e x = createSFR "legendre_Q0_e" $ gsl_sf_legendre_Q0_e x+foreign import ccall "legendre.h gsl_sf_legendre_Q0_e" gsl_sf_legendre_Q0_e :: Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_legendre_Q0(double x); --@@ -124,12 +142,12 @@ legendre_Q0 = gsl_sf_legendre_Q0 foreign import ccall "legendre.h gsl_sf_legendre_Q0" gsl_sf_legendre_Q0 :: Double -> Double --- | wrapper for int gsl_sf_legendre_Q1_e(double x,double* result);+-- | wrapper for int gsl_sf_legendre_Q1_e(double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_Q1_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_Q1_e :: Double -> Ptr Double -> CInt-legendre_Q1_e = gsl_sf_legendre_Q1_e-foreign import ccall "legendre.h gsl_sf_legendre_Q1_e" gsl_sf_legendre_Q1_e :: Double -> Ptr Double -> CInt+legendre_Q1_e :: Double -> (Double,Double)+legendre_Q1_e x = createSFR "legendre_Q1_e" $ gsl_sf_legendre_Q1_e x+foreign import ccall "legendre.h gsl_sf_legendre_Q1_e" gsl_sf_legendre_Q1_e :: Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_legendre_Q1(double x); --@@ -138,12 +156,12 @@ legendre_Q1 = gsl_sf_legendre_Q1 foreign import ccall "legendre.h gsl_sf_legendre_Q1" gsl_sf_legendre_Q1 :: Double -> Double --- | wrapper for int gsl_sf_legendre_Ql_e(int l,double x,double* result);+-- | wrapper for int gsl_sf_legendre_Ql_e(int l,double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_Ql_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_Ql_e :: CInt -> Double -> Ptr Double -> CInt-legendre_Ql_e = gsl_sf_legendre_Ql_e-foreign import ccall "legendre.h gsl_sf_legendre_Ql_e" gsl_sf_legendre_Ql_e :: CInt -> Double -> Ptr Double -> CInt+legendre_Ql_e :: CInt -> Double -> (Double,Double)+legendre_Ql_e l x = createSFR "legendre_Ql_e" $ gsl_sf_legendre_Ql_e l x+foreign import ccall "legendre.h gsl_sf_legendre_Ql_e" gsl_sf_legendre_Ql_e :: CInt -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_legendre_Ql(int l,double x); --@@ -152,12 +170,12 @@ legendre_Ql = gsl_sf_legendre_Ql foreign import ccall "legendre.h gsl_sf_legendre_Ql" gsl_sf_legendre_Ql :: CInt -> Double -> Double --- | wrapper for int gsl_sf_legendre_Plm_e(int l,int m,double x,double* result);+-- | wrapper for int gsl_sf_legendre_Plm_e(int l,int m,double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_Plm_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_Plm_e :: CInt -> CInt -> Double -> Ptr Double -> CInt-legendre_Plm_e = gsl_sf_legendre_Plm_e-foreign import ccall "legendre.h gsl_sf_legendre_Plm_e" gsl_sf_legendre_Plm_e :: CInt -> CInt -> Double -> Ptr Double -> CInt+legendre_Plm_e :: CInt -> CInt -> Double -> (Double,Double)+legendre_Plm_e l m x = createSFR "legendre_Plm_e" $ gsl_sf_legendre_Plm_e l m x+foreign import ccall "legendre.h gsl_sf_legendre_Plm_e" gsl_sf_legendre_Plm_e :: CInt -> CInt -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_legendre_Plm(int l,int m,double x); --@@ -180,12 +198,12 @@ legendre_Plm_deriv_array = gsl_sf_legendre_Plm_deriv_array foreign import ccall "legendre.h gsl_sf_legendre_Plm_deriv_array" gsl_sf_legendre_Plm_deriv_array :: CInt -> CInt -> Double -> Ptr Double -> Ptr Double -> CInt --- | wrapper for int gsl_sf_legendre_sphPlm_e(int l,int m,double x,double* result);+-- | wrapper for int gsl_sf_legendre_sphPlm_e(int l,int m,double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_sphPlm_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_sphPlm_e :: CInt -> CInt -> Double -> Ptr Double -> CInt-legendre_sphPlm_e = gsl_sf_legendre_sphPlm_e-foreign import ccall "legendre.h gsl_sf_legendre_sphPlm_e" gsl_sf_legendre_sphPlm_e :: CInt -> CInt -> Double -> Ptr Double -> CInt+legendre_sphPlm_e :: CInt -> CInt -> Double -> (Double,Double)+legendre_sphPlm_e l m x = createSFR "legendre_sphPlm_e" $ gsl_sf_legendre_sphPlm_e l m x+foreign import ccall "legendre.h gsl_sf_legendre_sphPlm_e" gsl_sf_legendre_sphPlm_e :: CInt -> CInt -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_legendre_sphPlm(int l,int m,double x); --@@ -215,12 +233,12 @@ legendre_array_size = gsl_sf_legendre_array_size foreign import ccall "legendre.h gsl_sf_legendre_array_size" gsl_sf_legendre_array_size :: CInt -> CInt -> CInt --- | wrapper for int gsl_sf_conicalP_half_e(double lambda,double x,double* result);+-- | wrapper for int gsl_sf_conicalP_half_e(double lambda,double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_conicalP_half_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-conicalP_half_e :: Double -> Double -> Ptr Double -> CInt-conicalP_half_e = gsl_sf_conicalP_half_e-foreign import ccall "legendre.h gsl_sf_conicalP_half_e" gsl_sf_conicalP_half_e :: Double -> Double -> Ptr Double -> CInt+conicalP_half_e :: Double -> Double -> (Double,Double)+conicalP_half_e lambda x = createSFR "conicalP_half_e" $ gsl_sf_conicalP_half_e lambda x+foreign import ccall "legendre.h gsl_sf_conicalP_half_e" gsl_sf_conicalP_half_e :: Double -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_conicalP_half(double lambda,double x); --@@ -229,12 +247,12 @@ conicalP_half = gsl_sf_conicalP_half foreign import ccall "legendre.h gsl_sf_conicalP_half" gsl_sf_conicalP_half :: Double -> Double -> Double --- | wrapper for int gsl_sf_conicalP_mhalf_e(double lambda,double x,double* result);+-- | wrapper for int gsl_sf_conicalP_mhalf_e(double lambda,double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_conicalP_mhalf_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-conicalP_mhalf_e :: Double -> Double -> Ptr Double -> CInt-conicalP_mhalf_e = gsl_sf_conicalP_mhalf_e-foreign import ccall "legendre.h gsl_sf_conicalP_mhalf_e" gsl_sf_conicalP_mhalf_e :: Double -> Double -> Ptr Double -> CInt+conicalP_mhalf_e :: Double -> Double -> (Double,Double)+conicalP_mhalf_e lambda x = createSFR "conicalP_mhalf_e" $ gsl_sf_conicalP_mhalf_e lambda x+foreign import ccall "legendre.h gsl_sf_conicalP_mhalf_e" gsl_sf_conicalP_mhalf_e :: Double -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_conicalP_mhalf(double lambda,double x); --@@ -243,12 +261,12 @@ conicalP_mhalf = gsl_sf_conicalP_mhalf foreign import ccall "legendre.h gsl_sf_conicalP_mhalf" gsl_sf_conicalP_mhalf :: Double -> Double -> Double --- | wrapper for int gsl_sf_conicalP_0_e(double lambda,double x,double* result);+-- | wrapper for int gsl_sf_conicalP_0_e(double lambda,double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_conicalP_0_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-conicalP_0_e :: Double -> Double -> Ptr Double -> CInt-conicalP_0_e = gsl_sf_conicalP_0_e-foreign import ccall "legendre.h gsl_sf_conicalP_0_e" gsl_sf_conicalP_0_e :: Double -> Double -> Ptr Double -> CInt+conicalP_0_e :: Double -> Double -> (Double,Double)+conicalP_0_e lambda x = createSFR "conicalP_0_e" $ gsl_sf_conicalP_0_e lambda x+foreign import ccall "legendre.h gsl_sf_conicalP_0_e" gsl_sf_conicalP_0_e :: Double -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_conicalP_0(double lambda,double x); --@@ -257,12 +275,12 @@ conicalP_0 = gsl_sf_conicalP_0 foreign import ccall "legendre.h gsl_sf_conicalP_0" gsl_sf_conicalP_0 :: Double -> Double -> Double --- | wrapper for int gsl_sf_conicalP_1_e(double lambda,double x,double* result);+-- | wrapper for int gsl_sf_conicalP_1_e(double lambda,double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_conicalP_1_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-conicalP_1_e :: Double -> Double -> Ptr Double -> CInt-conicalP_1_e = gsl_sf_conicalP_1_e-foreign import ccall "legendre.h gsl_sf_conicalP_1_e" gsl_sf_conicalP_1_e :: Double -> Double -> Ptr Double -> CInt+conicalP_1_e :: Double -> Double -> (Double,Double)+conicalP_1_e lambda x = createSFR "conicalP_1_e" $ gsl_sf_conicalP_1_e lambda x+foreign import ccall "legendre.h gsl_sf_conicalP_1_e" gsl_sf_conicalP_1_e :: Double -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_conicalP_1(double lambda,double x); --@@ -271,12 +289,12 @@ conicalP_1 = gsl_sf_conicalP_1 foreign import ccall "legendre.h gsl_sf_conicalP_1" gsl_sf_conicalP_1 :: Double -> Double -> Double --- | wrapper for int gsl_sf_conicalP_sph_reg_e(int l,double lambda,double x,double* result);+-- | wrapper for int gsl_sf_conicalP_sph_reg_e(int l,double lambda,double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_conicalP_sph_reg_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-conicalP_sph_reg_e :: CInt -> Double -> Double -> Ptr Double -> CInt-conicalP_sph_reg_e = gsl_sf_conicalP_sph_reg_e-foreign import ccall "legendre.h gsl_sf_conicalP_sph_reg_e" gsl_sf_conicalP_sph_reg_e :: CInt -> Double -> Double -> Ptr Double -> CInt+conicalP_sph_reg_e :: CInt -> Double -> Double -> (Double,Double)+conicalP_sph_reg_e l lambda x = createSFR "conicalP_sph_reg_e" $ gsl_sf_conicalP_sph_reg_e l lambda x+foreign import ccall "legendre.h gsl_sf_conicalP_sph_reg_e" gsl_sf_conicalP_sph_reg_e :: CInt -> Double -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_conicalP_sph_reg(int l,double lambda,double x); --@@ -285,12 +303,12 @@ conicalP_sph_reg = gsl_sf_conicalP_sph_reg foreign import ccall "legendre.h gsl_sf_conicalP_sph_reg" gsl_sf_conicalP_sph_reg :: CInt -> Double -> Double -> Double --- | wrapper for int gsl_sf_conicalP_cyl_reg_e(int m,double lambda,double x,double* result);+-- | wrapper for int gsl_sf_conicalP_cyl_reg_e(int m,double lambda,double x,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_conicalP_cyl_reg_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-conicalP_cyl_reg_e :: CInt -> Double -> Double -> Ptr Double -> CInt-conicalP_cyl_reg_e = gsl_sf_conicalP_cyl_reg_e-foreign import ccall "legendre.h gsl_sf_conicalP_cyl_reg_e" gsl_sf_conicalP_cyl_reg_e :: CInt -> Double -> Double -> Ptr Double -> CInt+conicalP_cyl_reg_e :: CInt -> Double -> Double -> (Double,Double)+conicalP_cyl_reg_e m lambda x = createSFR "conicalP_cyl_reg_e" $ gsl_sf_conicalP_cyl_reg_e m lambda x+foreign import ccall "legendre.h gsl_sf_conicalP_cyl_reg_e" gsl_sf_conicalP_cyl_reg_e :: CInt -> Double -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_conicalP_cyl_reg(int m,double lambda,double x); --@@ -299,12 +317,12 @@ conicalP_cyl_reg = gsl_sf_conicalP_cyl_reg foreign import ccall "legendre.h gsl_sf_conicalP_cyl_reg" gsl_sf_conicalP_cyl_reg :: CInt -> Double -> Double -> Double --- | wrapper for int gsl_sf_legendre_H3d_0_e(double lambda,double eta,double* result);+-- | wrapper for int gsl_sf_legendre_H3d_0_e(double lambda,double eta,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_H3d_0_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_H3d_0_e :: Double -> Double -> Ptr Double -> CInt-legendre_H3d_0_e = gsl_sf_legendre_H3d_0_e-foreign import ccall "legendre.h gsl_sf_legendre_H3d_0_e" gsl_sf_legendre_H3d_0_e :: Double -> Double -> Ptr Double -> CInt+legendre_H3d_0_e :: Double -> Double -> (Double,Double)+legendre_H3d_0_e lambda eta = createSFR "legendre_H3d_0_e" $ gsl_sf_legendre_H3d_0_e lambda eta+foreign import ccall "legendre.h gsl_sf_legendre_H3d_0_e" gsl_sf_legendre_H3d_0_e :: Double -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_legendre_H3d_0(double lambda,double eta); --@@ -313,12 +331,12 @@ legendre_H3d_0 = gsl_sf_legendre_H3d_0 foreign import ccall "legendre.h gsl_sf_legendre_H3d_0" gsl_sf_legendre_H3d_0 :: Double -> Double -> Double --- | wrapper for int gsl_sf_legendre_H3d_1_e(double lambda,double eta,double* result);+-- | wrapper for int gsl_sf_legendre_H3d_1_e(double lambda,double eta,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_H3d_1_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_H3d_1_e :: Double -> Double -> Ptr Double -> CInt-legendre_H3d_1_e = gsl_sf_legendre_H3d_1_e-foreign import ccall "legendre.h gsl_sf_legendre_H3d_1_e" gsl_sf_legendre_H3d_1_e :: Double -> Double -> Ptr Double -> CInt+legendre_H3d_1_e :: Double -> Double -> (Double,Double)+legendre_H3d_1_e lambda eta = createSFR "legendre_H3d_1_e" $ gsl_sf_legendre_H3d_1_e lambda eta+foreign import ccall "legendre.h gsl_sf_legendre_H3d_1_e" gsl_sf_legendre_H3d_1_e :: Double -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_legendre_H3d_1(double lambda,double eta); --@@ -327,12 +345,12 @@ legendre_H3d_1 = gsl_sf_legendre_H3d_1 foreign import ccall "legendre.h gsl_sf_legendre_H3d_1" gsl_sf_legendre_H3d_1 :: Double -> Double -> Double --- | wrapper for int gsl_sf_legendre_H3d_e(int l,double lambda,double eta,double* result);+-- | wrapper for int gsl_sf_legendre_H3d_e(int l,double lambda,double eta,gsl_sf_result* result); -- -- <http://www.google.com/search?q=gsl_sf_legendre_H3d_e&as_sitesearch=www.gnu.org/software/gsl/manual&btnI=Lucky>-legendre_H3d_e :: CInt -> Double -> Double -> Ptr Double -> CInt-legendre_H3d_e = gsl_sf_legendre_H3d_e-foreign import ccall "legendre.h gsl_sf_legendre_H3d_e" gsl_sf_legendre_H3d_e :: CInt -> Double -> Double -> Ptr Double -> CInt+legendre_H3d_e :: CInt -> Double -> Double -> (Double,Double)+legendre_H3d_e l lambda eta = createSFR "legendre_H3d_e" $ gsl_sf_legendre_H3d_e l lambda eta+foreign import ccall "legendre.h gsl_sf_legendre_H3d_e" gsl_sf_legendre_H3d_e :: CInt -> Double -> Double -> Ptr () -> IO CInt -- | wrapper for double gsl_sf_legendre_H3d(int l,double lambda,double eta); --
lib/Numeric/GSL/Special/auto.hs view
@@ -1,3 +1,5 @@+#!/usr/bin/env runhaskell+ -- automatic generation of wrappers for simple GSL special functions import Text.ParserCombinators.Parsec
lib/Numeric/GSL/Special/autoall.sh view
@@ -1,41 +1,40 @@ #!/bin/bash -#Airy.hs-# include Precision (..) in the export list-#Exp.hs-# remove extern inline definition, qualify name-#Coupling-# remove deprecated INCORRECT-#Trig.hs-# qualify names-#Legendre.hs-# remove extern inline-#Log.hs-# remove extern inline, qualify name+function rep {+ ./replace.hs "$1" "$2" < $3 > /tmp/tmp-rep+ cp /tmp/tmp-rep $3+} -#runhaskell auto airy-runhaskell auto bessel-runhaskell auto clausen-runhaskell auto coulomb-#runhaskell auto coupling-runhaskell auto dawson-runhaskell auto debye-runhaskell auto dilog-runhaskell auto elementary-runhaskell auto ellint-runhaskell auto erf-#runhaskell auto exp-runhaskell auto expint-runhaskell auto fermi_dirac-runhaskell auto gamma-runhaskell auto gegenbauer-runhaskell auto hyperg-runhaskell auto laguerre-runhaskell auto lambert-#runhaskell auto legendre legendre.h-#runhaskell auto log-runhaskell auto pow_int-runhaskell auto psi-runhaskell auto synchrotron-#runhaskell auto trig-runhaskell auto zeta++./auto.hs airy+rep ') where' ', Precision(..)\n) where' Airy.hs+./auto.hs bessel+./auto.hs clausen+./auto.hs coulomb+./auto.hs coupling+rep ', coupling_6j_INCORRECT_e\n, coupling_6j_INCORRECT\n' '' Coupling.hs+./auto.hs dawson+./auto.hs debye+./auto.hs dilog+./auto.hs elementary+./auto.hs ellint+./auto.hs erf+./auto.hs exp+rep ', exp\n' ', Numeric.GSL.Special.Exp.exp\n' Exp.hs+./auto.hs expint+./auto.hs fermi_dirac+./auto.hs gamma+./auto.hs gegenbauer+./auto.hs hyperg+./auto.hs laguerre+./auto.hs lambert+./auto.hs legendre gsl_sf_legendre.h+./auto.hs log+rep ', log\n' ', Numeric.GSL.Special.Log.log\n' Log.hs+./auto.hs pow_int+./auto.hs psi+./auto.hs synchrotron+./auto.hs trig+rep ', sin\n' ', Numeric.GSL.Special.Trig.sin\n' Trig.hs+rep ', cos\n' ', Numeric.GSL.Special.Trig.cos\n' Trig.hs+./auto.hs zeta
− lib/Numeric/GSL/Special/gsl_sf_exp.h
@@ -1,146 +0,0 @@-/* specfunc/gsl_sf_exp.h- * - * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman- * - * This program is free software; you can redistribute it and/or modify- * it under the terms of the GNU General Public License as published by- * the Free Software Foundation; either version 2 of the License, or (at- * your option) any later version.- * - * This program is distributed in the hope that it will be useful, but- * WITHOUT ANY WARRANTY; without even the implied warranty of- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU- * General Public License for more details.- * - * You should have received a copy of the GNU General Public License- * along with this program; if not, write to the Free Software- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.- */--/* Author: G. Jungman */--#ifndef __GSL_SF_EXP_H__-#define __GSL_SF_EXP_H__--#include <gsl/gsl_sf_result.h>-#include <gsl/gsl_precision.h>--#undef __BEGIN_DECLS-#undef __END_DECLS-#ifdef __cplusplus-# define __BEGIN_DECLS extern "C" {-# define __END_DECLS }-#else-# define __BEGIN_DECLS /* empty */-# define __END_DECLS /* empty */-#endif--__BEGIN_DECLS--/* Provide an exp() function with GSL semantics,- * i.e. with proper error checking, etc.- *- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW- */-int gsl_sf_exp_e(const double x, gsl_sf_result * result);-double gsl_sf_exp(const double x);---/* Exp(x)- *- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW- */-int gsl_sf_exp_e10_e(const double x, gsl_sf_result_e10 * result);---/* Exponentiate and multiply by a given factor: y * Exp(x)- *- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW- */-int gsl_sf_exp_mult_e(const double x, const double y, gsl_sf_result * result);-double gsl_sf_exp_mult(const double x, const double y);---/* Exponentiate and multiply by a given factor: y * Exp(x)- *- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW- */-int gsl_sf_exp_mult_e10_e(const double x, const double y, gsl_sf_result_e10 * result);---/* exp(x)-1- *- * exceptions: GSL_EOVRFLW- */-int gsl_sf_expm1_e(const double x, gsl_sf_result * result);-double gsl_sf_expm1(const double x);---/* (exp(x)-1)/x = 1 + x/2 + x^2/(2*3) + x^3/(2*3*4) + ...- *- * exceptions: GSL_EOVRFLW- */-int gsl_sf_exprel_e(const double x, gsl_sf_result * result);-double gsl_sf_exprel(const double x);---/* 2(exp(x)-1-x)/x^2 = 1 + x/3 + x^2/(3*4) + x^3/(3*4*5) + ...- *- * exceptions: GSL_EOVRFLW- */-int gsl_sf_exprel_2_e(double x, gsl_sf_result * result);-double gsl_sf_exprel_2(const double x);---/* Similarly for the N-th generalization of- * the above. The so-called N-relative exponential- *- * exprel_N(x) = N!/x^N (exp(x) - Sum[x^k/k!, {k,0,N-1}])- * = 1 + x/(N+1) + x^2/((N+1)(N+2)) + ...- * = 1F1(1,1+N,x)- */-int gsl_sf_exprel_n_e(const int n, const double x, gsl_sf_result * result);-double gsl_sf_exprel_n(const int n, const double x);---/* Exponentiate a quantity with an associated error.- */-int gsl_sf_exp_err_e(const double x, const double dx, gsl_sf_result * result);--/* Exponentiate a quantity with an associated error.- */-int gsl_sf_exp_err_e10_e(const double x, const double dx, gsl_sf_result_e10 * result);---/* Exponentiate and multiply by a given factor: y * Exp(x),- * for quantities with associated errors.- *- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW- */-int gsl_sf_exp_mult_err_e(const double x, const double dx, const double y, const double dy, gsl_sf_result * result);---/* Exponentiate and multiply by a given factor: y * Exp(x),- * for quantities with associated errors.- *- * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW- */-int gsl_sf_exp_mult_err_e10_e(const double x, const double dx, const double y, const double dy, gsl_sf_result_e10 * result);--__END_DECLS---#ifdef HAVE_INLINE-#include <gsl/gsl_math.h>-#include <gsl/gsl_errno.h>--__BEGIN_DECLS----__END_DECLS--#endif /* HAVE_INLINE */---#endif /* __GSL_SF_EXP_H__ */
+ lib/Numeric/GSL/Special/gsl_sf_legendre.h view
@@ -0,0 +1,319 @@+/* specfunc/gsl_sf_legendre.h+ * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2004 Gerard Jungman+ * + * This program is free software; you can redistribute it and/or modify+ * it under the terms of the GNU General Public License as published by+ * the Free Software Foundation; either version 2 of the License, or (at+ * your option) any later version.+ * + * This program is distributed in the hope that it will be useful, but+ * WITHOUT ANY WARRANTY; without even the implied warranty of+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU+ * General Public License for more details.+ * + * You should have received a copy of the GNU General Public License+ * along with this program; if not, write to the Free Software+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.+ */++/* Author: G. Jungman */++#ifndef __GSL_SF_LEGENDRE_H__+#define __GSL_SF_LEGENDRE_H__++#include <gsl/gsl_sf_result.h>++#undef __BEGIN_DECLS+#undef __END_DECLS+#ifdef __cplusplus+# define __BEGIN_DECLS extern "C" {+# define __END_DECLS }+#else+# define __BEGIN_DECLS /* empty */+# define __END_DECLS /* empty */+#endif++__BEGIN_DECLS+++/* P_l(x) l >= 0; |x| <= 1+ *+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_Pl_e(const int l, const double x, gsl_sf_result * result);+double gsl_sf_legendre_Pl(const int l, const double x);+++/* P_l(x) for l=0,...,lmax; |x| <= 1+ *+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_Pl_array(+ const int lmax, const double x,+ double * result_array+ );+++/* P_l(x) and P_l'(x) for l=0,...,lmax; |x| <= 1+ *+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_Pl_deriv_array(+ const int lmax, const double x,+ double * result_array,+ double * result_deriv_array+ );+++/* P_l(x), l=1,2,3+ *+ * exceptions: none+ */+int gsl_sf_legendre_P1_e(double x, gsl_sf_result * result);+int gsl_sf_legendre_P2_e(double x, gsl_sf_result * result);+int gsl_sf_legendre_P3_e(double x, gsl_sf_result * result);+double gsl_sf_legendre_P1(const double x);+double gsl_sf_legendre_P2(const double x);+double gsl_sf_legendre_P3(const double x);+++/* Q_0(x), x > -1, x != 1+ *+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_Q0_e(const double x, gsl_sf_result * result);+double gsl_sf_legendre_Q0(const double x);+++/* Q_1(x), x > -1, x != 1+ *+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_Q1_e(const double x, gsl_sf_result * result);+double gsl_sf_legendre_Q1(const double x);+++/* Q_l(x), x > -1, x != 1, l >= 0+ *+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_Ql_e(const int l, const double x, gsl_sf_result * result);+double gsl_sf_legendre_Ql(const int l, const double x);+++/* P_l^m(x) m >= 0; l >= m; |x| <= 1.0+ *+ * Note that this function grows combinatorially with l.+ * Therefore we can easily generate an overflow for l larger+ * than about 150.+ *+ * There is no trouble for small m, but when m and l are both large,+ * then there will be trouble. Rather than allow overflows, these+ * functions refuse to calculate when they can sense that l and m are+ * too big.+ *+ * If you really want to calculate a spherical harmonic, then DO NOT+ * use this. Instead use legendre_sphPlm() below, which uses a similar+ * recursion, but with the normalized functions.+ *+ * exceptions: GSL_EDOM, GSL_EOVRFLW+ */+int gsl_sf_legendre_Plm_e(const int l, const int m, const double x, gsl_sf_result * result);+double gsl_sf_legendre_Plm(const int l, const int m, const double x);+++/* P_l^m(x) m >= 0; l >= m; |x| <= 1.0+ * l=|m|,...,lmax+ *+ * exceptions: GSL_EDOM, GSL_EOVRFLW+ */+int gsl_sf_legendre_Plm_array(+ const int lmax, const int m, const double x,+ double * result_array+ );+++/* P_l^m(x) and d(P_l^m(x))/dx; m >= 0; lmax >= m; |x| <= 1.0+ * l=|m|,...,lmax+ *+ * exceptions: GSL_EDOM, GSL_EOVRFLW+ */+int gsl_sf_legendre_Plm_deriv_array(+ const int lmax, const int m, const double x,+ double * result_array,+ double * result_deriv_array+ );+++/* P_l^m(x), normalized properly for use in spherical harmonics+ * m >= 0; l >= m; |x| <= 1.0+ *+ * There is no overflow problem, as there is for the+ * standard normalization of P_l^m(x).+ *+ * Specifically, it returns:+ *+ * sqrt((2l+1)/(4pi)) sqrt((l-m)!/(l+m)!) P_l^m(x)+ *+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_sphPlm_e(const int l, int m, const double x, gsl_sf_result * result);+double gsl_sf_legendre_sphPlm(const int l, const int m, const double x);+++/* sphPlm(l,m,x) values+ * m >= 0; l >= m; |x| <= 1.0+ * l=|m|,...,lmax+ *+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_sphPlm_array(+ const int lmax, int m, const double x,+ double * result_array+ );+++/* sphPlm(l,m,x) and d(sphPlm(l,m,x))/dx values+ * m >= 0; l >= m; |x| <= 1.0+ * l=|m|,...,lmax+ *+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_sphPlm_deriv_array(+ const int lmax, const int m, const double x,+ double * result_array,+ double * result_deriv_array+ );++++/* size of result_array[] needed for the array versions of Plm+ * (lmax - m + 1)+ */+int gsl_sf_legendre_array_size(const int lmax, const int m);+++/* Irregular Spherical Conical Function+ * P^{1/2}_{-1/2 + I lambda}(x)+ *+ * x > -1.0+ * exceptions: GSL_EDOM+ */+int gsl_sf_conicalP_half_e(const double lambda, const double x, gsl_sf_result * result);+double gsl_sf_conicalP_half(const double lambda, const double x);+++/* Regular Spherical Conical Function+ * P^{-1/2}_{-1/2 + I lambda}(x)+ *+ * x > -1.0+ * exceptions: GSL_EDOM+ */+int gsl_sf_conicalP_mhalf_e(const double lambda, const double x, gsl_sf_result * result);+double gsl_sf_conicalP_mhalf(const double lambda, const double x);+++/* Conical Function+ * P^{0}_{-1/2 + I lambda}(x)+ *+ * x > -1.0+ * exceptions: GSL_EDOM+ */+int gsl_sf_conicalP_0_e(const double lambda, const double x, gsl_sf_result * result);+double gsl_sf_conicalP_0(const double lambda, const double x);+++/* Conical Function+ * P^{1}_{-1/2 + I lambda}(x)+ *+ * x > -1.0+ * exceptions: GSL_EDOM+ */+int gsl_sf_conicalP_1_e(const double lambda, const double x, gsl_sf_result * result);+double gsl_sf_conicalP_1(const double lambda, const double x);+++/* Regular Spherical Conical Function+ * P^{-1/2-l}_{-1/2 + I lambda}(x)+ *+ * x > -1.0, l >= -1+ * exceptions: GSL_EDOM+ */+int gsl_sf_conicalP_sph_reg_e(const int l, const double lambda, const double x, gsl_sf_result * result);+double gsl_sf_conicalP_sph_reg(const int l, const double lambda, const double x);+++/* Regular Cylindrical Conical Function+ * P^{-m}_{-1/2 + I lambda}(x)+ *+ * x > -1.0, m >= -1+ * exceptions: GSL_EDOM+ */+int gsl_sf_conicalP_cyl_reg_e(const int m, const double lambda, const double x, gsl_sf_result * result);+double gsl_sf_conicalP_cyl_reg(const int m, const double lambda, const double x);+++/* The following spherical functions are specializations+ * of Legendre functions which give the regular eigenfunctions+ * of the Laplacian on a 3-dimensional hyperbolic space.+ * Of particular interest is the flat limit, which is+ * Flat-Lim := {lambda->Inf, eta->0, lambda*eta fixed}.+ */+ +/* Zeroth radial eigenfunction of the Laplacian on the+ * 3-dimensional hyperbolic space.+ *+ * legendre_H3d_0(lambda,eta) := sin(lambda*eta)/(lambda*sinh(eta))+ * + * Normalization:+ * Flat-Lim legendre_H3d_0(lambda,eta) = j_0(lambda*eta)+ *+ * eta >= 0.0+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_H3d_0_e(const double lambda, const double eta, gsl_sf_result * result);+double gsl_sf_legendre_H3d_0(const double lambda, const double eta);+++/* First radial eigenfunction of the Laplacian on the+ * 3-dimensional hyperbolic space.+ *+ * legendre_H3d_1(lambda,eta) :=+ * 1/sqrt(lambda^2 + 1) sin(lam eta)/(lam sinh(eta))+ * (coth(eta) - lambda cot(lambda*eta))+ * + * Normalization:+ * Flat-Lim legendre_H3d_1(lambda,eta) = j_1(lambda*eta)+ *+ * eta >= 0.0+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_H3d_1_e(const double lambda, const double eta, gsl_sf_result * result);+double gsl_sf_legendre_H3d_1(const double lambda, const double eta);+++/* l'th radial eigenfunction of the Laplacian on the+ * 3-dimensional hyperbolic space.+ *+ * Normalization:+ * Flat-Lim legendre_H3d_l(l,lambda,eta) = j_l(lambda*eta)+ *+ * eta >= 0.0, l >= 0+ * exceptions: GSL_EDOM+ */+int gsl_sf_legendre_H3d_e(const int l, const double lambda, const double eta, gsl_sf_result * result);+double gsl_sf_legendre_H3d(const int l, const double lambda, const double eta);+++/* Array of H3d(ell), 0 <= ell <= lmax+ */+int gsl_sf_legendre_H3d_array(const int lmax, const double lambda, const double eta, double * result_array);++++++__END_DECLS++#endif /* __GSL_SF_LEGENDRE_H__ */
− lib/Numeric/GSL/Special/gsl_sf_log.h
@@ -1,90 +0,0 @@-/* specfunc/gsl_sf_log.h- * - * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman- * - * This program is free software; you can redistribute it and/or modify- * it under the terms of the GNU General Public License as published by- * the Free Software Foundation; either version 2 of the License, or (at- * your option) any later version.- * - * This program is distributed in the hope that it will be useful, but- * WITHOUT ANY WARRANTY; without even the implied warranty of- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU- * General Public License for more details.- * - * You should have received a copy of the GNU General Public License- * along with this program; if not, write to the Free Software- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.- */--/* Author: G. Jungman */--#ifndef __GSL_SF_LOG_H__-#define __GSL_SF_LOG_H__--#include <gsl/gsl_sf_result.h>--#undef __BEGIN_DECLS-#undef __END_DECLS-#ifdef __cplusplus-# define __BEGIN_DECLS extern "C" {-# define __END_DECLS }-#else-# define __BEGIN_DECLS /* empty */-# define __END_DECLS /* empty */-#endif--__BEGIN_DECLS---/* Provide a logarithm function with GSL semantics.- *- * exceptions: GSL_EDOM- */-int gsl_sf_log_e(const double x, gsl_sf_result * result);-double gsl_sf_log(const double x);---/* Log(|x|)- *- * exceptions: GSL_EDOM- */-int gsl_sf_log_abs_e(const double x, gsl_sf_result * result);-double gsl_sf_log_abs(const double x);---/* Complex Logarithm- * exp(lnr + I theta) = zr + I zi- * Returns argument in [-pi,pi].- *- * exceptions: GSL_EDOM- */-int gsl_sf_complex_log_e(const double zr, const double zi, gsl_sf_result * lnr, gsl_sf_result * theta);---/* Log(1 + x)- *- * exceptions: GSL_EDOM- */-int gsl_sf_log_1plusx_e(const double x, gsl_sf_result * result);-double gsl_sf_log_1plusx(const double x);---/* Log(1 + x) - x- *- * exceptions: GSL_EDOM- */-int gsl_sf_log_1plusx_mx_e(const double x, gsl_sf_result * result);-double gsl_sf_log_1plusx_mx(const double x);---#ifdef HAVE_INLINE-#include <gsl/gsl_math.h>-#include <gsl/gsl_errno.h>--#endif /* HAVE_INLINE */---__END_DECLS--#endif /* __GSL_SF_LOG_H__ */
+ lib/Numeric/GSL/Special/replace.hs view
@@ -0,0 +1,14 @@+#!/usr/bin/env runhaskell++import Data.List(isPrefixOf)+import System(getArgs)++rep (c,r) [] = []+rep (c,r) f@(x:xs) + | c `isPrefixOf` f = r ++ rep (c,r) (drop (length c) f)+ | otherwise = x:(rep (c,r) xs)++main = do+ args <- getArgs+ let [p',r'] = map (rep ("\\n","\n")) args+ interact $ rep (p',r')
lib/Numeric/LinearAlgebra/Algorithms.hs view
@@ -69,7 +69,7 @@ import Data.Array -- | Auxiliary typeclass used to define generic computations for both real and complex matrices.-class (Normed (Matrix t), Linear Matrix t) => Field t where+class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where -- | Singular value decomposition using lapack's dgesvd or zgesvd. svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) luPacked :: Matrix t -> (Matrix t, [Int])@@ -142,7 +142,7 @@ -- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf. ----- If @c = chol m@ then @m == c \<> ctrans c@.+-- If @c = chol m@ then @m == ctrans c \<> c@. chol :: Field t => Matrix t -> Matrix t chol m | m `equal` ctrans m = cholSH m | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
lib/Numeric/LinearAlgebra/Linear.hs view
@@ -1,4 +1,4 @@-{-# OPTIONS_GHC -fglasgow-exts #-}+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-} ----------------------------------------------------------------------------- {- | Module : Numeric.LinearAlgebra.Linear@@ -60,18 +60,7 @@ divide = vectorZipC Div equal u v = dim u == dim v && vectorMax (liftVector magnitude (sub u v)) == 0.0 -instance Linear Matrix Double where- scale x = liftMatrix (scale x)- scaleRecip x = liftMatrix (scaleRecip x)- addConstant x = liftMatrix (addConstant x)- add = liftMatrix2 add- sub = liftMatrix2 sub- mul = liftMatrix2 mul- divide = liftMatrix2 divide- equal a b = cols a == cols b && flatten a `equal` flatten b---instance Linear Matrix (Complex Double) where+instance (Linear Vector a, Container Matrix a) => (Linear Matrix a) where scale x = liftMatrix (scale x) scaleRecip x = liftMatrix (scaleRecip x) addConstant x = liftMatrix (addConstant x)
+ lib/Numeric/LinearAlgebra/Tests.hs view
@@ -0,0 +1,197 @@+-----------------------------------------------------------------------------+{- |+Module : Numeric.LinearAlgebra.Tests+Copyright : (c) Alberto Ruiz 2007+License : GPL-style++Maintainer : Alberto Ruiz (aruiz at um dot es)+Stability : provisional+Portability : portable++Some tests.++-}++module Numeric.LinearAlgebra.Tests(+-- module Numeric.LinearAlgebra.Tests.Instances,+-- module Numeric.LinearAlgebra.Tests.Properties,+ qCheck, runTests+--, runBigTests+) where++import Numeric.LinearAlgebra+import Numeric.LinearAlgebra.Tests.Instances+import Numeric.LinearAlgebra.Tests.Properties+import Test.QuickCheck+import Test.HUnit hiding ((~:),test)+import System.Info+import Data.List(foldl1')+import Numeric.GSL hiding (sin,cos,exp,choose)++qCheck n = check defaultConfig {configSize = const n}++utest str b = TestCase $ assertBool str b++feye n = flipud (ident n) :: Matrix Double++detTest1 = det m == 26+ && det mc == 38 :+ (-3)+ && det (feye 2) == -1+ where+ m = (3><3) + [ 1, 2, 3+ , 4, 5, 7+ , 2, 8, 4 :: Double+ ]+ mc = (3><3)+ [ 1, 2, 3+ , 4, 5, 7+ , 2, 8, i+ ]++--------------------------------------------------------------------++polyEval cs x = foldr (\c ac->ac*x+c) 0 cs++polySolveProp p = length p <2 || last p == 0|| 1E-8 > maximum (map magnitude $ map (polyEval (map (:+0) p)) (polySolve p))++---------------------------------------------------------------------++quad f a b = fst $ integrateQAGS 1E-9 100 f a b++-- A multiple integral can be easily defined using partial application+quad2 f a b g1 g2 = quad h a b+ where h x = quad (f x) (g1 x) (g2 x)++volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) + 0 r (const 0) (\x->sqrt (r*r-x*x))++---------------------------------------------------------------------++besselTest = utest "bessel_J0_e" ( abs (r-expected) < e )+ where (r,e) = bessel_J0_e 5.0+ expected = -0.17759677131433830434739701++exponentialTest = utest "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 )+ where (v,e,err) = exp_e10_e 30.0+ expected = exp 30.0++---------------------------------------------------------------------++nd1 = (3><3) [ 1/2, 1/4, 1/4+ , 0/1, 1/2, 1/4+ , 1/2, 1/4, 1/2 :: Double]++nd2 = (2><2) [1, 0, 1, 1:: Complex Double]++expmTest1 = expm nd1 :~14~: (3><3)+ [ 1.762110887278176+ , 0.478085470590435+ , 0.478085470590435+ , 0.104719410945666+ , 1.709751181805343+ , 0.425725765117601+ , 0.851451530235203+ , 0.530445176063267+ , 1.814470592751009 ]++expmTest2 = expm nd2 :~15~: (2><2)+ [ 2.718281828459045+ , 0.000000000000000+ , 2.718281828459045+ , 2.718281828459045 ]+++---------------------------------------------------------------------++rot :: Double -> Matrix Double+rot a = (3><3) [ c,0,s+ , 0,1,0+ ,-s,0,c ]+ where c = cos a+ s = sin a++rotTest = fun (10^5) :~12~: rot 5E4+ where fun n = foldl1' (<>) (map rot angles)+ where angles = toList $ linspace n (0,1)+++-- | All tests must pass with a maximum dimension of about 20+-- (some tests may fail with bigger sizes due to precision loss).+runTests :: Int -- ^ maximum dimension+ -> IO ()+runTests n = do+ setErrorHandlerOff+ let test p = qCheck n p+ putStrLn "------ lu"+ test (luProp . rM)+ test (luProp . cM)+ putStrLn "------ inv"+ test (invProp . rSqWC)+ test (invProp . cSqWC)+ putStrLn "------ pinv"+ test (pinvProp . rM)+ if os == "mingw32"+ then putStrLn "complex pinvTest skipped in this OS"+ else test (pinvProp . cM)+ putStrLn "------ det"+ test (detProp . rSqWC)+ test (detProp . cSqWC)+ putStrLn "------ svd"+ test (svdProp1 . rM)+ test (svdProp1 . cM)+ test (svdProp2 . rM)+ test (svdProp2 . cM)+ putStrLn "------ eig"+ test (eigSHProp . rHer)+ test (eigSHProp . cHer)+ test (eigProp . rSq)+ test (eigProp . cSq)+ putStrLn "------ nullSpace"+ test (nullspaceProp . rM)+ test (nullspaceProp . cM)+ putStrLn "------ qr"+ test (qrProp . rM)+ test (qrProp . cM)+ putStrLn "------ hess"+ test (hessProp . rSq)+ test (hessProp . cSq)+ putStrLn "------ schur"+ test (schurProp2 . rSq)+ if os == "mingw32"+ then putStrLn "complex schur skipped in this OS"+ else test (schurProp1 . cSq)+ putStrLn "------ chol"+ test (cholProp . rPosDef)+ test (cholProp . cPosDef)+ putStrLn "------ expm"+ test (expmDiagProp . rSqWC)+ test (expmDiagProp . cSqWC)+ putStrLn "------ fft"+ test (\v -> ifft (fft v) |~| v)+ putStrLn "------ vector operations"+ test (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM))+ test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))+ test (\u -> cos u * tan u |~| sin (u::RM))+ test (\u -> (cos u * tan u) |~| sin (u::CM))+ putStrLn "------ some unit tests"+ runTestTT $ TestList+ [ utest "1E5 rots" rotTest+ , utest "det1" detTest1+ , utest "expm1" (expmTest1)+ , utest "expm2" (expmTest2)+ , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)+ , utest "arith2" $ (((1+i) .* ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( (140*i-51).*1 :: CM)+ , utest "arith3" $ exp (i.*ones(10,10)*pi) + 1 |~| 0+ , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3]+ , utest "gamma" (gamma 5 == 24.0)+ , besselTest+ , exponentialTest+ , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < 1E-8)+ , utest "polySolve" (polySolveProp [1,2,3,4])+ ]+ return ()++-- -- | Some additional tests on big matrices. They take a few minutes.+-- runBigTests :: IO ()+-- runBigTests = undefined
+ lib/Numeric/LinearAlgebra/Tests/Instances.hs view
@@ -0,0 +1,142 @@+{-# OPTIONS #-}+-----------------------------------------------------------------------------+{- |+Module : Numeric.LinearAlgebra.Tests.Instances+Copyright : (c) Alberto Ruiz 2008+License : GPL-style++Maintainer : Alberto Ruiz (aruiz at um dot es)+Stability : provisional+Portability : portable++Arbitrary instances for vectors, matrices.++-}++module Numeric.LinearAlgebra.Tests.Instances(+ Sq(..), rSq,cSq,+ Rot(..), rRot,cRot,+ Her(..), rHer,cHer,+ WC(..), rWC,cWC,+ SqWC(..), rSqWC, cSqWC,+ PosDef(..), rPosDef, cPosDef,+ RM,CM, rM,cM+) where++import Numeric.LinearAlgebra+import Test.QuickCheck+import Control.Monad(replicateM)++instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where+ arbitrary = do+ r <- arbitrary+ i <- arbitrary+ return (r:+i)+ coarbitrary = undefined++chooseDim = sized $ \m -> choose (1,max 1 m)++instance (Field a, Arbitrary a) => Arbitrary (Vector a) where + arbitrary = do m <- chooseDim+ l <- vector m+ return $ fromList l+ coarbitrary = undefined++instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where + arbitrary = do+ m <- chooseDim+ n <- chooseDim+ l <- vector (m*n)+ return $ (m><n) l+ coarbitrary = undefined++-- a square matrix+newtype (Sq a) = Sq (Matrix a) deriving Show+instance (Element a, Arbitrary a) => Arbitrary (Sq a) where+ arbitrary = do+ n <- chooseDim+ l <- vector (n*n)+ return $ Sq $ (n><n) l+ coarbitrary = undefined++-- a unitary matrix+newtype (Rot a) = Rot (Matrix a) deriving Show+instance (Field a, Arbitrary a) => Arbitrary (Rot a) where+ arbitrary = do+ Sq m <- arbitrary+ let (q,_) = qr m+ return (Rot q)+ coarbitrary = undefined++-- a complex hermitian or real symmetric matrix+newtype (Her a) = Her (Matrix a) deriving Show+instance (Field a, Arbitrary a) => Arbitrary (Her a) where+ arbitrary = do+ Sq m <- arbitrary+ let m' = m/2+ return $ Her (m' + ctrans m')+ coarbitrary = undefined++-- a well-conditioned general matrix (the singular values are between 1 and 100)+newtype (WC a) = WC (Matrix a) deriving Show+instance (Field a, Arbitrary a) => Arbitrary (WC a) where+ arbitrary = do+ m <- arbitrary+ let (u,_,v) = svd m+ r = rows m+ c = cols m+ n = min r c+ sv <- replicateM n (choose (1,100))+ let s = diagRect (fromList sv) r c+ return $ WC (u <> real s <> trans v)+ coarbitrary = undefined++-- a well-conditioned square matrix (the singular values are between 1 and 100)+newtype (SqWC a) = SqWC (Matrix a) deriving Show+instance (Field a, Arbitrary a) => Arbitrary (SqWC a) where+ arbitrary = do+ Sq m <- arbitrary+ let (u,_,v) = svd m+ n = rows m+ sv <- replicateM n (choose (1,100))+ let s = diag (fromList sv)+ return $ SqWC (u <> real s <> trans v)+ coarbitrary = undefined++-- a positive definite square matrix (the eigenvalues are between 0 and 100)+newtype (PosDef a) = PosDef (Matrix a) deriving Show+instance (Field a, Arbitrary a) => Arbitrary (PosDef a) where+ arbitrary = do+ Her m <- arbitrary+ let (_,v) = eigSH m+ n = rows m+ l <- replicateM n (choose (0,100))+ let s = diag (fromList l)+ p = v <> real s <> ctrans v+ return $ PosDef (0.5 .* p + 0.5 .* ctrans p)+ coarbitrary = undefined++type RM = Matrix Double+type CM = Matrix (Complex Double)++rM m = m :: RM+cM m = m :: CM++rHer (Her m) = m :: RM+cHer (Her m) = m :: CM++rRot (Rot m) = m :: RM+cRot (Rot m) = m :: CM++rSq (Sq m) = m :: RM+cSq (Sq m) = m :: CM++rWC (WC m) = m :: RM+cWC (WC m) = m :: CM++rSqWC (SqWC m) = m :: RM+cSqWC (SqWC m) = m :: CM++rPosDef (PosDef m) = m :: RM+cPosDef (PosDef m) = m :: CM+
+ lib/Numeric/LinearAlgebra/Tests/Properties.hs view
@@ -0,0 +1,149 @@+{-# OPTIONS #-}+-----------------------------------------------------------------------------+{- |+Module : Numeric.LinearAlgebra.Tests.Properties+Copyright : (c) Alberto Ruiz 2008+License : GPL-style++Maintainer : Alberto Ruiz (aruiz at um dot es)+Stability : provisional+Portability : portable++Testing properties.++-}++module Numeric.LinearAlgebra.Tests.Properties (+ dist, (|~|), (~:), Aprox((:~)),+ zeros, ones,+ square,+ unitary,+ hermitian,+ wellCond,+ positiveDefinite,+ upperTriang,+ upperHessenberg,+ luProp,+ invProp,+ pinvProp,+ detProp,+ nullspaceProp,+ svdProp1, svdProp2,+ eigProp, eigSHProp,+ qrProp,+ hessProp,+ schurProp1, schurProp2,+ cholProp,+ expmDiagProp+) where++import Numeric.LinearAlgebra+import Numeric.LinearAlgebra.Tests.Instances(Sq(..),Her(..),Rot(..))+import Test.QuickCheck++-- relative error+dist :: (Normed t, Num t) => t -> t -> Double+dist a b = r+ where norm = pnorm Infinity+ na = norm a+ nb = norm b+ nab = norm (a-b)+ mx = max na nb+ mn = min na nb+ r = if mn < eps+ then mx+ else nab/mx++infixl 4 |~|+a |~| b = a :~10~: b+--a |~| b = dist a b < 10^^(-10)++data Aprox a = (:~) a Int+(~:) :: (Normed a, Num a) => Aprox a -> a -> Bool+a :~n~: b = dist a b < 10^^(-n)++------------------------------------------------------++square m = rows m == cols m++unitary m = square m && m <> ctrans m |~| ident (rows m)++hermitian m = square m && m |~| ctrans m++wellCond m = rcond m > 1/100++positiveDefinite m = minimum (toList e) > 0+ where (e,v) = eigSH m++upperTriang m = rows m == 1 || down == z+ where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))+ z = constant 0 (dim down)++upperHessenberg m = rows m < 3 || down == z+ where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))+ z = constant 0 (dim down)++zeros (r,c) = reshape c (constant 0 (r*c))++ones (r,c) = zeros (r,c) + 1++-----------------------------------------------------++luProp m = m |~| p <> l <> u && det p == s+ where (l,u,p,s) = lu m++invProp m = m <> inv m |~| ident (rows m)++pinvProp m = m <> p <> m |~| m+ && p <> m <> p |~| p+ && hermitian (m<>p)+ && hermitian (p<>m)+ where p = pinv m++detProp m = s d1 |~| s d2+ where d1 = det m+ d2 = det' m * det q+ det' m = product $ toList $ takeDiag r+ (q,r) = qr m+ s x = fromList [x]++nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c))+ where nl = nullspacePrec 1 m+ n = fromColumns nl+ r = rows m+ c = cols m - rank m++svdProp1 m = u <> real d <> trans v |~| m+ && unitary u && unitary v+ where (u,d,v) = full svd m++svdProp2 m = (m |~| 0) `trivial` ((m |~| 0) || u <> real (diag s) <> trans v |~| m)+ where (u,s,v) = economy svd m++eigProp m = complex m <> v |~| v <> diag s+ where (s, v) = eig m++eigSHProp m = m <> v |~| v <> real (diag s)+ && unitary v+ && m |~| v <> real (diag s) <> ctrans v+ where (s, v) = eigSH m++qrProp m = q <> r |~| m && unitary q && upperTriang r+ where (q,r) = qr m++hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h+ where (p,h) = hess m++schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s+ where (u,s) = schur m++schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme+ where (u,s) = schur m++cholProp m = m |~| ctrans c <> c && upperTriang c+ where c = chol m+ pos = positiveDefinite m++expmDiagProp m = expm (logm m) :~ 7 ~: complex m+ where logm m = matFunc log m+