packages feed

numeric-prelude 0.4.0.3 → 0.4.1

raw patch · 15 files changed

+216/−50 lines, 15 filesdep ~array

Dependency ranges changed: array

Files

numeric-prelude.cabal view
@@ -1,5 +1,5 @@ Name:           numeric-prelude-Version:        0.4.0.3+Version:        0.4.1 License:        BSD3 License-File:   LICENSE Author:         Dylan Thurston <dpt@math.harvard.edu>, Henning Thielemann <numericprelude@henning-thielemann.de>, Mikael Johansson@@ -151,7 +151,7 @@   default:     False  Source-Repository this-  Tag:         0.4.0.3+  Tag:         0.4.1   Type:        darcs   Location:    http://code.haskell.org/numeric-prelude/ @@ -170,7 +170,7 @@    -- splitBase   Build-Depends:-    array >=0.1 && <0.5,+    array >=0.1 && <0.6,     containers >=0.1 && <0.6,     random >=1.0 && <1.1,     base >= 2 && <5
src/Algebra/Additive.hs view
@@ -58,6 +58,7 @@ -}  class C a where+    {-# MINIMAL zero, (+), ((-) | negate) #-}     -- | zero element of the vector space     zero     :: a     -- | add and subtract elements
src/Algebra/Algebraic.hs view
@@ -21,6 +21,7 @@ {- | Minimal implementation: 'root' or '(^\/)'. -}  class (Field.C a) => C a where+    {-# MINIMAL root | (^/) #-}     sqrt :: a -> a     sqrt = root 2     -- sqrt x  =  x ** (1/2)
src/Algebra/Field.hs view
@@ -68,6 +68,7 @@ -}  class (Ring.C a) => C a where+    {-# MINIMAL recip | (/) #-}     (/)           :: a -> a -> a     recip         :: a -> a     fromRational' :: Rational -> a
src/Algebra/IntegralDomain.hs view
@@ -95,6 +95,7 @@ Minimal definition: 'divMod' or ('div' and 'mod') -} class (Ring.C a) => C a where+    {-# MINIMAL divMod | (div, mod) #-}     div, mod :: a -> a -> a     divMod :: a -> a -> (a,a)
src/Algebra/RealRing.hs view
@@ -115,6 +115,7 @@ -}  class (Absolute.C a, Ord a) => C a where+    {-# MINIMAL splitFraction | floor #-}     splitFraction    :: (Ring.C b) => a -> (b,a)     fraction         ::               a -> a     ceiling, floor   :: (Ring.C b) => a -> b
src/Algebra/Ring.hs view
@@ -65,6 +65,7 @@ -}  class (Additive.C a) => C a where+    {-# MINIMAL (*), (one | fromInteger) #-}     (*)         :: a -> a -> a     one         :: a     fromInteger :: Integer -> a
src/Algebra/Transcendental.hs view
@@ -31,9 +31,10 @@ branch cuts, etc.  Minimal complete definition:-     pi, exp, log, sin, cos, asin, acos, atan+     pi, exp, (log or logBase), sin, cos, atan -} class (Algebraic.C a) => C a where+    {-# MINIMAL pi, exp, (log | logBase), sin, cos, atan #-}     pi                  :: a     exp, log            :: a -> a     logBase, (**)       :: a -> a -> a@@ -56,6 +57,7 @@      x ** y           =  exp (log x * y)     logBase x y      =  log y / log x+    log              =  logBase (exp 1)      tan  x           =  sin x / cos x
src/Algebra/Units.hs view
@@ -70,6 +70,7 @@ -}  class (Integral.C a) => C a where+  {-# MINIMAL isUnit, (stdUnit | stdUnitInv) #-}   isUnit :: a -> Bool   stdAssociate, stdUnit, stdUnitInv :: a -> a
src/MathObj/Polynomial/Core.hs view
@@ -152,21 +152,26 @@ {- snd $ Poly.divMod (repeat (1::Double)) [1,1] -}+{- |+The modulus will always have one element less than the divisor.+This means that the modulus will be denormalized in some cases,+e.g. @mod [2,1,1] [1,1,1] == [1,0]@ instead of @[1]@.+-} divModRev :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a]) divModRev x y =-   let (y0:ys) = dropWhile isZero y-       -- the second parameter represents lazily (length x - length y)-       aux xs' =-         forcePair .-         switchL-           ([], xs')-           (P.const $-              let (x0:xs) = xs'-                  q0      = x0/y0-              in  mapFst (q0:) . aux (sub xs (scale q0 ys)))-   in  if isZero y-         then error "MathObj.Polynomial: division by zero"-         else aux x (drop (length y - 1) x)+   case dropWhile isZero y of+      [] -> error "MathObj.Polynomial: division by zero"+      y0:ys ->+         let -- the second parameter represents lazily (length x - length (normalize y))+             aux xs' =+               forcePair .+               switchL+                 ([], xs')+                 (P.const $+                    let (x0:xs) = xs'+                        q0      = x0/y0+                    in  mapFst (q0:) . aux (sub xs (scale q0 ys)))+         in  aux x (drop (length ys) x)  {-# INLINE stdUnit #-} stdUnit :: (ZeroTestable.C a, Ring.C a) => [a] -> a
src/MathObj/PowerSeries/Core.hs view
@@ -76,6 +76,16 @@    zipWith id (cycle [id, P.const zero, NP.negate, P.const zero])  +{- |+For power series of @f x@, compute the power series of @f(x^n)@.+-}+insertHoles :: Additive.C a => Int -> [a] -> [a]+insertHoles n =+   if n<=0+     then error $ "insertHoles requires positive exponent, but got " ++ show n+     else concatMap (\x -> x : replicate (n-1) zero)++ {- * Series arithmetic -}  add, sub :: (Additive.C a) => [a] -> [a] -> [a]@@ -159,18 +169,23 @@ {- pow alpha t = t^alpha (pow alpha . x)' = alpha * (pow (alpha-1) . x) * x'-alpha * (pow alpha . x) = x * x' * (pow alpha . x)'+(pow alpha . x)' * x = alpha * (pow alpha . x) * x'+ y = pow alpha . x-alpha * y = x * x' * y'+y' * x = alpha * y * x'++This yields an implementation that is a fused+exp (alpha * log x) -}  {- |-Input series must start with non-zero term.+Input series must start with a non-zero term,+even better with a positive one. -} pow :: (Field.C a) => (a -> a) -> a -> [a] -> [a] pow f0 expon x =    let y  = integrate (f0 (head x)) y'-       y' = scale expon (divide y (mul x (differentiate x)))+       y' = scale expon (mul y (derivedLog x))    in  y  @@ -257,22 +272,56 @@ composeTaylor x []     = x 0  +{-+X(t) = t*x(t)+R(t) = t*r(t) +r(t) = 1 / (x(r(t)*t))+R(t)/t+   = 1 / (x(R(t)))+   = 1 / (X(R(t)) / R(t))+   = 1 / (t / R(t))+-}++{- |+This function returns the series of the inverse function in the form:+(point of the expansion, power series).++That is, say we have the equation:++> y = a + f(x)++where function f is given by a power series with f(0) = 0.+We want to solve for x:++> x = f^-1(y-a)++If you pass the power series of @a+f(x)@ to 'inv',+you get @(a, f^-1)@ as answer, where @f^-1@ is a power series.++The linear term of @f@ (the coefficient of @x@) must be non-zero.++This needs cubic run-time and thus is exceptionally slow.+Computing inverse series for special power series might be faster.+-}+-- how about NonEmpty.T here?+inv :: (Eq a, Field.C a) => [a] -> (a, [a])+inv [] = error "inv: power series must be non-zero"+inv (x:xs) =+   (x, let r = divide [1] (compose xs r) in 0 : r)++ {- (x . y) = id (x' . y) * y' = 1 y' = 1 / (x' . y) -} -{- |-This function returns the series of the function in the form:-(point of the expansion, power series)--This is exceptionally slow and needs cubic run-time.+{-+Like 'inv' but with a slightly cumbersome implementation. -}--inv :: (Field.C a) => [a] -> (a, [a])-inv x =+invDiff :: (Field.C a) => [a] -> (a, [a])+invDiff x =    let y' = divide [1] (compose (differentiate x) (tail y))        y  = integrate 0 y'             -- the first term is zero, which is required for composition
src/MathObj/PowerSeries2.hs view
@@ -86,6 +86,11 @@ const x = lift0 [[x]]  +{-# INLINE truncate #-}+truncate :: Int -> T a -> T a+truncate n = lift1 (take n)++ instance Functor T where    fmap f (Cons xs) = Cons (map (map f) xs) @@ -124,5 +129,4 @@  instance (Algebraic.C a) => Algebraic.C (T a) where    sqrt   = lift1 (Core.sqrt Algebraic.sqrt)---   x ^/ y = lift1 (Core.pow (Algebraic.^/ y)---                       (fromRational' y)) x+   x ^/ y = lift1 (Core.pow (Algebraic.^/ y) (fromRational' y)) x
src/MathObj/PowerSeries2/Core.hs view
@@ -59,6 +59,11 @@    lift1fromPowerSeries $    PSCore.sqrt (PS.const . (\[x] -> fSqRt x) . PS.coeffs) +pow :: (Field.C a) =>+   (a -> a) -> a -> T a -> T a+pow fPow expon =+   lift1fromPowerSeries $+   PSCore.pow (PS.const . (\[x] -> fPow x) . PS.coeffs) (PS.const expon)   swapVariables :: T a -> T a
test/Test/MathObj/Polynomial.hs view
@@ -11,6 +11,7 @@ import qualified Algebra.Laws as Laws  import qualified Data.List as List+import Data.Tuple.HT (mapPair, mapSnd, )  import Test.NumericPrelude.Utility (testUnit) import Test.QuickCheck (Property, quickCheck, (==>), Testable, )@@ -30,7 +31,35 @@ mul :: (Ring.C a, Eq a, ZeroTestable.C a) => [a] -> [a] -> Bool mul xs ys  =  PolyCore.equal (PolyCore.mul xs ys) (PolyCore.mulShear xs ys) +divNormal :: [Rational] -> [Rational] -> Property+divNormal x y =+   case (PolyCore.normalize x, PolyCore.normalize y) of+      (nx, ny) ->+         not (null ny) ==>+            mapSnd PolyCore.normalize (PolyCore.divMod nx ny)+            ==+            mapPair+               (PolyCore.normalize, PolyCore.normalize)+               (PolyCore.divMod x y) +normalizedQuotient :: [Rational] -> [Rational] -> Property+normalizedQuotient x y =+   case PolyCore.normalize x of+      nx ->+         not (isZero y) ==>+            let z = fst $ PolyCore.divMod nx y+            in  PolyCore.normalize z == z++modulusSize :: [Rational] -> [Rational] -> Property+modulusSize x y =+   case PolyCore.normalize y of+      ny ->+         not (null ny) ==>+            List.length (snd $ PolyCore.divMod x y)+            <+            List.length ny++ test :: Testable a => (Poly.T Integer -> a) -> IO () test = quickCheck @@ -52,5 +81,8 @@       ("multiplication, commutative",  test (Laws.commutative (*))) :       ("multiplication, associative",  test (Laws.associative (*))) :       ("multiplication and addition, distributive",   test (Laws.leftDistributive (*) (+))) :-      ("division",       testRat (Integral.propInverse)) :+      ("division",            testRat Integral.propInverse) :+      ("division, normalize", quickCheck divNormal) :+      ("normalized quotient", quickCheck normalizedQuotient) :+      ("modulus size",        quickCheck modulusSize) :       []
test/Test/MathObj/PowerSeries.hs view
@@ -1,12 +1,13 @@ {-# LANGUAGE NoImplicitPrelude #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE FlexibleInstances #-} module Test.MathObj.PowerSeries where +import qualified MathObj.PowerSeries         as PST import qualified MathObj.PowerSeries.Core    as PS import qualified MathObj.PowerSeries.Example as PSE -import Test.NumericPrelude.Utility (equalInfLists {- , testUnit -} )+import qualified Test.QuickCheck.Modifiers as Mod+import Test.NumericPrelude.Utility (equalInfLists, testUnit)+import Test.QuickCheck (quickCheck) -- import Test.QuickCheck (Property, quickCheck, (==>)) import qualified Test.HUnit as HUnit @@ -58,6 +59,20 @@    ("sqrt",  500, 1:1:repeat 0, PS.sqrt (\1 -> 1) (PS.mul [1,1] [1,1])) :    [] +identitiesHoles :: [(String, Int, [Rational] -> [Rational], Rational)]+identitiesHoles =+   ("exp",    30, PS.exp  (\0 -> 1), 0) :+   ("log",    30, PS.log  (\1 -> 0), 1) :+   ("tan",    20, PS.tan  (\0 -> (0,1)), 0) :+   ("atan",   20, PS.atan (\0 -> 0), 0) :+   ("sin",    20, PS.sin  (\0 -> (0,1)), 0) :+   ("cos",    20, PS.cos  (\0 -> (0,1)), 0) :+   ("asin",   30, PS.asin (\1 -> 1) (\0 -> 0), 0) :+   ("sqrt",   50, PS.sqrt (\1 -> 1), 1) :+   ("pow13",  30, PS.pow  (\1 -> 1) (1/3), 1) :+   ("pow25",  30, PS.pow  (\1 -> 1) (2/5), 1) :+   []+ testSeriesIdentity :: (String, Int, [Rational], [Rational]) -> HUnit.Test testSeriesIdentity (label, len, x, y) =    HUnit.test (HUnit.assertBool label (equalInfLists len [x,y]))@@ -68,36 +83,83 @@    HUnit.TestLabel label $      HUnit.TestList $ map testSeriesIdentity ids -checkSeriesIdentities ::+_checkSeriesIdentities ::    [(String, Int, [Rational], [Rational])] -> [(String,Bool)]-checkSeriesIdentities =+_checkSeriesIdentities =    map (\(label, len, x, y) -> (label, equalInfLists len [x,y]))  +holesMultiplicative :: Int -> Int -> Int -> [Rational] -> Bool+holesMultiplicative trunc expon0 expon1 xs =+   let n0 = 1 + mod expon0 10+       n1 = 1 + mod expon1 10+   in  equalInfLists trunc+          [PS.insertHoles n0 $ PS.insertHoles n1 xs,+           PS.insertHoles n1 $ PS.insertHoles n0 xs,+           PS.insertHoles (n0*n1) xs] +testHolesIdentity ::+   (String, Int, [Rational] -> [Rational], Rational) -> HUnit.Test+testHolesIdentity (label, len, f, x0) =+   HUnit.test $ testUnit $ (,) ("holes in " ++ label) $+   quickCheck $ \expon0 xs -> checkHoles len expon0 f x0 xs -powerMult :: Rational -> Rational -> Bool-powerMult exp0 exp1 =-   PS.mul (PSE.pow exp0) (PSE.pow exp1)  ==  PSE.pow (exp0+exp1) -powerExplODE :: Rational -> Bool-powerExplODE expon =-   PSE.powODE expon == PSE.powExpl expon+checkHoles ::+   Int -> Int -> ([Rational] -> [Rational]) ->+   Rational -> [Rational] -> Bool+checkHoles trunc expon0 f x xs =+   let expon = 1 + mod expon0 10+   in  equalInfLists trunc+          [(f $ PS.insertHoles expon (x:xs)) ++ repeat zero,+           (PS.insertHoles expon $ f $ x:xs) ++ repeat zero]  +powerMultSeries :: Int -> Integer -> Mod.Positive Rational -> [Rational] -> Bool+powerMultSeries trunc expon0 xp xs =+   let expon = 1 + mod expon0 10+       x = Mod.getPositive xp+       xt = x:xs+   in  equalInfLists trunc+          [PS.pow+              (const x) (1 % expon)+              (PST.coeffs (PST.fromCoeffs xt ^ expon))+            ++ repeat zero,+           xt ++ repeat zero]++powerMult :: Int -> Rational -> Rational -> Bool+powerMult trunc exp0 exp1 =+   equalInfLists trunc+      [PS.mul (PSE.pow exp0) (PSE.pow exp1), PSE.pow (exp0+exp1)]++powerExplODE :: Int -> Rational -> Bool+powerExplODE trunc expon =+   equalInfLists trunc [PSE.powODE expon, PSE.powExpl expon]++invDiff :: Int -> Rational -> Mod.NonZero Rational -> [Rational] -> Bool+invDiff trunc x0 x1 xs_ =+   let xs = x0 : Mod.getNonZero x1 : xs_+       (y,ys) = PS.inv xs+       (z,zs) = PS.invDiff xs+   in  y==z && equalInfLists trunc [ys, zs]++ tests :: HUnit.Test tests =    HUnit.TestLabel "power series" $    HUnit.TestList [       testSeriesIdentities "explicit vs. ODE solution" identitiesExplODE,       testSeriesIdentities "transcendent functions of series" identitiesSeriesFunction,-      testSeriesIdentities "inverses of some series" identitiesInverses-{-+      testSeriesIdentities "inverses of some series" identitiesInverses,       HUnit.TestLabel "laws" $       HUnit.TestList $-         map testUnit $-            ("products of powers",     quickCheck (powerMult)) :-            ("power explicit vs. ODE", quickCheck (powerExplODE)) :-            []--}+         map testHolesIdentity identitiesHoles+         +++         (map testUnit $+            ("multiplicative holes",   quickCheck (holesMultiplicative 100)) :+            ("powers of series",       quickCheck (powerMultSeries 15)) :+            ("products of powers",     quickCheck (powerMult 30)) :+            ("power explicit vs. ODE", quickCheck (powerExplODE 50)) :+            ("inv vs. invDiff",        quickCheck (invDiff 15)) :+            [])     ]