diff --git a/numeric-prelude.cabal b/numeric-prelude.cabal --- a/numeric-prelude.cabal +++ b/numeric-prelude.cabal @@ -1,5 +1,5 @@ Name: numeric-prelude -Version: 0.4.0.3 +Version: 0.4.1 License: BSD3 License-File: LICENSE Author: Dylan Thurston , Henning Thielemann , 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 diff --git a/src/Algebra/Additive.hs b/src/Algebra/Additive.hs --- a/src/Algebra/Additive.hs +++ b/src/Algebra/Additive.hs @@ -58,6 +58,7 @@ -} class C a where + {-# MINIMAL zero, (+), ((-) | negate) #-} -- | zero element of the vector space zero :: a -- | add and subtract elements diff --git a/src/Algebra/Algebraic.hs b/src/Algebra/Algebraic.hs --- a/src/Algebra/Algebraic.hs +++ b/src/Algebra/Algebraic.hs @@ -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) diff --git a/src/Algebra/Field.hs b/src/Algebra/Field.hs --- a/src/Algebra/Field.hs +++ b/src/Algebra/Field.hs @@ -68,6 +68,7 @@ -} class (Ring.C a) => C a where + {-# MINIMAL recip | (/) #-} (/) :: a -> a -> a recip :: a -> a fromRational' :: Rational -> a diff --git a/src/Algebra/IntegralDomain.hs b/src/Algebra/IntegralDomain.hs --- a/src/Algebra/IntegralDomain.hs +++ b/src/Algebra/IntegralDomain.hs @@ -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) diff --git a/src/Algebra/RealRing.hs b/src/Algebra/RealRing.hs --- a/src/Algebra/RealRing.hs +++ b/src/Algebra/RealRing.hs @@ -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 diff --git a/src/Algebra/Ring.hs b/src/Algebra/Ring.hs --- a/src/Algebra/Ring.hs +++ b/src/Algebra/Ring.hs @@ -65,6 +65,7 @@ -} class (Additive.C a) => C a where + {-# MINIMAL (*), (one | fromInteger) #-} (*) :: a -> a -> a one :: a fromInteger :: Integer -> a diff --git a/src/Algebra/Transcendental.hs b/src/Algebra/Transcendental.hs --- a/src/Algebra/Transcendental.hs +++ b/src/Algebra/Transcendental.hs @@ -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 diff --git a/src/Algebra/Units.hs b/src/Algebra/Units.hs --- a/src/Algebra/Units.hs +++ b/src/Algebra/Units.hs @@ -70,6 +70,7 @@ -} class (Integral.C a) => C a where + {-# MINIMAL isUnit, (stdUnit | stdUnitInv) #-} isUnit :: a -> Bool stdAssociate, stdUnit, stdUnitInv :: a -> a diff --git a/src/MathObj/Polynomial/Core.hs b/src/MathObj/Polynomial/Core.hs --- a/src/MathObj/Polynomial/Core.hs +++ b/src/MathObj/Polynomial/Core.hs @@ -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 diff --git a/src/MathObj/PowerSeries/Core.hs b/src/MathObj/PowerSeries/Core.hs --- a/src/MathObj/PowerSeries/Core.hs +++ b/src/MathObj/PowerSeries/Core.hs @@ -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 diff --git a/src/MathObj/PowerSeries2.hs b/src/MathObj/PowerSeries2.hs --- a/src/MathObj/PowerSeries2.hs +++ b/src/MathObj/PowerSeries2.hs @@ -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 diff --git a/src/MathObj/PowerSeries2/Core.hs b/src/MathObj/PowerSeries2/Core.hs --- a/src/MathObj/PowerSeries2/Core.hs +++ b/src/MathObj/PowerSeries2/Core.hs @@ -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 diff --git a/test/Test/MathObj/Polynomial.hs b/test/Test/MathObj/Polynomial.hs --- a/test/Test/MathObj/Polynomial.hs +++ b/test/Test/MathObj/Polynomial.hs @@ -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) : [] diff --git a/test/Test/MathObj/PowerSeries.hs b/test/Test/MathObj/PowerSeries.hs --- a/test/Test/MathObj/PowerSeries.hs +++ b/test/Test/MathObj/PowerSeries.hs @@ -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)) : + []) ]