{-# LANGUAGE RebindableSyntax #-}
{- |
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.Transcendental as Trans
import qualified MathObj.PowerSeries.Example as PSE

import qualified Data.List.Reverse.StrictElement as Rev
import NumericPrelude.List (mapLast, )
import Data.Function.HT (powerAssociative, )
import Data.List.HT (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 :: Integer -> a -> Integer
fromFloat Integer
den a
x =
   a -> Integer
forall a b. (C a, C b) => a -> b
round (a
x a -> a -> a
forall a. C a => a -> a -> a
* Integer -> a
forall a. C a => Integer -> a
NP.fromInteger Integer
den)

-- | denominator conversion
fromFixedPoint :: Integer -> Integer -> Integer -> Integer
fromFixedPoint :: Integer -> Integer -> Integer -> Integer
fromFixedPoint Integer
denDst Integer
denSrc Integer
x = Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div (Integer
xInteger -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
denDst) Integer
denSrc


{- ** text -}

{- |
very efficient because it can make use of the decimal output of 'show'
-}
showPositionalDec :: Integer -> Integer -> String
showPositionalDec :: Integer -> Integer -> String
showPositionalDec Integer
den = (Integer -> String) -> Integer -> String
liftShowPosToInt ((Integer -> String) -> Integer -> String)
-> (Integer -> String) -> Integer -> String
forall a b. (a -> b) -> a -> b
$ \Integer
x ->
   let packetSize :: Int
packetSize = Int
50  -- process digits in packets of this size
       basis :: Integer
basis = Int -> Integer -> Integer
forall a b. (C a, C b) => b -> a -> a
ringPower Int
packetSize Integer
10
       (Integer
int,[Integer]
frac) = Integer -> Integer -> Integer -> (Integer, [Integer])
toPositional Integer
basis Integer
den Integer
x
   in  Integer -> String
forall a. Show a => a -> String
show Integer
int String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"." String -> String -> String
forall a. [a] -> [a] -> [a]
++
          [String] -> String
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ((String -> String) -> [String] -> [String]
forall a. (a -> a) -> [a] -> [a]
mapLast ((Char -> Bool) -> String -> String
forall a. (a -> Bool) -> [a] -> [a]
Rev.dropWhile (Char
'0'Char -> Char -> Bool
forall a. Eq a => a -> a -> Bool
==))
             ((Integer -> String) -> [Integer] -> [String]
forall a b. (a -> b) -> [a] -> [b]
map (Char -> Int -> String -> String
forall a. a -> Int -> [a] -> [a]
padLeft Char
'0' Int
packetSize (String -> String) -> (Integer -> String) -> Integer -> String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> String
forall a. Show a => a -> String
show) [Integer]
frac))

showPositionalHex :: Integer -> Integer -> String
showPositionalHex :: Integer -> Integer -> String
showPositionalHex = Integer -> Integer -> Integer -> String
showPositionalBasis Integer
16

showPositionalBin :: Integer -> Integer -> String
showPositionalBin :: Integer -> Integer -> String
showPositionalBin = Integer -> Integer -> Integer -> String
showPositionalBasis Integer
2

showPositionalBasis :: Integer -> Integer -> Integer -> String
showPositionalBasis :: Integer -> Integer -> Integer -> String
showPositionalBasis Integer
basis Integer
den = (Integer -> String) -> Integer -> String
liftShowPosToInt ((Integer -> String) -> Integer -> String)
-> (Integer -> String) -> Integer -> String
forall a b. (a -> b) -> a -> b
$ \Integer
x ->
   let (Integer
int,[Integer]
frac) = Integer -> Integer -> Integer -> (Integer, [Integer])
toPositional Integer
basis Integer
den Integer
x
   in  Integer -> String
forall a. Show a => a -> String
show Integer
int String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"." String -> String -> String
forall a. [a] -> [a] -> [a]
++ (Integer -> Char) -> [Integer] -> String
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Char
intToDigit (Int -> Char) -> (Integer -> Int) -> Integer -> Char
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Int
forall a. C a => Integer -> a
fromInteger) [Integer]
frac

liftShowPosToInt :: (Integer -> String) -> (Integer -> String)
liftShowPosToInt :: (Integer -> String) -> Integer -> String
liftShowPosToInt Integer -> String
f Integer
n =
   if Integer
nInteger -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>=Integer
0
     then       Integer -> String
f   Integer
n
     else Char
'-' Char -> String -> String
forall a. a -> [a] -> [a]
: Integer -> String
f (-Integer
n)

toPositional :: Integer -> Integer -> Integer -> (Integer, [Integer])
toPositional :: Integer -> Integer -> Integer -> (Integer, [Integer])
toPositional Integer
basis Integer
den Integer
x =
   let (Integer
int, Integer
frac) = Integer -> Integer -> (Integer, Integer)
forall a. C a => a -> a -> (a, a)
divMod Integer
x Integer
den
   in  (Integer
int, (Integer -> Maybe (Integer, Integer)) -> Integer -> [Integer]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr (\Integer
rm -> Bool -> (Integer, Integer) -> Maybe (Integer, Integer)
forall a. Bool -> a -> Maybe a
toMaybe (Integer
rmInteger -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/=Integer
0) (Integer -> Integer -> (Integer, Integer)
forall a. C a => a -> a -> (a, a)
divMod (Integer
basisInteger -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
rm) Integer
den)) Integer
frac)


{- * Additive -}

add :: Integer -> Integer -> Integer -> Integer
add :: Integer -> Integer -> Integer -> Integer
add Integer
_ = Integer -> Integer -> Integer
forall a. C a => a -> a -> a
(+)

sub :: Integer -> Integer -> Integer -> Integer
sub :: Integer -> Integer -> Integer -> Integer
sub Integer
_ = (-)


{- * Ring -}

mul :: Integer -> Integer -> Integer -> Integer
mul :: Integer -> Integer -> Integer -> Integer
mul Integer
den Integer
x Integer
y = Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div (Integer
xInteger -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
y) Integer
den


{- * Field -}

divide :: Integer -> Integer -> Integer -> Integer
divide :: Integer -> Integer -> Integer -> Integer
divide Integer
den Integer
x Integer
y = Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div (Integer
xInteger -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
den) Integer
y

recip :: Integer -> Integer -> Integer
recip :: Integer -> Integer -> Integer
recip Integer
den Integer
x = Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div (Integer
denInteger -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
2) Integer
x


{- * Algebra -}

{-
Newton's method for computing roots.
-}

magnitudes :: [Integer]
magnitudes :: [Integer]
magnitudes =
   [[Integer]] -> [Integer]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[Integer]] -> [[Integer]]
forall a. [[a]] -> [[a]]
transpose [(Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
2) Integer
4, (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
2) Integer
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 :: Integer -> Integer -> Integer
sqrt Integer
den Integer
x =
   let xden :: Integer
xden     = Integer
xInteger -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
den
       initial :: Integer
initial  = (Integer, Integer) -> Integer
forall a b. (a, b) -> a
fst ([(Integer, Integer)] -> (Integer, Integer)
forall a. [a] -> a
head (((Integer, Integer) -> Bool)
-> [(Integer, Integer)] -> [(Integer, Integer)]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile ((Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
xden) (Integer -> Bool)
-> ((Integer, Integer) -> Integer) -> (Integer, Integer) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer, Integer) -> Integer
forall a b. (a, b) -> b
snd)
                                ([Integer] -> [Integer] -> [(Integer, Integer)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer]
magnitudes ([Integer] -> [Integer]
forall a. [a] -> [a]
tail ([Integer] -> [Integer]
forall a. [a] -> [a]
tail [Integer]
magnitudes)))))
       approxs :: [Integer]
approxs  = (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (\Integer
y -> Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div (Integer
y Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+ Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div Integer
xden Integer
y) Integer
2) Integer
initial
       isRoot :: Integer -> Bool
isRoot Integer
y = Integer
yInteger -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
2 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
xden Bool -> Bool -> Bool
&& Integer
xden Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< (Integer
yInteger -> Integer -> Integer
forall a. C a => a -> a -> a
+Integer
1)Integer -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
2
   in  [Integer] -> Integer
forall a. [a] -> a
head ((Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Bool -> Bool
not (Bool -> Bool) -> (Integer -> Bool) -> Integer -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Bool
isRoot) [Integer]
approxs)

-- bug: needs too long:  root (12::Int) (fromIntegerBase 10 1000 2)
root :: Integer -> Integer -> Integer -> Integer
root :: Integer -> Integer -> Integer -> Integer
root Integer
n Integer
den Integer
x =
   let n1 :: Integer
n1       = Integer
nInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
1
       xden :: Integer
xden     = Integer
x Integer -> Integer -> Integer
forall a. C a => a -> a -> a
* Integer
denInteger -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
n1
       initial :: Integer
initial  = (Integer, Integer) -> Integer
forall a b. (a, b) -> a
fst ([(Integer, Integer)] -> (Integer, Integer)
forall a. [a] -> a
head (((Integer, Integer) -> Bool)
-> [(Integer, Integer)] -> [(Integer, Integer)]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile ((\Integer
y -> Integer
yInteger -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
xden) (Integer -> Bool)
-> ((Integer, Integer) -> Integer) -> (Integer, Integer) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer, Integer) -> Integer
forall a b. (a, b) -> b
snd)
                                ([Integer] -> [Integer] -> [(Integer, Integer)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer]
magnitudes ([Integer] -> [Integer]
forall a. [a] -> [a]
tail [Integer]
magnitudes))))
       approxs :: [Integer]
approxs  = (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (\Integer
y -> Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div (Integer
n1Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
y Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+ Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div Integer
xden (Integer
yInteger -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
n1)) Integer
n) Integer
initial
       isRoot :: Integer -> Bool
isRoot Integer
y = Integer
yInteger -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
xden Bool -> Bool -> Bool
&& Integer
xden Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< (Integer
yInteger -> Integer -> Integer
forall a. C a => a -> a -> a
+Integer
1)Integer -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
n
   in  [Integer] -> Integer
forall a. [a] -> a
head ((Integer -> Bool) -> [Integer] -> [Integer]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Bool -> Bool
not (Bool -> Bool) -> (Integer -> Bool) -> Integer -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Bool
isRoot) [Integer]
approxs)



{- * Transcendental -}

-- very simple evaluation by power series with lots of rounding errors
evalPowerSeries :: [Rational] -> Integer -> Integer -> Integer
evalPowerSeries :: [Rational] -> Integer -> Integer -> Integer
evalPowerSeries [Rational]
series Integer
den Integer
x =
   let powers :: [Integer]
powers   = (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer -> Integer -> Integer -> Integer
mul Integer
den Integer
x) Integer
den
       summands :: [Integer]
summands = (Rational -> Integer -> Integer)
-> [Rational] -> [Integer] -> [Integer]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Rational
c Integer
p -> Rational -> Integer
forall a b. (C a, C b) => a -> b
round (Rational
c Rational -> Rational -> Rational
forall a. C a => a -> a -> a
* Integer -> Rational
forall a. C a => Integer -> a
fromInteger Integer
p)) [Rational]
series [Integer]
powers
   in  [Integer] -> Integer
