-- Modified by Lennart Augustsson to fit into Haskell numerical hierarchy.
--
-- 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 Data.Number.FixedFunctions where
import Prelude hiding (pi, sqrt, tan, atan, exp, log)
import Data.Ratio

approx      :: Rational -> Rational -> Rational
approx :: Rational -> Rational -> Rational
approx Rational
eps Rational
x = forall a. RealFrac a => a -> a -> Rational
approxRational Rational
x Rational
eps

------------------------------------------------------------------
--              Category: Conversion
--      from continued fraction to fraction and vice versa,
--      from Taylor series to continued fraction.
-------------------------------------------------------------------
type CF = [(Rational, Rational)]

fromCF :: CF -> Rational
fromCF :: CF -> Rational
fromCF CF
x =
        --
        -- Convert finite continued fraction to fraction
        -- evaluating from right to left. This is used
        -- mainly for testing in conjunction with "toCF".
        --
        forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (Rational, Rational) -> Rational -> Rational
g Rational
1 CF
x
        where
            g :: (Rational, Rational) -> Rational -> Rational
            g :: (Rational, Rational) -> Rational -> Rational
g (Rational, Rational)
u Rational
v = (forall a b. (a, b) -> a
fst (Rational, Rational)
u) forall a. Num a => a -> a -> a
+ (forall a b. (a, b) -> b
snd (Rational, Rational)
u) forall a. Fractional a => a -> a -> a
/ Rational
v

toCF :: Rational -> CF
toCF :: Rational -> CF
toCF Rational
x =
        --
        -- Convert fraction to finite continued fraction
        --
        forall {a} {a}.
(Integral a, Integral a) =>
Ratio a -> [(Ratio a, Ratio a)] -> [(Ratio a, Ratio a)]
toCF' Rational
x []
        where
            toCF' :: Ratio a -> [(Ratio a, Ratio a)] -> [(Ratio a, Ratio a)]
toCF' Ratio a
u [(Ratio a, Ratio a)]
lst =
                case a
r of
                a
0 -> forall a. [a] -> [a]
reverse (((a
qforall a. Integral a => a -> a -> Ratio a
%a
1),(a
0forall a. Integral a => a -> a -> Ratio a
%a
1))forall a. a -> [a] -> [a]
:[(Ratio a, Ratio a)]
lst)
                a
_ -> Ratio a -> [(Ratio a, Ratio a)] -> [(Ratio a, Ratio a)]
toCF' (a
bforall a. Integral a => a -> a -> Ratio a
%a
r) (((a
qforall a. Integral a => a -> a -> Ratio a
%a
1),(a
1forall a. Integral a => a -> a -> Ratio a
%a
1))forall a. a -> [a] -> [a]
:[(Ratio a, Ratio a)]
lst)
                where
                    a :: a
a = forall a. Ratio a -> a
numerator Ratio a
u
                    b :: a
b = forall a. Ratio a -> a
denominator Ratio a
u
                    (a
q,a
r) = forall a. Integral a => a -> a -> (a, a)
quotRem a
a a
b


approxCF :: Rational -> CF -> Rational
approxCF :: Rational -> CF -> Rational
approxCF Rational
eps [] = Rational
0
approxCF Rational
eps CF
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.
        --
        = Rational
-> CF
-> Rational
-> Rational
-> Rational
-> Rational
-> Rational
-> Int
-> Rational
approxCF' Rational
eps CF
x Rational
0 Rational
1 Rational
1 Rational
q' Rational
p' Int
1
            where
                h :: Rational
