{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2006 Maintainer : numericprelude@henning-thielemann.de Stability : provisional Portability : requires multi-parameter type classes Fixed point numbers. They are implemented as ratios with fixed denominator. Many routines fail for some arguments. When they work, they can be useful for obtaining approximations of some constants. We have not paid attention to rounding errors and thus some of the trailing digits may be wrong. -} module Number.FixedPoint where import qualified Algebra.RealField as RealField import qualified Algebra.Additive as Additive -- import qualified Algebra.ZeroTestable as ZeroTestable import qualified Algebra.Transcendental as Trans import qualified MathObj.PowerSeries.Example as PSE import NumericPrelude.List (mapLast, ) import Data.Function.HT (powerAssociative, ) import Data.List.HT (dropWhileRev, padLeft, ) import Data.Maybe.HT (toMaybe, ) import Data.List (transpose, unfoldr, ) import Data.Char (intToDigit, ) import PreludeBase import NumericPrelude hiding (recip, sqrt, exp, sin, cos, tan, fromRational') import qualified NumericPrelude as NP {- ** Conversion -} {- ** other number types -} fromFloat :: RealField.C a => Integer -> a -> Integer fromFloat den x = round (x * NP.fromInteger den) -- | denominator conversion fromFixedPoint :: Integer -> Integer -> Integer -> Integer fromFixedPoint denDst denSrc x = div (x*denDst) denSrc {- ** text -} {- | very efficient because it can make use of the decimal output of 'show' -} showPositionalDec :: Integer -> Integer -> String showPositionalDec den = liftShowPosToInt $\x -> let packetSize = 50 -- process digits in packets of this size basis = ringPower packetSize 10 (int,frac) = toPositional basis den x in show int ++ "." ++ concat (mapLast (dropWhileRev ('0'==)) (map (padLeft '0' packetSize . show) frac)) showPositionalHex :: Integer -> Integer -> String showPositionalHex = showPositionalBasis 16 showPositionalBin :: Integer -> Integer -> String showPositionalBin = showPositionalBasis 2 showPositionalBasis :: Integer -> Integer -> Integer -> String showPositionalBasis basis den = liftShowPosToInt$ \x -> let (int,frac) = toPositional basis den x in show int ++ "." ++ map (intToDigit . fromInteger) frac liftShowPosToInt :: (Integer -> String) -> (Integer -> String) liftShowPosToInt f n = if n>=0 then f n else '-' : f (-n) toPositional :: Integer -> Integer -> Integer -> (Integer, [Integer]) toPositional basis den x = let (int, frac) = divMod x den in (int, unfoldr (\rm -> toMaybe (rm/=0) (divMod (basis*rm) den)) frac) {- * Additive -} add :: Integer -> Integer -> Integer -> Integer add _ = (+) sub :: Integer -> Integer -> Integer -> Integer sub _ = (-) {- * Ring -} mul :: Integer -> Integer -> Integer -> Integer mul den x y = div (x*y) den {- * Field -} divide :: Integer -> Integer -> Integer -> Integer divide den x y = div (x*den) y recip :: Integer -> Integer -> Integer recip den x = div (den^2) x {- * Algebra -} {- Newton's method for computing roots. -} magnitudes :: [Integer] magnitudes = concat (transpose [iterate (^2) 4, iterate (^2) 8]) {- Maybe we can speed up the algorithm by calling sqrt recursively on deflated arguments. -} sqrt :: Integer -> Integer -> Integer sqrt den x = let xden = x*den initial = fst (head (dropWhile ((<= xden) . snd) (zip magnitudes (tail (tail magnitudes))))) approxs = iterate (\y -> div (y + div xden y) 2) initial isRoot y = y^2 <= xden && xden < (y+1)^2 in head (dropWhile (not . isRoot) approxs) -- bug: needs too long: root (12::Int) (fromIntegerBase 10 1000 2) root :: Integer -> Integer -> Integer -> Integer root n den x = let n1 = n-1 xden = x * den^n1 initial = fst (head (dropWhile ((\y -> y^n <= xden) . snd) (zip magnitudes (tail magnitudes)))) approxs = iterate (\y -> div (n1*y + div xden (y^n1)) n) initial isRoot y = y^n <= xden && xden < (y+1)^n in head (dropWhile (not . isRoot) approxs) {- * Transcendental -} -- very simple evaluation by power series with lots of rounding errors evalPowerSeries :: [Rational] -> Integer -> Integer -> Integer evalPowerSeries series den x = let powers = iterate (mul den x) den summands = zipWith (\c p -> round (c * fromInteger p)) series powers in sum (map snd (takeWhile (\(c,s) -> s/=0 || c==0) (zip series summands))) cos, sin, tan :: Integer -> Integer -> Integer cos = evalPowerSeries PSE.cos sin = evalPowerSeries PSE.sin -- tan will suffer from inaccuracies for small cosine tan den x = divide den (sin den x) (cos den x) -- it must abs x <= den arctanSmall :: Integer -> Integer -> Integer arctanSmall = evalPowerSeries PSE.atan -- will fail for large inputs arctan :: Integer -> Integer -> Integer arctan den x = let estimate = fromFloat den (Trans.atan (NP.fromRational' (x % den)) :: Double) tanEst = tan den estimate residue = divide den (x-tanEst) (den + mul den x tanEst) in estimate + arctanSmall den residue piConst :: Integer -> Integer piConst den = let den4 = 4*den stArcTan k x = let d = k*den4 in arctanSmall d (div d x) in {- formula 4 * (8 * arctan (1/10) - arctan (1/239) - 4 * arctan (1/515)) from "Bartsch: Mathematische Formeln" -} -- (stArcTan 8 10 - stArcTan 1 239 - stArcTan 4 515) -- formula by Stoermer (stArcTan 44 57 + stArcTan 7 239 - stArcTan 12 682 + stArcTan 24 12943) expSmall :: Integer -> Integer -> Integer expSmall = evalPowerSeries PSE.exp eConst :: Integer -> Integer eConst den = expSmall den den recipEConst :: Integer -> Integer recipEConst den = expSmall den (-den) exp :: Integer -> Integer -> Integer exp den x = let den2 = div den 2 (int,frac) = divMod (x + den2) den expFrac = expSmall den (frac-den2) in case compare int 0 of EQ -> expFrac GT -> powerAssociative (mul den) expFrac (eConst den) int LT -> powerAssociative (mul den) expFrac (recipEConst den) (-int) -- LT -> nest (-int) (divide den e) expFrac approxLogBase :: Integer -> Integer -> (Int, Integer) approxLogBase base x = until ((<=base) . snd) (\(xE,xM) -> (succ xE, div xM base)) (0,x) lnSmall :: Integer -> Integer -> Integer lnSmall den x = evalPowerSeries PSE.log den (x-den) -- uses Double's log for an estimate and dramatic speed up ln :: Integer -> Integer -> Integer ln den x = let fac = 10^50 {- A constant which is representable by Double and which will quickly split our number it pieces small enough for Double. -} (denE, denM) = approxLogBase fac den (xE, xM) = approxLogBase fac x approxDouble :: Double approxDouble = log (NP.fromInteger fac) * fromIntegral (xE-denE) + log (NP.fromInteger xM / NP.fromInteger denM) {- We convert first with respect to @fac@ in order to keep in the range of Double values. -} approxFac = round (approxDouble * NP.fromInteger fac) approx = fromFixedPoint den fac approxFac xSmall = divide den x (exp den approx) in add den approx (lnSmall den xSmall)