forall a. C a => [a] -> a
sum (((Rational, Integer) -> Integer)
-> [(Rational, Integer)] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Rational, Integer) -> Integer
forall a b. (a, b) -> b
snd (((Rational, Integer) -> Bool)
-> [(Rational, Integer)] -> [(Rational, Integer)]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (\(Rational
c,Integer
s) -> Integer
sInteger -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/=Integer
0 Bool -> Bool -> Bool
|| Rational
cRational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
==Rational
0)
                               ([Rational] -> [Integer] -> [(Rational, Integer)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Rational]
series [Integer]
summands)))

cos, sin, tan :: Integer -> Integer -> Integer
cos :: Integer -> Integer -> Integer
cos = [Rational] -> Integer -> Integer -> Integer
evalPowerSeries [Rational]
forall a. C a => [a]
PSE.cos
sin :: Integer -> Integer -> Integer
sin = [Rational] -> Integer -> Integer -> Integer
evalPowerSeries [Rational]
forall a. C a => [a]
PSE.sin
-- tan will suffer from inaccuracies for small cosine
tan :: Integer -> Integer -> Integer
tan Integer
den Integer
x = Integer -> Integer -> Integer -> Integer
divide Integer
den (Integer -> Integer -> Integer
sin Integer
den Integer
x) (Integer -> Integer -> Integer
cos Integer
den Integer
x)