h = forall a b. (a, b) -> a
fst (CF
xforall a. [a] -> Int -> a
!!Int
0)
                (Rational
q', Rational
p') = CF
xforall a. [a] -> Int -> a
!!Int
0
                approxCF' :: Rational
-> CF
-> Rational
-> Rational
-> Rational
-> Rational
-> Rational
-> Int
-> Rational
approxCF' Rational
eps CF
x Rational
v2 Rational
v1 Rational
u2 Rational
u1 Rational
a' Int
n
                    | forall a. Num a => a -> a
abs (Rational
1 forall a. Num a => a -> a -> a
- Rational
f1forall a. Fractional a => a -> a -> a
/Rational
f) forall a. Ord a => a -> a -> Bool
< Rational
eps = Rational -> Rational -> Rational
approx Rational
eps Rational
f
                    | Rational
a forall a. Eq a => a -> a -> Bool
== Rational
0    = Rational -> Rational -> Rational
approx Rational
eps Rational
f
                    | Bool
otherwise = Rational
-> CF
-> Rational
-> Rational
-> Rational
-> Rational
-> Rational
-> Int
-> Rational
approxCF' Rational
eps CF
x Rational
v1 Rational
v Rational
u1 Rational
u Rational
a (Int
nforall a. Num a => a -> a -> a
+Int
1)
                    where
                        (Rational
b, Rational
a) = CF
xforall a. [a] -> Int -> a
!!Int
n
                        u :: Rational
u  = Rational
bforall a. Num a => a -> a -> a
*Rational
u1 forall a. Num a => a -> a -> a
+ Rational
a'forall a. Num a => a -> a -> a
*Rational
u2
                        v :: Rational
v  = Rational
bforall a. Num a => a -> a -> a
*Rational
v1 forall a. Num a => a -> a -> a
+ Rational
a'forall a. Num a => a -> a -> a
*Rational
v2
                        f :: Rational
f  = Rational
uforall a. Fractional a => a -> a -> a
/Rational
v
                        f1 :: Rational
f1 = Rational
u1forall a. Fractional a => a -> a -> a
/Rational
v1


-- Type signature determined by GHC.
fromTaylorToCF :: Fractional a => [a] -> a -> [(a, a)]
fromTaylorToCF :: forall a. Fractional a => [a] -> a -> [(a, a)]
fromTaylorToCF [a]
s a
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
        --
        (a, a)
zeroforall a. a -> [a] -> [a]
:(a, a)
oneforall a. a -> [a] -> [a]
:[Int -> (a, a)
higher Int
m | Int
m <- [Int
2..]]
        where
            zero :: (a, a)
zero      = ([a]
sforall a. [a] -> Int -> a
!!Int
0, [a]
sforall a. [a] -> Int -> a
!!Int
1 forall a. Num a => a -> a -> a
* a
x)
            one :: (a, a)
one       = (a
1, -[a]
sforall a. [a] -> Int -> a
!!Int
2forall a. Fractional a => a -> a -> a
/[a]
sforall a. [a] -> Int -> a
!!Int
1 forall a. Num a => a -> a -> a
* a
x)
            higher :: Int -> (a, a)
higher Int
m  = (a
1 forall a. Num a => a -> a -> a
+ [a]
sforall a. [a] -> Int -> a
!!Int
mforall a. Fractional a => a -> a -> a
/[a]
sforall a. [a] -> Int -> a
!!(Int
mforall a. Num a => a -> a -> a
-Int
1) forall a. Num a => a -> a -> a
* a
x, -[a]
sforall a. [a] -> Int -> a
!!(Int
mforall a. Num a => a -> a -> a
+Int
1)forall a. Fractional a => a -> a -> a
/[a]
sforall a. [a] -> Int -> a
!!Int
m forall a. Num a => a -> a -> a
* a
x)


------------------------------------------------------------------
--                Category: Auxiliaries
------------------------------------------------------------------

fac :: Integer -> Integer
fac :: Integer -> Integer
fac = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Enum a => a -> a -> [a]
enumFromTo Integer
1

integerRoot2 :: Integer -> Integer
integerRoot2 :: Integer -> Integer
integerRoot2 Integer
1 = Integer
1
integerRoot2 Integer
x =
        --
        -- Biggest integer m, such that x - m^2 >= 0,
        -- where x is a positive integer
        --
        forall {a}. Integral a => a -> a -> a -> a -> a
integerRoot2' Integer
0 Integer
x (Integer
x forall a. Integral a => a -> a -> a
`div` Integer
2) Integer
x
        where
            integerRoot2' :: a -> a -> a -> a -> a
integerRoot2' a
lo a
hi a
r a
y
                | a
c forall a. Ord a => a -> a -> Bool
> a
y      = a -> a -> a -> a -> a
integerRoot2' a
lo a
r ((a
r forall a. Num a => a -> a -> a
+ a
lo) forall a. Integral a => a -> a -> a
`div` a
2) a
y
                | a
c forall a. Eq a => a -> a -> Bool
== a
y     = a
r
                | Bool
otherwise  =
                    if (a
rforall a. Num a => a -> a -> a
+a
1)forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 forall a. Ord a => a -> a -> Bool
> a
y then
                        a
r
                    else
                        a -> a -> a -> a -> a
integerRoot2' a
r a
hi ((a
r forall a. Num a => a -> a -> a
+ a
hi) forall a. Integral a => a -> a -> a
`div` a
2) a
y
                    where c :: a
c = a
rforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2

-------------------------------------------------------------------
-- Everything below is the instantiation of class Transcendental
-- for type Rational. See also modules Cofra and Numerus.
--
--                Category: Constants
-------------------------------------------------------------------

pi :: Rational -> Rational
pi :: Rational -> Rational
pi Rational
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
        --
        (Rational -> Rational -> Rational
sqrt Rational
eps Rational
d) forall a. Fractional a => a -> a -> a
/ (Rational -> CF -> Rational
approxCF Rational
eps (forall a. Fractional a => [a] -> a -> [(a, a)]
fromTaylorToCF [Rational]
s Rational
x))
        where
            x :: Rational
x = Integer
1forall a. Integral a => a -> a -> Ratio a
%(Integer
640320forall a b. (Num a, Integral b) => a -> b -> a
^Integer
3)::Rational
            s :: [Rational]
s = [((-Integer
1)forall a b. (Num a, Integral b) => a -> b -> a
^Integer
kforall a. Num a => a -> a -> a
*(Integer -> Integer
fac (Integer
6forall a. Num a => a -> a -> a
*Integer
k))forall a. Integral a => a -> a -> Ratio a
%((Integer -> Integer
fac Integer
k)forall a b. (Num a, Integral b) => a -> b -> a
^Integer
3forall a. Num a => a -> a -> a
*(Integer -> Integer
fac (Integer
3forall a. Num a => a -> a -> a
*Integer
k))))forall a. Num a => a -> a -> a
*((Integer
aforall a. Num a => a -> a -> a
*Integer
kforall a. Num a => a -> a -> a
+Integer
b)forall a. Integral a => a -> a -> Ratio a
%Integer
c) | Integer
k<-[Integer
0..]]
            a :: Integer
a = Integer
545140134
            b :: Integer
b = Integer
13591409
            c :: Integer
c = Integer
426880
            d :: Rational
d = Rational
10005

---------------------------------------------------------------------
--                Category: Trigonometry
---------------------------------------------------------------------

tan :: Rational -> Rational -> Rational
tan :: Rational -> Rational -> Rational
tan Rational
eps Rational
0  = Rational
0
tan Rational
eps Rational
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]
        --
        | Rational
