-- Module: -- -- Fraction.hs -- -- Language: -- -- 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. -- -- 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. -- -- 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. -- -- -- Required: -- -- Haskell Prelude -- -- Author: -- -- Jan Skibinski, Numeric Quest Inc. -- -- Date: -- -- 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. ------------------------------------------------------------------- module Fraction where import Data.Ratio infix 7 :-: ------------------------------------------------------------------- -- Category: Basics ------------------------------------------------------------------- data Fraction = Integer :-: Integer deriving (Eq) num, den :: Fraction -> Integer num (x:-:_) = x den (_:-:y) = y reduce :: Fraction -> Fraction reduce (x:-:0) | x < 0 = (-1):-:0 | otherwise = 1:-:0 reduce (x:-:y) = (u `quot` d) :-: (v `quot` d) where d = gcd u v (u,v) | y < 0 = (-x,-y) | otherwise = (x,y) (//) :: Integer -> Integer -> Fraction x // y = reduce (x:-:y) approx :: Fraction -> Fraction -> Fraction approx _ (x:-:0) = x//0 approx eps x = simplest (x-eps) (x+eps) where simplest y z | z < y = simplest z y | y == z = y | y > 0 = simplest' (num y) (den y) (num z) (den z) | z < 0 = - simplest' (-(num z)) (den z) (-(num y)) (den y) | otherwise = 0 :-: 1 simplest' n d n' d' -- assumes 0 < n//d < n'//d' | r == 0 = q :-: 1 | q /= q' = (q+1) :-: 1 | otherwise = (q*n''+d'') :-: n'' where (q,r) = quotRem n d (q',r') = quotRem n' d' (n'':-:d'') = simplest' d' r' d r ------------------------------------------------------------------- -- Category: Instantiation of some Prelude classes ------------------------------------------------------------------- instance Read Fraction where readsPrec p = readParen (p > 7) (\r -> [(x//y,u) | (x,s) <- reads r, ("//",t) <- lex s, (y,u) <- reads t ]) instance Show Fraction where showsPrec p (x:-:y) | y == 1 = showsPrec p x | otherwise = showParen (p > 7) (shows x . showString "/" . shows y) instance Ord Fraction where compare (x:-:y) (x':-:y') = compare (x*y') (x'*y) instance Num Fraction where (x:-:y) + (x':-:y') = reduce ((x*y' + x'*y):-:(y*y')) (x:-:y) - (x':-:y') = reduce ((x*y' - x'*y):-:(y*y')) (x:-:y) * (x':-:y') = reduce ((x*x') :-: (y*y')) negate (x:-:y) = negate x :-: y abs (x:-:y) = abs x :-: y signum (x:-:_) = signum x :-: 1 fromInteger n = fromInteger n :-: 1 instance Fractional Fraction where (x:-:0) / (x':-:0) = ((signum x * signum x'):-:0) (_:-:_) / (_:-:0) = (0:-:1) (x:-:0) / (_:-:_) = (x:-:0) (x:-:y) / (x':-:y') = reduce ((x*y') :-: (y*x')) recip (x:-:y) = if x < 0 then (-y) :-: (-x) else y :-: x fromRational a = x :-: y where x = numerator a y = denominator a instance Real Fraction where toRational (_ :-: 0) = toRational ((0::Int)%(1::Int)) -- or shoud we return some huge number instead? toRational (x :-: y) = toRational (x % y) instance RealFrac Fraction where properFraction (x :-: y) = (fromInteger q, r :-: y) where (q,r) = quotRem x y instance Enum Fraction where toEnum = fromIntegral fromEnum = truncate -- dubious enumFrom = numericEnumFrom enumFromTo = numericEnumFromTo enumFromThen = numericEnumFromThen enumFromThenTo = numericEnumFromThenTo numericEnumFrom :: Real a => a -> [a] numericEnumFromThen :: Real a => a -> a -> [a] numericEnumFromTo :: Real a => a -> a -> [a] numericEnumFromThenTo :: Real a => a -> a -> a -> [a] -- -- Prelude does not export these, so here are the copies numericEnumFrom n = n : (numericEnumFrom $! (n+1)) numericEnumFromThen n m = iterate ((m-n)+) n numericEnumFromTo n m = takeWhile (<= m) (numericEnumFrom n) numericEnumFromThenTo n n' m = takeWhile p (numericEnumFromThen n n') where p | n' >= n = (<= m) | otherwise = (>= m) ------------------------------------------------------------------ -- Category: Conversion -- from continued fraction to fraction and vice versa, -- from Taylor series to continued fraction. ------------------------------------------------------------------- type CF = [(Fraction, Fraction)] fromCF :: CF -> Fraction 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//1) x where g :: (Fraction, Fraction) -> Fraction -> Fraction g u v = (fst u) + (snd u)/v toCF :: Fraction -> CF toCF (u:-:0) = [(u//0,0//1)] toCF x = -- -- 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 = num u b = den u (q,r) = quotRem a b approxCF :: Fraction -> CF -> Fraction approxCF _ [] = 0//1 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. -- | den h == 0 = h | otherwise = approxCF' eps x 0 1 1 q' p' 1 where h = fst (x!!0) (q', p') = x!!0 approxCF' ept y v2 v1 u2 u1 a' n | abs (1 - f1/f) < ept = approx ept f | a == 0 = approx ept f | otherwise = approxCF' ept y v1 v u1 u a (n+1) where (b, a) = y!!n u = b*u1 + a'*u2 v = b*v1 + a'*v2 f = u/v f1 = u1/v1 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 -- -- It is an application of Euler's continued fraction formula. -- 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) fromFraction :: Fraction -> Double fromFraction = fromRational . toRational ------------------------------------------------------------------ -- Category: Auxiliaries ------------------------------------------------------------------ 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 -- 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::Int) > y then r else integerRoot2' r hi ((r + hi) `div` 2) y where c = r^(2::Int) ------------------------------------------------------------------ -- Category: Class Transcendental -- -- This class declares functions for three data types: -- Fraction, Cofraction (complex fraction) and Numerus -- - a generalization of Integer, Fraction and Cofraction. ------------------------------------------------------------------ class Transcendental a where pi' :: Fraction -> a tan' :: Fraction -> a -> a sin' :: Fraction -> a -> a cos' :: Fraction -> a -> a atan' :: Fraction -> a -> a asin' :: Fraction -> a -> a acos' :: Fraction -> a -> a sqrt' :: Fraction -> a -> a root' :: Fraction -> a-> Integer -> a power' :: Fraction -> a -> a -> a exp' :: Fraction -> a -> a tanh' :: Fraction -> a -> a sinh' :: Fraction -> a -> a cosh' :: Fraction -> a -> a atanh' :: Fraction -> a -> a asinh' :: Fraction -> a -> a acosh' :: Fraction -> a -> a log' :: Fraction -> a -> a decimal :: Integer -> a -> IO () ------------------------------------------------------------------- -- Everything below is the instantiation of class Transcendental -- for type Fraction. See also modules Cofra and Numerus. -- -- Category: Constants ------------------------------------------------------------------- instance Transcendental Fraction where 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::Int))::Fraction s = [((-1)^k*(fac (6*k))//((fac k)^(3::Int)*(fac (3*k))))*((a*k+b)//c) | k<-[0..]] a = 545140134 b = 13591409 c = 426880 d = 10005 --------------------------------------------------------------------- -- Category: Trigonometry --------------------------------------------------------------------- tan' _ 0 = 0 tan' _ (_:-:0) = 1//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)*p) | x <= -half_pi' = tan' eps (x + ((1+m)//1)*p) --- | absx > 1 = 2 * t/(1 - t^2) | otherwise = approxCF eps (cf x) where absx = abs x _ = tan' eps (x/2) m = floor ((absx - half_pi)/ p) p = pi' eps half_pi'= 158//100 half_pi = p * (1//2) cf u = ((0//1,1//1):[((2*r + 1)/u, -1) | r <- [0..]]) sin' _ 0 = 0 sin' _ (_:-:0)= 1//0 sin' eps x = 2*t/(1 + t*t) where t = tan' eps (x/2) cos' _ 0 = 1 cos' _ (_:-:0)= 1//0 cos' eps x = (1 - p)/(1 + p) where t = tan' eps (x/2) p = t*t atan' eps x -- -- Inverse tangent of x with approximation eps -- | x == 1//0 = (pi' eps)/2 | x == (-1//0) = -(pi' eps)/2 | 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::Int)) | m<- [1..]]) asin' eps x -- -- Inverse sine of x with approximation eps -- | x == 0 = 0//1 | abs x > 1 = 1//0 | x == 1 = (pi' eps) *(1//2) | x == -1 = (pi' eps) * ((-1)//2) | otherwise = atan' eps (x / (sqrt' eps (1 - x^(2::Int)))) acos' eps x -- -- Inverse cosine of x with approximation eps -- | x == 0 = (pi' eps)*(1//2) | abs x > 1 = 1//0 | x == 1 = 0//1 | x == -1 = pi' eps | otherwise = atan' eps ((sqrt' eps (1 - x^(2::Int))) / x) --------------------------------------------------------------------- -- Category: Roots --------------------------------------------------------------------- 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 == 1//0 = 1//0 | x < 0 = 1//0 | x == 0 = 0 | x < 1 = 1/(sqrt' eps (1/x)) | otherwise = approxCF eps ((m,x-m^(2::Int)):[(2*m,x-m^(2::Int)) | _<-[(0::Integer)..]]) where m = (integerRoot2 (floor x))//1 root' eps x k -- -- k-th root of positive number x with approximation eps -- | x == (1//0) = 1//0 | x < 0 = 1//0 | x == 0 = 0 | k == 0 = 1//0 | otherwise = exp' eps ((log' eps x) * (1//k)) --------------------------------------------------------------------- -- Category: Powers --------------------------------------------------------------------- power' eps x y -- -- x to power of y with approximation eps -- | x == (1//0) = 1//0 | x < 0 = 1//0 | x == 0 = 0 | y == 0 = 1 | y == (1//0) = 1//0 | y == (-1//0) = 0 | otherwise = exp' eps (y * (log' eps x)) --------------------------------------------------------------------- -- Category: Exponentials and hyperbolics --------------------------------------------------------------------- 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 == 1//0 = 1//0 | x == (-1//0) = 0 | 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::Int)/12):[(1,y^(2::Int)/(16*n^(2::Int)-4)) | n<-[2..]] cosh' eps x = -- -- Hyperbolic cosine with approximation eps -- (a + b)*(1//2) where a = exp' eps x b = 1/a sinh' eps x = -- -- Hyperbolic sine with approximation eps -- (a - b)*(1//2) where a = exp' eps x b = 1/a tanh' eps x = -- -- Hyperbolic tangent with approximation eps -- (a - b)/ (a + b) where a = exp' eps x b = 1/a 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' eps x -- -- Inverse hyperbolic sine -- | x == 1//0 = 1//0 | x == -1//0 = -1//0 | otherwise = log' eps (x + (sqrt' eps (x^(2::Int) + 1))) acosh' eps x -- -- Inverse hyperbolic cosine -- | x == 1//0 = 1//0 | x < 1 = 1//0 | otherwise = log' eps (x + (sqrt' eps (x^(2::Int) - 1))) --------------------------------------------------------------------- -- Category: Logarithms --------------------------------------------------------------------- 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 == 1//0 = 1//0 | x <= 0 = -1//0 | 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 :: Fraction -> CF series u = (0,u):(1,u/2):[(1,u*((m+n)//(4*m + 2)))|m<-[1..],n<-[0,1]] scaled :: (Fraction,Integer) -> (Fraction, Integer) scaled (y, n) | y == 2 = (1,n+1) | y < 2 = (y, n) | otherwise = scaled (y*(1//2), n+1) --------------------------------------------------------------------- -- Category: IO --------------------------------------------------------------------- decimal _ (u:-:0) = putStr (show u++"//0") decimal n x -- -- Print Fraction with an accuracy to n decimal places, -- or symbols +/- 1//0 for infinities. | n <= 0 = decimal 1 x | x < 0 = putStr (g (-v*10) (den x) n ("-"++show (-u) ++".")) | otherwise = putStr (g (v*10) (den x) n (show u++".")) where (u, v) = quotRem (num x) (den x) g _ _ 0 str = str g y z m str = case (p, q) of (_,0) -> str ++ show p (_,_) -> g (q*10) z (m-1) (str ++ show p) where (p, q) = quotRem y z --------------------------------------------------------------------------- -- References: -- -- 1. Classical Gosper notes on continued fraction arithmetic: -- http://www.inwap.com/pdp10/hbaker/hakmem/cf.html -- 2. Pages on numerical constants represented as continued fractions: -- 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 -------------------------------------------------------------------------- -------------------------------------------------------------------------- -- The following representation of continued fractions is used: -- -- Continued fraction: CF representation: -- ================== ==================== -- b0 + a0 -- ------- ==> [(b0, a0), (b1, a1), (b2, a2).....] -- b1 + a1 -- ------- -- b2 + ... -- -- where "a's" and "b's" are Fractions. -- -- 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 -- -- <jans@numeric-quest.com> -- -- http://www.numeric-quest.com -- -- License: -- -- GNU General Public License, GPL -- -----------------------------------------------------------------------------