-- it must abs x <= den
arctanSmall :: Integer -> Integer -> Integer
arctanSmall :: Integer -> Integer -> Integer
arctanSmall = [Rational] -> Integer -> Integer -> Integer
evalPowerSeries [Rational]
forall a. C a => [a]
PSE.atan

-- will fail for large inputs
arctan :: Integer -> Integer -> Integer
arctan :: Integer -> Integer -> Integer
arctan Integer
den Integer
x =
   let estimate :: Integer
estimate = Integer -> Double -> Integer
forall a. C a => Integer -> a -> Integer
fromFloat Integer
den
                     (Double -> Double
forall a. C a => a -> a
Trans.atan (Rational -> Double
forall a. C a => Rational -> a
NP.fromRational' (Integer
x Integer -> Integer -> Rational
forall a. C a => a -> a -> T a
% Integer
den)) :: Double)
       tanEst :: Integer
tanEst   = Integer -> Integer -> Integer
tan Integer
den Integer
estimate
       residue :: Integer
residue  = Integer -> Integer -> Integer -> Integer
divide Integer
den (Integer
xInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
tanEst) (Integer
den Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+ Integer -> Integer -> Integer -> Integer
mul Integer
den Integer
x Integer
tanEst)
   in  Integer
estimate Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+ Integer -> Integer -> Integer
arctanSmall Integer
den Integer
residue