x forall a. Ord a => a -> a -> Bool
>= Rational
half_pi'  = Rational -> Rational -> Rational
tan Rational
eps (Rational
x forall a. Num a => a -> a -> a
- ((Integer
1forall a. Num a => a -> a -> a
+Integer
m)forall a. Integral a => a -> a -> Ratio a
%Integer
1)forall a. Num a => a -> a -> a
*Rational
xpi)
        | Rational
x forall a. Ord a => a -> a -> Bool
<= -Rational
half_pi' = Rational -> Rational -> Rational
tan Rational
eps (Rational
x forall a. Num a => a -> a -> a
+ ((Integer
1forall a. Num a => a -> a -> a
+Integer
m)forall a. Integral a => a -> a -> Ratio a
%Integer
1)forall a. Num a => a -> a -> a
*Rational
xpi)
        --- | absx > 1       = 2 * t/(1 - t^2)
        | Bool
otherwise      = Rational -> CF -> Rational
approxCF Rational
eps (forall {a} {a}.
(Integral a, Integral a) =>
Ratio a -> [(Ratio a, Ratio a)]
cf Rational
x)
        where
            absx :: Rational
absx    = forall a. Num a => a -> a
abs Rational
x
            t :: Rational
t       = Rational -> Rational -> Rational
tan Rational
eps (Rational
xforall a. Fractional a => a -> a -> a
/Rational
2)
            m :: Integer
m       = forall a b. (RealFrac a, Integral b) => a -> b
floor ((Rational
absx forall a. Num a => a -> a -> a
- Rational
half_pi)forall a. Fractional a => a -> a -> a
/ Rational
xpi)
            xpi :: Rational
