{-# 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.RealRing    as RealRing
-- 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 NumericPrelude.Base
import NumericPrelude.Numeric hiding (recip, sqrt, exp, sin, cos, tan,
                              fromRational')

import qualified NumericPrelude.Numeric as NP


{- ** Conversion -}

{- ** other number types -}

fromFloat :: RealRing.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.

ToDo:
The algorithm just computes floor(sqrt(den*x)).
We might factor out the algorithm for (floor.sqrt)
and move it to a different module
together with Fermat factors and so on.
-}
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)