piConst :: Integer -> Integer
piConst :: Integer -> Integer
piConst Integer
den =
   let den4 :: Integer
den4 = Integer
4Integer -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
den
       stArcTan :: Integer -> Integer -> Integer
stArcTan Integer
k Integer
x = let d :: Integer
d = Integer
kInteger -> Integer -> Integer
forall a. C a => a -> a -> a
*Integer
den4 in Integer -> Integer -> Integer
arctanSmall Integer
d (Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div Integer
d Integer
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
       (Integer -> Integer -> Integer
stArcTan Integer
44 Integer
57 Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+ Integer -> Integer -> Integer
stArcTan Integer
7 Integer
239 Integer -> Integer -> Integer
forall a. C a => a -> a -> a
- Integer -> Integer -> Integer
stArcTan Integer
12 Integer
682 Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+ Integer -> Integer -> Integer
stArcTan Integer
24 Integer
12943)


expSmall :: Integer -> Integer -> Integer
expSmall :: Integer -> Integer -> Integer
expSmall = [Rational] -> Integer -> Integer -> Integer
evalPowerSeries [Rational]
forall a. C a => [a]
PSE.exp

eConst :: Integer -> Integer
eConst :: Integer -> Integer
eConst Integer
den = Integer -> Integer -> Integer
expSmall Integer
den Integer
den

recipEConst :: Integer -> Integer
recipEConst :: Integer -> Integer
recipEConst Integer
den = Integer -> Integer -> Integer
expSmall Integer
den (-Integer
den)

exp :: Integer -> Integer -> Integer
exp :: Integer -> Integer -> Integer
exp Integer
den Integer
x =
   let den2 :: Integer
den2 = Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div Integer
den Integer
2
       (Integer
int,Integer
frac) = Integer -> Integer -> (Integer, Integer)
forall a. C a => a -> a -> (a, a)
divMod (Integer
x Integer -> Integer -> Integer
forall a. C a => a -> a -> a
+ Integer
den2) Integer
den
       expFrac :: Integer
expFrac = Integer -> Integer -> Integer
expSmall Integer
den (Integer
fracInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
den2)
   in  case Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Integer
int Integer
0 of
          Ordering
EQ -> Integer
expFrac
          Ordering
GT -> (Integer -> Integer -> Integer)
-> Integer -> Integer -> Integer -> Integer
forall a. (a -> a -> a) -> a -> a -> Integer -> a
powerAssociative (Integer -> Integer -> Integer -> Integer
mul Integer
den) Integer
expFrac (Integer -> Integer
eConst      Integer
den)   Integer
int
          Ordering
LT -> (Integer -> Integer -> Integer)
-> Integer -> Integer -> Integer -> Integer
forall a. (a -> a -> a) -> a -> a -> Integer -> a
powerAssociative (Integer -> Integer -> Integer -> Integer
mul Integer
den) Integer
expFrac (Integer -> Integer
recipEConst Integer
den) (-Integer
int)
          -- LT -> nest (-int) (divide den e) expFrac


approxLogBase :: Integer -> Integer -> (Int, Integer)
approxLogBase :: Integer -> Integer -> (Int, Integer)
approxLogBase Integer
base Integer
x =
   ((Int, Integer) -> Bool)
-> ((Int, Integer) -> (Int, Integer))
-> (Int, Integer)
-> (Int, Integer)
forall a. (a -> Bool) -> (a -> a) -> a -> a
until ((Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<=Integer
base) (Integer -> Bool)
-> ((Int, Integer) -> Integer) -> (Int, Integer) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, Integer) -> Integer
forall a b. (a, b) -> b
snd) (\(Int
xE,Integer
xM) -> (Int -> Int
forall a. Enum a => a -> a
succ Int
xE, Integer -> Integer -> Integer
forall a. C a => a -> a -> a
div Integer
xM Integer
base)) (Int
0,Integer
x)