xpi     = Rational -> Rational
pi Rational
eps
            half_pi' :: Rational
half_pi'= Integer
158forall a. Integral a => a -> a -> Ratio a
%Integer
100
            half_pi :: Rational
half_pi = Rational
xpi forall a. Num a => a -> a -> a
* (Integer
1forall a. Integral a => a -> a -> Ratio a
%Integer
2)
            cf :: Ratio a -> [(Ratio a, Ratio a)]
cf Ratio a
u    = ((a
0forall a. Integral a => a -> a -> Ratio a
%a
1,a
1forall a. Integral a => a -> a -> Ratio a
%a
1)forall a. a -> [a] -> [a]
:[((Ratio a
2forall a. Num a => a -> a -> a
*Ratio a
r forall a. Num a => a -> a -> a
+ Ratio a
1)forall a. Fractional a => a -> a -> a
/Ratio a
u, -Ratio a
1) | Ratio a
r <- [Ratio a
0..]])

sin :: Rational -> Rational -> Rational
sin :: Rational -> Rational -> Rational
sin Rational
eps Rational
0      = Rational
0
sin Rational
eps Rational
x      = Rational
2forall a. Num a => a -> a -> a
*Rational
tforall a. Fractional a => a -> a -> a
/(Rational
1 forall a. Num a => a -> a -> a
+ Rational
tforall a. Num a => a -> a -> a
*Rational
t)
        where
            t :: Rational
t = Rational -> Rational -> Rational
tan Rational
eps (Rational
xforall a. Fractional a => a -> a -> a
/Rational
2)

cos :: Rational -> Rational -> Rational
cos :: Rational -> Rational -> Rational
cos Rational
eps Rational
0      = Rational
1
cos Rational
eps Rational
x      = (Rational
1 forall a. Num a => a -> a -> a
- Rational
p)forall a. Fractional a => a -> a -> a
/(Rational
1 forall a. Num a => a -> a -> a
+ Rational
p)
        where
            t :: Rational
t = Rational -> Rational -> Rational
tan Rational
eps (Rational
xforall a. Fractional a => a -> a -> a
/Rational
2)
            p :: Rational
p = Rational
tforall a. Num a => a -> a -> a
*Rational
t

atan :: Rational -> Rational -> Rational
atan :: Rational -> Rational -> Rational
atan Rational
eps Rational
x
        --
        -- Inverse tangent of x with approximation eps
        --
        | Rational
x forall a. Eq a => a -> a -> Bool
== Rational
0       = Rational
0
        | Rational
x forall a. Ord a => a -> a -> Bool
> Rational
1        =  (Rational -> Rational
pi Rational
eps)forall a. Fractional a => a -> a -> a
/Rational
2 forall a. Num a => a -> a -> a
- Rational -> Rational -> Rational
atan Rational
eps (Rational
1forall a. Fractional a => a -> a -> a
/Rational
x)
        | Rational
x forall a. Ord a => a -> a -> Bool
< -Rational
1       = -(Rational -> Rational
pi Rational
eps)forall a. Fractional a => a -> a -> a
/Rational
2 forall a. Num a => a -> a -> a
- Rational -> Rational -> Rational
atan Rational
eps (Rational
1forall a. Fractional a => a -> a -> a
/Rational
x)
        | Bool
otherwise    = Rational -> CF -> Rational
approxCF Rational
eps ((Rational
0,Rational
x)forall a. a -> [a] -> [a]
:[((Rational
2forall a. Num a => a -> a -> a
*Rational
m forall a. Num a => a -> a -> a
- Rational
1),(Rational
mforall a. Num a => a -> a -> a
*Rational
x)forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2) | Rational
m<- [Rational
1..]])


asin :: Rational -> Rational -> Rational
asin :: Rational -> Rational -> Rational
asin Rational
eps Rational
x
        --
        -- Inverse sine of x with approximation eps
        --
        | Rational
x forall a. Eq a => a -> a -> Bool
== Rational
0    = Rational
0
        | forall a. Num a => a -> a
abs Rational
x forall a. Ord a => a -> a -> Bool
> Rational
1 = forall a. HasCallStack => [Char] -> a
error [Char]
"Fraction.asin"
        | Rational
