numbers 3000.0.0.0 → 3000.1.0.0
raw patch · 11 files changed
+402/−356 lines, 11 filesdep +QuickCheckdep +test-frameworkdep +test-framework-quickcheck2PVP ok
version bump matches the API change (PVP)
Dependencies added: QuickCheck, test-framework, test-framework-quickcheck2
API changes (from Hackage documentation)
- Data.Number.BigFloat: instance Ord (BigFloat e)
+ Data.Number.BigFloat: instance Epsilon e => Ord (BigFloat e)
Files
- Data/Number/BigFloat.hs +15/−10
- Data/Number/CReal.hs +5/−6
- Data/Number/Dif.hs +7/−3
- Data/Number/Fixed.hs +22/−10
- Data/Number/FixedFunctions.hs +314/−312
- Data/Number/Interval.hs +1/−1
- Data/Number/Natural.hs +3/−3
- Data/Number/Symbolic.hs +5/−6
- Data/Number/Vectorspace.hs +3/−1
- LICENSE +1/−1
- numbers.cabal +26/−3
Data/Number/BigFloat.hs view
@@ -5,10 +5,10 @@ -- The numbers are stored in base 10. module Data.Number.BigFloat( BigFloat,- Epsilon, Eps1, EpsDiv10, Prec10, Prec50, PrecPlus20,+ Epsilon, Eps1, EpsDiv10, Prec10, Prec50, PrecPlus20 ) where-import Numeric(showSigned) +import Numeric(showSigned) import Data.Number.Fixed import qualified Data.Number.FixedFunctions as F @@ -19,7 +19,7 @@ -- but is more work. -- | Floating point number where the precision is determined by the type /e/. data BigFloat e = BF (Fixed e) Integer- deriving (Eq, Ord)+ deriving (Eq) instance (Epsilon e) => Show (BigFloat e) where showsPrec = showSigned showBF@@ -29,7 +29,7 @@ instance (Epsilon e) => Num (BigFloat e) where BF m1 e1 + BF m2 e2 = bf (m1' + m2') e where (m1', m2') = if e == e1 then (m1, m2 / base^(e-e2))- else (m1 / base^(e-e1), m2)+ else (m1 / base^(e-e1), m2) e = e1 `max` e2 -- Do - via negate BF m1 e1 * BF m2 e2 = bf (m1 * m2) (e1 + e2)@@ -41,12 +41,17 @@ instance (Epsilon e) => Real (BigFloat e) where toRational (BF e m) = toRational e * base^^m +instance (Epsilon e) => Ord (BigFloat e) where+ compare x y = compare (toRational x) (toRational y)+ instance (Epsilon e) => Fractional (BigFloat e) where recip (BF m e) = bf (base / m) (-(e + 1)) -- Take care not to lose precision for small numbers- fromRational x = if abs x < 1 then recip $ bf (fromRational (recip x)) 0- else bf (fromRational x) 0+ fromRational x+ | x == 0 || abs x >= 1 = bf (fromRational x) 0+ | otherwise = recip $ bf (fromRational (recip x)) 0 + -- normalizing constructor -- XXX The scaling is very inefficient bf :: (Epsilon e) => Fixed e -> Integer -> BigFloat e@@ -59,8 +64,8 @@ instance (Epsilon e) => RealFrac (BigFloat e) where properFraction x@(BF m e) = if e < 0 then (0, x)- else let (i, f) = properFraction (m * base^^e)- in (i, bf f 0)+ else let (i, f) = properFraction (m * base^^e)+ in (i, bf f 0) instance (Epsilon e) => Floating (BigFloat e) where pi = bf pi 0@@ -87,7 +92,7 @@ floatRange _ = (minBound, maxBound) decodeFloat x@(BF m e) = let d = floatDigits x- in (round $ m * base^d, fromInteger e - d)+ in (round $ m * base^d, fromInteger e - d) encodeFloat m e = bf (fromInteger m) (toInteger e) exponent (BF _ e) = fromInteger e significand (BF m _) = BF m 0@@ -99,7 +104,7 @@ isIEEE _ = False toFloat1 :: (Epsilon e) => (Rational -> Rational -> Rational) ->- BigFloat e -> BigFloat e+ BigFloat e -> BigFloat e toFloat1 f x@(BF m e) = fromRational $ f (precision m * scl) (toRational m * scl) where scl = base^^e
Data/Number/CReal.hs view
@@ -1,4 +1,3 @@-{-# OPTIONS -fglasgow-exts #-} -- ERA: Exact Real Arithmetic (version 1.0) -- -- A tolerably efficient and possibly correct implementation of the computable @@ -144,7 +143,7 @@ xr = x' p'; xn = 2^p'; g yn = round_uk ((yn*xr) % (2^p')) in round_uk (accumulate (iterate g xn) (take t ps) % (2^l2t))) where accumulate _ [] = 0 - accumulate [] _ = error "CReal.power_series.accumulate" + accumulate [] _ = error "CReal.power_series.accumulate" accumulate (x:xs) (c:cs) = let t = round_uk (c*(x % 1)) in if t == 0 then 0 else t + accumulate xs cs @@ -169,8 +168,8 @@ enumFromTo n e = takeWhile (<= e) $ iterate (+ 1)n enumFromThen n m = iterate (+(m-n)) n enumFromThenTo n m e = if m >= n then takeWhile (<= e) $ iterate (+(m-n)) n - else takeWhile (>= e) $ iterate (+(m-n)) n - + else takeWhile (>= e) $ iterate (+(m-n)) n + instance Real CReal where -- toRational x@(CR x') = x' n % 2^n where n = digitsToBits digits toRational _ = error "CReal.toRational" @@ -206,8 +205,8 @@ (s,ds') = let sgn = head ds == '-' in (sgn, if sgn then tail ds else ds) ds'' = take (max (d+1-length ds') 0) (repeat '0') ++ ds' (zs,fs) = splitAt (length ds'' -d) ds'' - fs' = case reverse $ dropWhile (== '0') $ reverse fs of - "" -> "0" + fs' = case reverse $ dropWhile (== '0') $ reverse fs of + "" -> "0" xs -> xs digitsToBits :: Int -> Int
Data/Number/Dif.hs view
@@ -61,9 +61,9 @@ -- |The 'deriv' function is a simple utility to take the -- derivative of a (single argument) function. -- It is simply defined as--- +-- -- > deriv f = val . df . f . dVar--- +-- deriv :: (Num a, Num b, Eq a, Eq b) => (Dif a -> Dif b) -> (a -> b) deriv f = val . df . f . dVar @@ -118,7 +118,7 @@ instance (Fractional a, Eq a) => Fractional (Dif a) where recip (C x) = C (recip x) recip (D x x') = ip- where ip = D (recip x) (-x' * ip * ip)+ where ip = D (recip x) (-x' * ip * ip) fromRational r = C (fromRational r) lift :: (Num a, Eq a) => [a -> a] -> Dif a -> Dif a@@ -177,3 +177,7 @@ isDenormalized = isDenormalized . val isNegativeZero = isNegativeZero . val isIEEE = isIEEE . val+ -- Set these to undefined rather than omit them to avoid compiler+ -- warnings.+ decodeFloat = undefined+ encodeFloat = undefined
Data/Number/Fixed.hs view
@@ -1,4 +1,9 @@-{-# OPTIONS_GHC -fglasgow-exts -fscoped-type-variables #-}+{-# LANGUAGE+ EmptyDataDecls,+ GeneralizedNewtypeDeriving,+ ScopedTypeVariables,+ Rank2Types #-}+ -- | Numbers with a fixed number of decimals. module Data.Number.Fixed( Fixed,@@ -24,7 +29,7 @@ instance (Epsilon e) => Epsilon (EpsDiv10 e) where eps e = eps (un e) / 10 where un :: EpsDiv10 e -> e- un = undefined+ un = undefined -- | Ten decimals. data Prec10@@ -46,7 +51,7 @@ instance (Epsilon e) => Epsilon (PrecPlus20 e) where eps e = 1e-20 * eps (un e) where un :: PrecPlus20 e -> e- un = undefined+ un = undefined ----------- @@ -82,7 +87,7 @@ convertFixed e@(F x) = f where f = F $ if feps > eeps then approx x feps else x feps = getEps f- eeps = getEps e+ eeps = getEps e getEps :: (Epsilon e) => Fixed e -> Rational getEps = eps . un@@ -93,11 +98,11 @@ showsPrec = showSigned showFixed where showFixed f@(F x) = showString $ show q ++ "." ++ decimals r e where q :: Integer- (q, r) = properFraction (x + e/2)- e = getEps f- decimals a e | e >= 1 = ""- | otherwise = intToDigit b : decimals c (10 * e)- where (b, c) = properFraction (10 * a)+ (q, r) = properFraction (x + e/2)+ e = getEps f+ decimals a e | e >= 1 = ""+ | otherwise = intToDigit b : decimals c (10 * e)+ where (b, c) = properFraction (10 * a) instance (Epsilon e) => Read (Fixed e) where readsPrec _ = readSigned readFixed@@ -136,6 +141,13 @@ isDenormalized _ = False isNegativeZero _ = False isIEEE _ = False+ -- Explicitly undefine these rather than omitting them; this+ -- prevents a compiler warning at least.+ floatRadix = undefined+ floatDigits = undefined+ floatRange = undefined+ decodeFloat = undefined+ encodeFloat = undefined ----------- @@ -143,4 +155,4 @@ dynamicEps :: forall a . Rational -> (forall e . Epsilon e => Fixed e -> a) -> Rational -> a dynamicEps r f v = loop (undefined :: Eps1) where loop :: forall x . (Epsilon x) => x -> a- loop e = if eps e <= r then f (fromRational v :: Fixed x) else loop (undefined :: EpsDiv10 x)+ loop e = if eps e <= r then f (fromRational v :: Fixed x) else loop (undefined :: EpsDiv10 x)
Data/Number/FixedFunctions.hs view
@@ -2,70 +2,70 @@ -- -- Module: ----- Fraction.hs+-- Fraction.hs -- -- Language: ----- Haskell+-- Haskell -- -- Description: Rational with transcendental functionalities -- ----- This is a generalized Rational in disguise. Rational, as a type--- synonim, could not be directly made an instance of any new class--- at all.--- But we would like it to be an instance of Transcendental, where--- trigonometry, hyperbolics, logarithms, etc. are defined.--- So here we are tiptoe-ing around, re-defining everything from--- scratch, before designing the transcendental functions -- which--- is the main motivation for this module.+-- This is a generalized Rational in disguise. Rational, as a type+-- synonim, could not be directly made an instance of any new class+-- at all.+-- But we would like it to be an instance of Transcendental, where+-- trigonometry, hyperbolics, logarithms, etc. are defined.+-- So here we are tiptoe-ing around, re-defining everything from+-- scratch, before designing the transcendental functions -- which+-- is the main motivation for this module. ----- Aside from its ability to compute transcendentals, Fraction--- allows for denominators zero. Unlike Rational, Fraction does--- not produce run-time errors for zero denominators, but use such--- entities as indicators of invalid results -- plus or minus--- infinities. Operations on fractions never fail in principle.+-- Aside from its ability to compute transcendentals, Fraction+-- allows for denominators zero. Unlike Rational, Fraction does+-- not produce run-time errors for zero denominators, but use such+-- entities as indicators of invalid results -- plus or minus+-- infinities. Operations on fractions never fail in principle. -- -- However, some function may compute slowly when both numerators--- and denominators of their arguments are chosen to be huge.--- For example, periodicity relations are utilized with large--- arguments in trigonometric functions to reduce the arguments--- to smaller values and thus improve on the convergence--- of continued fractions. Yet, if pi number is chosen to--- be extremely accurate then the reduced argument would--- become a fraction with huge numerator and denominator--- -- thus slowing down the entire computation of a trigonometric--- function.+-- and denominators of their arguments are chosen to be huge.+-- For example, periodicity relations are utilized with large+-- arguments in trigonometric functions to reduce the arguments+-- to smaller values and thus improve on the convergence+-- of continued fractions. Yet, if pi number is chosen to+-- be extremely accurate then the reduced argument would+-- become a fraction with huge numerator and denominator+-- -- thus slowing down the entire computation of a trigonometric+-- function. -- -- Usage: ----- When computation speed is not an issue and accuracy is important--- this module replaces some of the functionalities typically handled--- by the floating point numbers: trigonometry, hyperbolics, roots--- and some special functions. All computations, including definitions--- of the basic constants pi and e, can be carried with any desired--- accuracy. One suggested usage is for mathematical servers, where--- safety might be more important than speed. See also the module--- Numerus, which supports mixed arithmetic between Integer,--- Fraction and Cofra (Complex fraction), and returns complex--- legal answers in some cases where Fraction would produce--- infinities: log (-5), sqrt (-1), etc. +-- When computation speed is not an issue and accuracy is important+-- this module replaces some of the functionalities typically handled+-- by the floating point numbers: trigonometry, hyperbolics, roots+-- and some special functions. All computations, including definitions+-- of the basic constants pi and e, can be carried with any desired+-- accuracy. One suggested usage is for mathematical servers, where+-- safety might be more important than speed. See also the module+-- Numerus, which supports mixed arithmetic between Integer,+-- Fraction and Cofra (Complex fraction), and returns complex+-- legal answers in some cases where Fraction would produce+-- infinities: log (-5), sqrt (-1), etc. ----- +-- -- Required: ----- Haskell Prelude --- +-- Haskell Prelude+-- -- Author: ----- Jan Skibinski, Numeric Quest Inc.+-- Jan Skibinski, Numeric Quest Inc. -- -- Date: ----- 1998.08.16, last modified 2000.05.31--- +-- 1998.08.16, last modified 2000.05.31+-- -- See also bottom of the page for description of the format used--- for continued fractions, references, etc. +-- for continued fractions, references, etc. ------------------------------------------------------------------- module Data.Number.FixedFunctions where@@ -76,162 +76,164 @@ approx eps x = approxRational x eps --------------------------------------------------------------------- Category: Conversion--- from continued fraction to fraction and vice versa,--- from Taylor series to continued fraction.+-- Category: Conversion+-- from continued fraction to fraction and vice versa,+-- from Taylor series to continued fraction. --------------------------------------------------------------------type CF = [(Rational, Rational)]+type CF = [(Rational, Rational)] fromCF :: CF -> Rational fromCF x =- --- -- Convert finite continued fraction to fraction- -- evaluating from right to left. This is used- -- mainly for testing in conjunction with "toCF".- --- foldr g 1 x- where- g :: (Rational, Rational) -> Rational -> Rational- g u v = (fst u) + (snd u) / v+ --+ -- Convert finite continued fraction to fraction+ -- evaluating from right to left. This is used+ -- mainly for testing in conjunction with "toCF".+ --+ foldr g 1 x+ where+ g :: (Rational, Rational) -> Rational -> Rational+ g u v = (fst u) + (snd u) / v -toCF :: Rational -> CF+toCF :: Rational -> CF toCF x =- --- -- Convert fraction to finite continued fraction- --- toCF' x []- where- toCF' u lst =+ --+ -- Convert fraction to finite continued fraction+ --+ toCF' x []+ where+ toCF' u lst = case r of 0 -> reverse (((q%1),(0%1)):lst)- _ -> toCF' (b%r) (((q%1),(1%1)):lst) - where- a = numerator u- b = denominator u - (q,r) = quotRem a b + _ -> toCF' (b%r) (((q%1),(1%1)):lst)+ where+ a = numerator u+ b = denominator u+ (q,r) = quotRem a b approxCF :: Rational -> CF -> Rational approxCF eps [] = 0-approxCF eps x - --- -- Approximate infinite continued fraction x by fraction,- -- evaluating from left to right, and stopping when- -- accuracy eps is achieved, or when a partial numerator- -- is zero -- as it indicates the end of CF.- --- -- This recursive function relates continued fraction- -- to rational approximation.- --- = approxCF' eps x 0 1 1 q' p' 1- where- h = fst (x!!0)- (q', p') = x!!0- approxCF' eps x v2 v1 u2 u1 a' n - | abs (1 - f1/f) < eps = approx eps f- | a == 0 = approx eps f - | otherwise = approxCF' eps x v1 v u1 u a (n+1)- where- (b, a) = x!!n- u = b*u1 + a'*u2- v = b*v1 + a'*v2- f = u/v- f1 = u1/v1- +approxCF eps x+ --+ -- Approximate infinite continued fraction x by fraction,+ -- evaluating from left to right, and stopping when+ -- accuracy eps is achieved, or when a partial numerator+ -- is zero -- as it indicates the end of CF.+ --+ -- This recursive function relates continued fraction+ -- to rational approximation.+ --+ = approxCF' eps x 0 1 1 q' p' 1+ where+ h = fst (x!!0)+ (q', p') = x!!0+ approxCF' eps x v2 v1 u2 u1 a' n+ | abs (1 - f1/f) < eps = approx eps f+ | a == 0 = approx eps f+ | otherwise = approxCF' eps x v1 v u1 u a (n+1)+ where+ (b, a) = x!!n+ u = b*u1 + a'*u2+ v = b*v1 + a'*v2+ f = u/v+ f1 = u1/v1 ++-- Type signature determined by GHC.+fromTaylorToCF :: Fractional a => [a] -> a -> [(a, a)] fromTaylorToCF s x =- --- -- Convert infinite number of terms of Taylor expansion of - -- a function f(x) to an infinite continued fraction,- -- where s = [s0,s1,s2,s3....] is a list of Taylor- -- series coefficients, such that f(x)=s0 + s1*x + s2*x^2.... - --- -- Require: No Taylor coefficient is zero- --- zero:one:[higher m | m <- [2..]]- where- zero = (s!!0, s!!1 * x) - one = (1, -s!!2/s!!1 * x)- higher m = (1 + s!!m/s!!(m-1) * x, -s!!(m+1)/s!!m * x)- + --+ -- Convert infinite number of terms of Taylor expansion of+ -- a function f(x) to an infinite continued fraction,+ -- where s = [s0,s1,s2,s3....] is a list of Taylor+ -- series coefficients, such that f(x)=s0 + s1*x + s2*x^2....+ --+ -- Require: No Taylor coefficient is zero+ --+ zero:one:[higher m | m <- [2..]]+ where+ zero = (s!!0, s!!1 * x)+ one = (1, -s!!2/s!!1 * x)+ higher m = (1 + s!!m/s!!(m-1) * x, -s!!(m+1)/s!!m * x) + --------------------------------------------------------------------- Category: Auxiliaries+-- Category: Auxiliaries ------------------------------------------------------------------ -fac :: Integer -> Integer +fac :: Integer -> Integer fac = product . enumFromTo 1 integerRoot2 :: Integer -> Integer integerRoot2 1 = 1 integerRoot2 x = --- -- Biggest integer m, such that x - m^2 >= 0,- -- where x is a positive integer+ -- Biggest integer m, such that x - m^2 >= 0,+ -- where x is a positive integer -- integerRoot2' 0 x (x `div` 2) x where- integerRoot2' lo hi r y - | c > y = integerRoot2' lo r ((r + lo) `div` 2) y- | c == y = r- | otherwise = - if (r+1)^2 > y then- r- else- integerRoot2' r hi ((r + hi) `div` 2) y- where c = r^2+ integerRoot2' lo hi r y+ | c > y = integerRoot2' lo r ((r + lo) `div` 2) y+ | c == y = r+ | otherwise =+ if (r+1)^2 > y then+ r+ else+ integerRoot2' r hi ((r + hi) `div` 2) y+ where c = r^2 ------------------------------------------------------------------- -- Everything below is the instantiation of class Transcendental -- for type Rational. See also modules Cofra and Numerus. ----- Category: Constants +-- Category: Constants ------------------------------------------------------------------- pi :: Rational -> Rational pi eps =- --- -- pi with accuracy eps- --- -- Based on Ramanujan formula, as described in Ref. 3- -- Accuracy: extremely good, 10^-19 for one term of continued- -- fraction- --- (sqrt eps d) / (approxCF eps (fromTaylorToCF s x))- where- x = 1%(640320^3)::Rational- s = [((-1)^k*(fac (6*k))%((fac k)^3*(fac (3*k))))*((a*k+b)%c) | k<-[0..]]+ --+ -- pi with accuracy eps+ --+ -- Based on Ramanujan formula, as described in Ref. 3+ -- Accuracy: extremely good, 10^-19 for one term of continued+ -- fraction+ --+ (sqrt eps d) / (approxCF eps (fromTaylorToCF s x))+ where+ x = 1%(640320^3)::Rational+ s = [((-1)^k*(fac (6*k))%((fac k)^3*(fac (3*k))))*((a*k+b)%c) | k<-[0..]] a = 545140134- b = 13591409- c = 426880- d = 10005- + b = 13591409+ c = 426880+ d = 10005+ ------------------------------------------------------------------------ Category: Trigonometry+-- Category: Trigonometry --------------------------------------------------------------------- tan :: Rational -> Rational -> Rational tan eps 0 = 0 tan eps x- --- -- Tangent x computed with accuracy of eps.- -- - -- Trigonometric identities are used first to reduce- -- the value of x to a value from within the range of [-pi/2,pi/2]- --- | x >= half_pi' = tan eps (x - ((1+m)%1)*xpi)- | x <= -half_pi' = tan eps (x + ((1+m)%1)*xpi)- --- | absx > 1 = 2 * t/(1 - t^2)- | otherwise = approxCF eps (cf x) - where- absx = abs x - t = tan eps (x/2)- m = floor ((absx - half_pi)/ xpi)- xpi = pi eps- half_pi'= 158%100- half_pi = xpi * (1%2)- cf u = ((0%1,1%1):[((2*r + 1)/u, -1) | r <- [0..]])- + --+ -- Tangent x computed with accuracy of eps.+ --+ -- Trigonometric identities are used first to reduce+ -- the value of x to a value from within the range of [-pi/2,pi/2]+ --+ | x >= half_pi' = tan eps (x - ((1+m)%1)*xpi)+ | x <= -half_pi' = tan eps (x + ((1+m)%1)*xpi)+ --- | absx > 1 = 2 * t/(1 - t^2)+ | otherwise = approxCF eps (cf x)+ where+ absx = abs x+ t = tan eps (x/2)+ m = floor ((absx - half_pi)/ xpi)+ xpi = pi eps+ half_pi'= 158%100+ half_pi = xpi * (1%2)+ cf u = ((0%1,1%1):[((2*r + 1)/u, -1) | r <- [0..]])+ sin :: Rational -> Rational -> Rational sin eps 0 = 0 sin eps x = 2*t/(1 + t*t)@@ -242,179 +244,179 @@ cos eps 0 = 1 cos eps x = (1 - p)/(1 + p) where- t = tan eps (x/2) + t = tan eps (x/2) p = t*t- + atan :: Rational -> Rational -> Rational atan eps x- --- -- Inverse tangent of x with approximation eps- --- | x == 0 = 0- | x > 1 = (pi eps)/2 - atan eps (1/x)- | x < -1 = -(pi eps)/2 - atan eps (1/x)- | otherwise = approxCF eps ((0,x):[((2*m - 1),(m*x)^2) | m<- [1..]])- - + --+ -- Inverse tangent of x with approximation eps+ --+ | x == 0 = 0+ | x > 1 = (pi eps)/2 - atan eps (1/x)+ | x < -1 = -(pi eps)/2 - atan eps (1/x)+ | otherwise = approxCF eps ((0,x):[((2*m - 1),(m*x)^2) | m<- [1..]])++ asin :: Rational -> Rational -> Rational-asin eps x - --- -- Inverse sine of x with approximation eps- --- | x == 0 = 0- | abs x > 1 = error "Fraction.asin"- | x == 1 = (pi eps) * (1%2)- | x == -1 = (pi eps) * (-1%2)- | otherwise = atan eps (x / (sqrt eps (1 - x^2)))+asin eps x+ --+ -- Inverse sine of x with approximation eps+ --+ | x == 0 = 0+ | abs x > 1 = error "Fraction.asin"+ | x == 1 = (pi eps) * (1%2)+ | x == -1 = (pi eps) * (-1%2)+ | otherwise = atan eps (x / (sqrt eps (1 - x^2))) - + acos :: Rational -> Rational -> Rational-acos eps x - --- -- Inverse cosine of x with approximation eps- --- | x == 0 = (pi eps)*(1%2)- | abs x > 1 = error "Fraction.sin"- | x == 1 = 0- | x == -1 = pi eps- | otherwise = atan eps ((sqrt eps (1 - x^2)) / x)- +acos eps x+ --+ -- Inverse cosine of x with approximation eps+ --+ | x == 0 = (pi eps)*(1%2)+ | abs x > 1 = error "Fraction.sin"+ | x == 1 = 0+ | x == -1 = pi eps+ | otherwise = atan eps ((sqrt eps (1 - x^2)) / x)+ ------------------------------------------------------------------------ Category: Roots+-- Category: Roots ---------------------------------------------------------------------- + sqrt :: Rational -> Rational -> Rational sqrt eps x --- -- Square root of x with approximation eps- --- -- The CF pattern is: [(m,x-m^2),(2m,x-m^2),(2m,x-m^2)....]- -- where m is the biggest integer such that x-m^2 >= 0- --- | x < 0 = error "Fraction.sqrt"- | x == 0 = 0- | x < 1 = 1/(sqrt eps (1/x))- | otherwise = approxCF eps ((m,x-m^2):[(2*m,x-m^2) | r<-[0..]]) - where- m = (integerRoot2 (floor x))%1- + -- Square root of x with approximation eps+ --+ -- The CF pattern is: [(m,x-m^2),(2m,x-m^2),(2m,x-m^2)....]+ -- where m is the biggest integer such that x-m^2 >= 0+ --+ | x < 0 = error "Fraction.sqrt"+ | x == 0 = 0+ | x < 1 = 1/(sqrt eps (1/x))+ | otherwise = approxCF eps ((m,x-m^2):[(2*m,x-m^2) | r<-[0..]])+ where+ m = (integerRoot2 (floor x))%1+ ------------------------------------------------------------------------ Category: Exponentials and hyperbolics+-- Category: Exponentials and hyperbolics --------------------------------------------------------------------- exp :: Rational -> Rational -> Rational-exp eps x - --- -- Exponent of x with approximation eps- --- -- Based on Jacobi type continued fraction for exponential,- -- with fractional terms:- -- n == 0 ==> (1,x) - -- n == 1 ==> (1 -x/2, x^2/12) - -- n >= 2 ==> (1, x^2/(16*n^2 - 4))- -- For x outside [-1,1] apply identity exp(x) = (exp(x/2))^2- --- | x == 0 = 1- | x > 1 = (approxCF eps (f (x*(1%p))))^p- | x < (-1) = (approxCF eps (f (x*(1%q))))^q- | otherwise = approxCF eps (f x)- where- p = ceiling x- q = -(floor x)- f y = (1,y):(1-y/2,y^2/12):[(1,y^2/(16*n^2-4)) | n<-[2..]]- - +exp eps x+ --+ -- Exponent of x with approximation eps+ --+ -- Based on Jacobi type continued fraction for exponential,+ -- with fractional terms:+ -- n == 0 ==> (1,x)+ -- n == 1 ==> (1 -x/2, x^2/12)+ -- n >= 2 ==> (1, x^2/(16*n^2 - 4))+ -- For x outside [-1,1] apply identity exp(x) = (exp(x/2))^2+ --+ | x == 0 = 1+ | x > 1 = (approxCF eps (f (x*(1%p))))^p+ | x < (-1) = (approxCF eps (f (x*(1%q))))^q+ | otherwise = approxCF eps (f x)+ where+ p = ceiling x+ q = -(floor x)+ f y = (1,y):(1-y/2,y^2/12):[(1,y^2/(16*n^2-4)) | n<-[2..]]++ cosh :: Rational -> Rational -> Rational cosh eps x =- --- -- Hyperbolic cosine with approximation eps- --- (a + b)*(1%2)- where- a = exp eps x- b = 1/a+ --+ -- Hyperbolic cosine with approximation eps+ --+ (a + b)*(1%2)+ where+ a = exp eps x+ b = 1/a sinh :: Rational -> Rational -> Rational sinh eps x =- --- -- Hyperbolic sine with approximation eps- --- (a - b)*(1%2)- where- a = exp eps x- b = 1/a+ --+ -- Hyperbolic sine with approximation eps+ --+ (a - b)*(1%2)+ where+ a = exp eps x+ b = 1/a tanh :: Rational -> Rational -> Rational tanh eps x =- --- -- Hyperbolic tangent with approximation eps- --- (a - b)/ (a + b)- where- a = exp eps x- b = 1/a+ --+ -- Hyperbolic tangent with approximation eps+ --+ (a - b)/ (a + b)+ where+ a = exp eps x+ b = 1/a atanh :: Rational -> Rational -> Rational-atanh eps x - --- -- Inverse hyperbolic tangent with approximation eps- --- --- | x >= 1 = 1%0--- | x <= -1 = -1%0- | otherwise = (1%2) * (log eps ((1 + x) / (1 - x)))- +atanh eps x+ --+ -- Inverse hyperbolic tangent with approximation eps+ --++-- | x >= 1 = 1%0+-- | x <= -1 = -1%0+ | otherwise = (1%2) * (log eps ((1 + x) / (1 - x)))+ asinh :: Rational -> Rational -> Rational-asinh eps x - --- -- Inverse hyperbolic sine- ----- | x == 1%0 = 1%0--- | x == -1%0 = -1%0- | otherwise = log eps (x + (sqrt eps (x^2 + 1)))- +asinh eps x+ --+ -- Inverse hyperbolic sine+ --+-- | x == 1%0 = 1%0+-- | x == -1%0 = -1%0+ | otherwise = log eps (x + (sqrt eps (x^2 + 1)))+ acosh :: Rational -> Rational -> Rational acosh eps x- --- -- Inverse hyperbolic cosine- ----- | x == 1%0 = 1%0--- | x < 1 = 1%0- | otherwise = log eps (x + (sqrt eps (x^2 - 1)))- + --+ -- Inverse hyperbolic cosine+ --+-- | x == 1%0 = 1%0+-- | x < 1 = 1%0+ | otherwise = log eps (x + (sqrt eps (x^2 - 1)))+ ------------------------------------------------------------------------ Category: Logarithms+-- Category: Logarithms --------------------------------------------------------------------- log :: Rational -> Rational -> Rational log eps x- -- - -- Natural logarithm of strictly positive x - --- -- Based on Stieltjes type continued fraction for log (1+y)- -- (0,y):(1,y/2):[(1,my/(4m+2)),(1,(m+1)y/(4m+2)),....- -- (m >= 1, two elements per m)- -- Efficient only for x close to one. For larger x we recursively- -- apply the identity log(x) = log(x/2) + log(2)- --- | x <= 0 = error "Fraction.log"- | x < 1 = -log eps (1/x)- | x == 1 = 0- | otherwise =- case (scaled (x,0)) of- (1,s) -> (s%1) * approxCF eps (series 1)- (y,0) -> approxCF eps (series (y-1)) - (y,s) -> approxCF eps (series (y-1)) + (s%1)*approxCF eps (series 1)- where + --+ -- Natural logarithm of strictly positive x+ --+ -- Based on Stieltjes type continued fraction for log (1+y)+ -- (0,y):(1,y/2):[(1,my/(4m+2)),(1,(m+1)y/(4m+2)),....+ -- (m >= 1, two elements per m)+ -- Efficient only for x close to one. For larger x we recursively+ -- apply the identity log(x) = log(x/2) + log(2)+ --+ | x <= 0 = error "Fraction.log"+ | x < 1 = -log eps (1/x)+ | x == 1 = 0+ | otherwise =+ case (scaled (x,0)) of+ (1,s) -> (s%1) * approxCF eps (series 1)+ (y,0) -> approxCF eps (series (y-1))+ (y,s) -> approxCF eps (series (y-1)) + (s%1)*approxCF eps (series 1)+ where series :: Rational -> CF series u = (0,u):(1,u/2):[(1,u*((m+n)%(4*m + 2)))|m<-[1..],n<-[0,1]]- scaled :: (Rational,Integer) -> (Rational, Integer)+ scaled :: (Rational,Integer) -> (Rational, Integer) scaled (x, n)- | x == 2 = (1,n+1)- | x < 2 = (x, n)- | otherwise = scaled (x*(1%2), n+1)+ | x == 2 = (1,n+1)+ | x < 2 = (x, n)+ | otherwise = scaled (x*(1%2), n+1) - + --------------------------------------------------------------------------- -- References: --@@ -424,46 +426,46 @@ -- http:%www.mathsoft.com/asolve/constant/cntfrc/cntfrc.html -- 3. "Efficient on-line computation of real functions using exact floating -- point", by Peter John Potts, Imperial College--- http:%theory.doc.ic.ac.uk/~pjp/ieee.html+-- http:%theory.doc.ic.ac.uk/~pjp/ieee.html -------------------------------------------------------------------------- -------------------------------------------------------------------------- --- The following representation of continued fractions is used:+-- The following representation of continued fractions is used: ----- Continued fraction: CF representation:--- ================== ====================--- b0 + a0+-- Continued fraction: CF representation:+-- ================== ====================+-- b0 + a0 -- ------- ==> [(b0, a0), (b1, a1), (b2, a2).....] -- b1 + a1 -- ------- -- b2 + ... ----- where "a's" and "b's" are Rationals.--- --- Many continued fractions could be represented by much simpler form--- [b1,b2,b3,b4..], where all coefficients "a" would have the same value 1--- and would not need to be explicitely listed; and the coefficients "b"--- could be chosen as integers.--- However, there are some useful continued fractions that are--- given with fraction coefficients: "a", "b" or both.--- A fractional form can always be converted to an integer form, but--- a conversion process is not always simple and such an effort is not--- always worth of the achieved savings in the storage space or the--- computational efficiency. +-- where "a's" and "b's" are Rationals. --+-- Many continued fractions could be represented by much simpler form+-- [b1,b2,b3,b4..], where all coefficients "a" would have the same value 1+-- and would not need to be explicitely listed; and the coefficients "b"+-- could be chosen as integers.+-- However, there are some useful continued fractions that are+-- given with fraction coefficients: "a", "b" or both.+-- A fractional form can always be converted to an integer form, but+-- a conversion process is not always simple and such an effort is not+-- always worth of the achieved savings in the storage space or the+-- computational efficiency.+-- ---------------------------------------------------------------------------- -- -- Copyright: ----- (C) 1998 Numeric Quest, All rights reserved+-- (C) 1998 Numeric Quest, All rights reserved -- -- <jans@numeric-quest.com> ----- http://www.numeric-quest.com +-- http://www.numeric-quest.com -- -- License: ----- GNU General Public License, GPL--- +-- GNU General Public License, GPL+-- -----------------------------------------------------------------------------
Data/Number/Interval.hs view
@@ -40,6 +40,6 @@ instance (Ord a, Fractional a) => Fractional (Interval a) where I l h / I l' h' | signum l' == signum h' && l' /= 0 = I (minimum xs) (maximum xs)- | otherwise = error "Interval: division by 0"+ | otherwise = error "Interval: division by 0" where xs = [l/l', l/h', h/l', h/h'] fromRational r = I l l where l = fromRational r
Data/Number/Natural.hs view
@@ -62,10 +62,10 @@ -- Not the most efficient version, but efficiency isn't the point of this module. :) quotRem x y = if x < y then- (0, x)+ (0, x) else- let (q, r) = quotRem (x-y) y- in (1+q, r)+ let (q, r) = quotRem (x-y) y+ in (1+q, r) div = quot mod = rem toInteger Z = 0
Data/Number/Symbolic.hs view
@@ -7,7 +7,6 @@ import Data.Char(isAlpha) import Data.Maybe(fromMaybe)-import Debug.Trace -- | Symbolic numbers over some base type for the literals. data Sym a = Con a | App String ([a]->a) [Sym a]@@ -36,7 +35,7 @@ subst :: (Num a, Eq a) => String -> Sym a -> Sym a -> Sym a subst _ _ e@(Con _) = e subst x v e@(App x' _ []) | x == x' = v- | otherwise = e+ | otherwise = e subst x v (App s f es) = case map (subst x v) es of [e] -> unOp (\ x -> f [x]) s e@@ -56,10 +55,10 @@ showParen (p>q) (showsPrec ql x . showString op . showsPrec qr y) where (ql, q, qr) = fromMaybe (9,9,9) $ lookup op [ ("**", (9,8,8)),- ("/", (7,7,8)),- ("*", (7,7,8)),- ("+", (6,6,7)),- ("-", (6,6,7))]+ ("/", (7,7,8)),+ ("*", (7,7,8)),+ ("+", (6,6,7)),+ ("-", (6,6,7))] showsPrec p (App "negate" _ [x]) = showParen (p>=6) (showString "-" . showsPrec 7 x) showsPrec p (App f _ xs) =
Data/Number/Vectorspace.hs view
@@ -1,4 +1,6 @@-{-# OPTIONS_GHC -fglasgow-exts #-}+{-# LANGUAGE+ FunctionalDependencies,+ MultiParamTypeClasses #-} module Data.Number.Vectorspace(Vectorspace(..)) where -- |Class of vector spaces /v/ with scalar /s/.
LICENSE view
@@ -1,6 +1,6 @@ Copyright (c) 2007-2012 Lennart Augustsson, Russell O'Connor, Richard Smith,-Daniel Wagner, Dan Burton+Daniel Wagner, Dan Burton, Michael Orlitzky All rights reserved.
numbers.cabal view
@@ -1,5 +1,5 @@ Name: numbers-Version: 3000.0.0.0+Version: 3000.1.0.0 License: BSD3 License-file: LICENSE Author: Lennart Augustsson@@ -27,11 +27,34 @@ tag: numbers-3000.0.0.0 Library- Build-Depends: base >= 3 && < 5+ Build-Depends:+ base >= 3 && < 5+ Exposed-modules: Data.Number.Symbolic Data.Number.Dif Data.Number.CReal Data.Number.Fixed Data.Number.Interval Data.Number.BigFloat Data.Number.Natural- Other-modules: Data.Number.Vectorspace Data.Number.FixedFunctions+ Other-modules:+ Data.Number.Vectorspace+ Data.Number.FixedFunctions + Ghc-Options:+ -Wall+ -fno-warn-name-shadowing+ -fno-warn-unused-binds+ -fno-warn-unused-matches+ -fno-warn-incomplete-patterns+ -fno-warn-overlapping-patterns+ -fno-warn-type-defaults++test-suite testsuite+ type: exitcode-stdio-1.0+ hs-source-dirs: . test+ main-is: TestSuite.hs+ build-depends:+ base >= 3 && < 5,+ -- Additional test dependencies.+ QuickCheck == 2.*,+ test-framework == 0.6.*,+ test-framework-quickcheck2 == 0.2.*