lnSmall :: Integer -> Integer -> Integer
lnSmall :: Integer -> Integer -> Integer
lnSmall Integer
den Integer
x =
   [Rational] -> Integer -> Integer -> Integer
evalPowerSeries [Rational]
forall a. C a => [a]
PSE.log Integer
den (Integer
xInteger -> Integer -> Integer
forall a. C a => a -> a -> a
-Integer
den)

-- uses Double's log for an estimate and dramatic speed up
ln :: Integer -> Integer -> Integer
ln :: Integer -> Integer -> Integer
ln Integer
den Integer
x =
   let fac :: Integer
fac = Integer
10Integer -> Integer -> Integer
forall a. C a => a -> Integer -> a
^Integer
50 {- A constant which is representable by Double
                      and which will quickly split our number it pieces
                      small enough for Double. -}
       (Int
denE, Integer
denM) = Integer -> Integer -> (Int, Integer)
approxLogBase Integer
fac Integer
den
       (Int
xE,   Integer
xM)   = Integer -> Integer -> (Int, Integer)
approxLogBase Integer
fac Integer
x
       approxDouble :: Double
       approxDouble :: Double
approxDouble =
          Double -> Double
forall a. C a => a -> a
log (Integer -> Double
forall a. C a => Integer -> a
NP.fromInteger Integer
fac) Double -> Double -> Double
forall a. C a => a -> a -> a
* Int -> Double
forall a b. (C a, C b) => a -> b
fromIntegral (Int
xEInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
denE) Double -> Double -> Double
forall a. C a => a -> a -> a
+
          Double -> Double
forall a. C a => a -> a
log (Integer -> Double
forall a. C a => Integer -> a
NP.fromInteger Integer
xM Double -> Double -> Double
forall a. C a => a -> a -> a
/ Integer -> Double
forall a. C a => Integer -> a
NP.fromInteger Integer
denM)
       {- We convert first with respect to @fac@
          in order to keep in the range of Double values. -}
       approxFac :: Integer
approxFac = Double -> Integer
forall a b. (C a, C b) => a -> b
round (Double
approxDouble Double -> Double -> Double
forall a. C a => a -> a -> a
* Integer -> Double
forall a. C a => Integer -> a
NP.fromInteger Integer
fac)
       approx :: Integer
approx    = Integer -> Integer -> Integer -> Integer
fromFixedPoint Integer
den Integer
fac Integer
approxFac
       xSmall :: Integer
xSmall    = Integer -> Integer -> Integer -> Integer
divide Integer
den Integer
x (Integer -> Integer -> Integer
exp Integer
den Integer
approx)
   in  Integer -> Integer -> Integer -> Integer
add Integer
den Integer
approx (Integer -> Integer -> Integer
lnSmall Integer
den Integer
xSmall)