x forall a. Eq a => a -> a -> Bool
== Rational
1    = (Rational -> Rational
pi Rational
eps) forall a. Num a => a -> a -> a
*  (Integer
1forall a. Integral a => a -> a -> Ratio a
%Integer
2)
        | Rational
x forall a. Eq a => a -> a -> Bool
== -Rational
1   = (Rational -> Rational
pi Rational
eps) forall a. Num a => a -> a -> a
* (-Integer
1forall a. Integral a => a -> a -> Ratio a
%Integer
2)
        | Bool
otherwise = Rational -> Rational -> Rational
atan Rational
eps (Rational
x forall a. Fractional a => a -> a -> a
/ (Rational -> Rational -> Rational
sqrt Rational
eps (Rational
1 forall a. Num a => a -> a -> a
- Rational
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)))


acos :: Rational -> Rational -> Rational
acos :: Rational -> Rational -> Rational
acos Rational
eps Rational
x
        --
        -- Inverse cosine of x with approximation eps
        --
        | Rational
x forall a. Eq a => a -> a -> Bool
== Rational
0    = (Rational -> Rational
pi Rational
eps)forall a. Num a => a -> a -> a
*(Integer
1forall a. Integral a => a -> a -> Ratio a
%Integer
2)
        | forall a. Num a => a -> a
abs Rational
x forall a. Ord a => a -> a -> Bool
> Rational
1 = forall a. HasCallStack => [Char] -> a
error [Char]
"Fraction.sin"
        | Rational
x forall a. Eq a => a -> a -> Bool
== Rational
1    = Rational
0
        | Rational
x forall a. Eq a => a -> a -> Bool
== -Rational
1   = Rational -> Rational
pi Rational
eps
        | Bool
otherwise = Rational -> Rational -> Rational
atan Rational
eps ((Rational -> Rational -> Rational
sqrt Rational
eps (Rational
1 forall a. Num a => a -> a -> a
- Rational
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)) forall a. Fractional a => a -> a -> a
/ Rational
x)

---------------------------------------------------------------------
--                Category: Roots
---------------------------------------------------------------------

sqrt :: Rational -> Rational -> Rational
sqrt :: Rational -> Rational -> Rational
sqrt Rational
eps Rational
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
        --
        | Rational
x forall a. Ord a => a -> a -> Bool
< Rational
0        = forall a. HasCallStack => [Char] -> a
error [Char]
"Fraction.sqrt"
        | Rational
x forall a. Eq a => a -> a -> Bool
== Rational
0       = Rational
0
        | Rational
x forall a. Ord a => a -> a -> Bool
< Rational
1        = Rational
1forall a. Fractional a => a -> a -> a
/(Rational -> Rational -> Rational
sqrt Rational
eps (Rational
1forall a. Fractional a => a -> a -> a
/Rational
x))
        | Bool
otherwise    = Rational -> CF -> Rational
approxCF Rational
eps ((Rational
m,Rational
xforall a. Num a => a -> a -> a
-Rational
mforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)forall a. a -> [a] -> [a]
:[(Rational
2forall a. Num a => a -> a -> a
*Rational
m,Rational
xforall a. Num a => a -> a -> a
-Rational
mforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2) | Integer
r<-[Integer
0..]])
        where
            m :: Rational
m = (Integer -> Integer
integerRoot2 (forall a b. (RealFrac a, Integral b) => a -> b
floor Rational
x))forall a. Integral a => a -> a -> Ratio a
%Integer
1

---------------------------------------------------------------------
--              Category: Exponentials and hyperbolics
---------------------------------------------------------------------

exp :: Rational -> Rational -> Rational
exp :: Rational -> Rational -> Rational
exp Rational
eps Rational
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
        --
        | Rational
x forall a. Eq a => a -> a -> Bool
== Rational
0       = Rational
1
        | Rational
x forall a. Ord a => a -> a -> Bool
> Rational
1        = (Rational -> CF -> Rational
approxCF Rational
eps (forall {a}. (Fractional a, Enum a) => a -> [(a, a)]
f (Rational
xforall a. Num a => a -> a -> a
*(Integer
1forall a. Integral a => a -> a -> Ratio a
%Integer
p))))forall a b. (Num a, Integral b) => a -> b -> a
^Integer
p
        | Rational
