module Number.FixedPoint where
import qualified Algebra.RealField as RealField
import qualified Algebra.Additive as Additive
import qualified Algebra.Transcendental as Trans
import qualified MathObj.PowerSeries.Example as PSE
import NumericPrelude.List (dropWhileRev, mapLast, padLeft)
import NumericPrelude.Condition (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
fromFloat :: RealField.C a => Integer -> a -> Integer
fromFloat den x =
round (x * NP.fromInteger den)
fromFixedPoint :: Integer -> Integer -> Integer -> Integer
fromFixedPoint denDst denSrc x = div (x*denDst) denSrc
showPositionalDec :: Integer -> Integer -> String
showPositionalDec den = liftShowPosToInt $ \x ->
let packetSize = 50
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)
add :: Integer -> Integer -> Integer -> Integer
add _ = (+)
sub :: Integer -> Integer -> Integer -> Integer
sub _ = ()
mul :: Integer -> Integer -> Integer -> Integer
mul den x y = div (x*y) den
divide :: Integer -> Integer -> Integer -> Integer
divide den x y = div (x*den) y
recip :: Integer -> Integer -> Integer
recip den x = div (den^2) x
magnitudes :: [Integer]
magnitudes =
concat (transpose [iterate (^2) 4, iterate (^2) 8])
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)
root :: Integer -> Integer -> Integer -> Integer
root n den x =
let n1 = n1
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)
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 den x = divide den (sin den x) (cos den x)
arctanSmall :: Integer -> Integer -> Integer
arctanSmall = evalPowerSeries PSE.atan
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 (xtanEst) (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
(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 (fracden2)
in case compare int 0 of
EQ -> expFrac
GT -> reduceRepeated (mul den) expFrac (eConst den) int
LT -> reduceRepeated (mul den) expFrac (recipEConst den) (int)
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 (xden)
ln :: Integer -> Integer -> Integer
ln den x =
let fac = 10^50
(denE, denM) = approxLogBase fac den
(xE, xM) = approxLogBase fac x
approxDouble :: Double
approxDouble =
log (NP.fromInteger fac) * fromIntegral (xEdenE) +
log (NP.fromInteger xM / NP.fromInteger denM)
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)