diff --git a/Data/Number/BigFloat.hs b/Data/Number/BigFloat.hs --- a/Data/Number/BigFloat.hs +++ b/Data/Number/BigFloat.hs @@ -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 diff --git a/Data/Number/CReal.hs b/Data/Number/CReal.hs --- a/Data/Number/CReal.hs +++ b/Data/Number/CReal.hs @@ -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 diff --git a/Data/Number/Dif.hs b/Data/Number/Dif.hs --- a/Data/Number/Dif.hs +++ b/Data/Number/Dif.hs @@ -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 diff --git a/Data/Number/Fixed.hs b/Data/Number/Fixed.hs --- a/Data/Number/Fixed.hs +++ b/Data/Number/Fixed.hs @@ -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) diff --git a/Data/Number/FixedFunctions.hs b/Data/Number/FixedFunctions.hs --- a/Data/Number/FixedFunctions.hs +++ b/Data/Number/FixedFunctions.hs @@ -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 -- -- -- --- http://www.numeric-quest.com +-- http://www.numeric-quest.com -- -- License: -- --- GNU General Public License, GPL --- +-- GNU General Public License, GPL +-- ----------------------------------------------------------------------------- diff --git a/Data/Number/Interval.hs b/Data/Number/Interval.hs --- a/Data/Number/Interval.hs +++ b/Data/Number/Interval.hs @@ -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 diff --git a/Data/Number/Natural.hs b/Data/Number/Natural.hs --- a/Data/Number/Natural.hs +++ b/Data/Number/Natural.hs @@ -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 diff --git a/Data/Number/Symbolic.hs b/Data/Number/Symbolic.hs --- a/Data/Number/Symbolic.hs +++ b/Data/Number/Symbolic.hs @@ -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) = diff --git a/Data/Number/Vectorspace.hs b/Data/Number/Vectorspace.hs --- a/Data/Number/Vectorspace.hs +++ b/Data/Number/Vectorspace.hs @@ -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/. diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -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. diff --git a/numbers.cabal b/numbers.cabal --- a/numbers.cabal +++ b/numbers.cabal @@ -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.*