x forall a. Ord a => a -> a -> Bool
< (-Rational
1)     = (Rational -> CF -> Rational
approxCF Rational
eps (forall {a}. (Fractional a, Enum a) => a -> [(a, a)]
f (Rational
xforall a. Num a => a -> a -> a
*(Integer
1forall a. Integral a => a -> a -> Ratio a
%Integer
q))))forall a b. (Num a, Integral b) => a -> b -> a
^Integer
q
        | Bool
otherwise    = Rational -> CF -> Rational
approxCF Rational
eps (forall {a}. (Fractional a, Enum a) => a -> [(a, a)]
f Rational
x)
        where
            p :: Integer
p = forall a b. (RealFrac a, Integral b) => a -> b
ceiling Rational
x
            q :: Integer
q = -(forall a b. (RealFrac a, Integral b) => a -> b
floor Rational
x)
            f :: a -> [(a, a)]
f a
y = (a
1,a
y)forall a. a -> [a] -> [a]
:(a
1forall a. Num a => a -> a -> a
-a
yforall a. Fractional a => a -> a -> a
/a
2,a
yforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2forall a. Fractional a => a -> a -> a
/a
12)forall a. a -> [a] -> [a]
:[(a
1,a
yforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2forall a. Fractional a => a -> a -> a
/(a
16forall a. Num a => a -> a -> a
*a
nforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2forall a. Num a => a -> a -> a
-a
4)) | a
n<-[a
2..]]


cosh :: Rational -> Rational -> Rational
cosh :: Rational -> Rational -> Rational
cosh Rational
eps Rational
x =
        --
        -- Hyperbolic cosine with approximation eps
        --
        (Rational
a forall a. Num a => a -> a -> a
+ Rational
b)forall a. Num a => a -> a -> a
*(Integer
1forall a. Integral a => a -> a -> Ratio a
%Integer
2)
        where
            a :: Rational
a = Rational -> Rational -> Rational
exp Rational
eps Rational
x
            b :: Rational
b = Rational
1forall a. Fractional a => a -> a -> a
/Rational
a

sinh :: Rational -> Rational -> Rational
sinh :: Rational -> Rational -> Rational
sinh Rational
eps Rational
x =
        --
        -- Hyperbolic sine with approximation eps
        --
        (Rational
a forall a. Num a => a -> a -> a
- Rational
b)forall a. Num a => a -> a -> a
*(Integer
1forall a. Integral a => a -> a -> Ratio a
%Integer
2)
        where
            a :: Rational
a = Rational -> Rational -> Rational
exp Rational
eps Rational
x
            b :: Rational
b = Rational
1forall a. Fractional a => a -> a -> a
/Rational
a

tanh :: Rational -> Rational -> Rational
tanh :: Rational -> Rational -> Rational
tanh Rational
eps Rational
x =
        --
        -- Hyperbolic tangent with approximation eps
        --
        (Rational
a forall a. Num a => a -> a -> a
- Rational
b)forall a. Fractional a => a -> a -> a
/ (Rational
a forall a. Num a => a -> a -> a
+ Rational
b)
        where
            a :: Rational
a = Rational -> Rational -> Rational
exp Rational
eps Rational
x
            b :: Rational
b = Rational
1forall a. Fractional a => a -> a -> a
/Rational
a

atanh :: Rational -> Rational -> Rational
atanh :: Rational -> Rational -> Rational
atanh Rational
eps Rational
x
        --
        -- Inverse hyperbolic tangent with approximation eps
        --

--      | x >= 1     = 1%0
--      | x <= -1    = -1%0
        | Bool
otherwise  = (Integer
1forall a. Integral a => a -> a -> Ratio a
%Integer
2) forall a. Num a => a -> a -> a
* (Rational -> Rational -> Rational
log Rational
eps ((Rational
1 forall a. Num a => a -> a -> a
+ Rational
x) forall a. Fractional a => a -> a -> a
/ (Rational
1 forall a. Num a => a -> a -> a
- Rational
x)))

asinh :: Rational -> Rational -> Rational
asinh :: Rational -> Rational -> Rational
asinh Rational
eps Rational
x
        --
        -- Inverse hyperbolic sine
        --
--      | x == 1%0  =  1%0
--      | x == -1%0 = -1%0
        | Bool
otherwise  = Rational -> Rational -> Rational
log Rational
eps (Rational
x forall a. Num a => a -> a -> a
+ (Rational -> Rational -> Rational
sqrt Rational
eps (Rational
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 forall a. Num a => a -> a -> a
+ Rational
1)))

acosh :: Rational -> Rational -> Rational
acosh :: Rational -> Rational -> Rational
acosh Rational
eps Rational
x
        --
        -- Inverse hyperbolic cosine
        --
--      | x == 1%0 = 1%0
--      | x < 1     = 1%0
        | Bool
otherwise = Rational -> Rational -> Rational
log Rational
eps (Rational
x forall a. Num a => a -> a -> a
+ (Rational -> Rational -> Rational
sqrt Rational
eps (Rational
xforall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 forall a. Num a => a -> a -> a
- Rational
1)))

---------------------------------------------------------------------
--                Category: Logarithms
---------------------------------------------------------------------

log :: Rational -> Rational -> Rational
log :: Rational -> Rational -> Rational
log Rational
eps Rational
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)
        --
        | Rational
x forall a. Ord a => a -> a -> Bool
<= Rational
0    = forall a. HasCallStack => [Char] -> a
error [Char]
"Fraction.log"
        | Rational
x forall a. Ord a => a -> a -> Bool
<  Rational
1    = -Rational -> Rational -> Rational
log Rational
eps (Rational
1forall a. Fractional a => a -> a -> a
/Rational
x)
        | Rational
x forall a. Eq a => a -> a -> Bool
== Rational
1    =  Rational
0
        | Bool
otherwise =
            case ((Rational, Integer) -> (Rational, Integer)
scaled (Rational
x,Integer
0)) of
            (Rational
1,Integer
s) -> (Integer
sforall a. Integral a => a -> a -> Ratio a
%Integer
1) forall a. Num a => a -> a -> a
* Rational -> CF -> Rational
approxCF Rational
eps (Rational -> CF
series Rational
1)
            (Rational
y,Integer
0) -> Rational -> CF -> Rational
approxCF Rational
eps (Rational -> CF
series (Rational
yforall a. Num a => a -> a -> a
-Rational
1))
            (Rational
y,Integer
s) -> Rational -> CF -> Rational
approxCF Rational
eps (Rational -> CF
series (Rational
yforall a. Num a => a -> a -> a
-Rational
1)) forall a. Num a => a -> a -> a
+ (Integer
sforall a. Integral a => a -> a -> Ratio a
%Integer
1)forall a. Num a => a -> a -> a
*Rational -> CF -> Rational
approxCF Rational
eps (Rational -> CF
series Rational
1)
        where
            series :: Rational -> CF
            series :: Rational -> CF
series Rational
u = (Rational
0,Rational
u)forall a. a -> [a] -> [a]
:(Rational
1,Rational
uforall a. Fractional a => a -> a -> a
/Rational
2)forall a. a -> [a] -> [a]
:[(Rational
1,Rational
uforall a. Num a => a -> a -> a
*((Integer
mforall a. Num a => a -> a -> a
+Integer
n)forall a. Integral a => a -> a -> Ratio a
%(Integer
4forall a. Num a => a -> a -> a
*Integer
m forall a. Num a => a -> a -> a
+ Integer
2)))|Integer
m<-[Integer
1..],Integer
n<-[Integer
0,Integer
1]]
            scaled :: (Rational,Integer) -> (Rational, Integer)
            scaled :: (Rational, Integer) -> (Rational, Integer)
scaled (Rational
x, Integer
n)
                | Rational
x forall a. Eq a => a -> a -> Bool
== Rational
2 = (Rational
1,Integer
nforall a. Num a => a -> a -> a
+Integer
1)
                | Rational
x forall a. Ord a => a -> a -> Bool
< Rational
2 = (Rational
x, Integer
n)
                | Bool
otherwise = (Rational, Integer) -> (Rational, Integer)
scaled (Rational
xforall a. Num a => a -> a -> a
*(Integer
1forall a. Integral a => a -> a -> Ratio a
%Integer
2), Integer
nforall a. Num a => a -> a -> a
+Integer
1)


---------------------------------------------------------------------------
-- 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 Rationals.
--
--      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
--
-----------------------------------------------------------------------------