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

toCF :: Rational -> CF
toCF :: Rational -> CF
toCF x :: Rational
x =
        --
        -- Convert fraction to finite continued fraction
        --
        Rational -> CF -> CF
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' u :: Ratio a
u lst :: [(Ratio a, Ratio a)]
lst =
                case a
r of
                0 -> [(Ratio a, Ratio a)] -> [(Ratio a, Ratio a)]
forall a. [a] -> [a]
reverse (((a
qa -> a -> Ratio a
forall a. Integral a => a -> a -> Ratio a
%1),(0a -> a -> Ratio a
forall a. Integral a => a -> a -> Ratio a
%1))(Ratio a, Ratio a) -> [(Ratio a, Ratio a)] -> [(Ratio a, Ratio a)]
forall a. a -> [a] -> [a]
:[(Ratio a, Ratio a)]
lst)
                _ -> Ratio a -> [(Ratio a, Ratio a)] -> [(Ratio a, Ratio a)]
toCF' (a
ba -> a -> Ratio a
forall a. Integral a => a -> a -> Ratio a
%a
r) (((a
qa -> a -> Ratio a
forall a. Integral a => a -> a -> Ratio a
%1),(1a -> a -> Ratio a
forall a. Integral a => a -> a -> Ratio a
%1))(Ratio a, Ratio a) -> [(Ratio a, Ratio a)] -> [(Ratio a, Ratio a)]
forall a. a -> [a] -> [a]
:[(Ratio a, Ratio a)]
lst)
                where
                    a :: a
a = Ratio a -> a
forall a. Ratio a -> a
numerator Ratio a
u
                    b :: a
b = Ratio a -> a
forall a. Ratio a -> a
denominator Ratio a
u
                    (q :: a
q,r :: a
r) = a -> a -> (a, a)
forall a. Integral a => a -> a -> (a, a)
quotRem a
a a
b


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


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


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

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

integerRoot2 :: Integer -> Integer
integerRoot2 :: Integer -> Integer
integerRoot2 1 = 1
integerRoot2 x :: Integer
x =
        --
        -- Biggest integer m, such that x - m^2 >= 0,
        -- where x is a positive integer
        --
        Integer -> Integer -> Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a -> a -> a
integerRoot2' 0 Integer
x (Integer
x Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` 2) Integer
x
        where
            integerRoot2' :: a -> a -> a -> a -> a
integerRoot2' lo :: a
lo hi :: a
hi r :: a
r y :: a
y
                | a
c a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
y      = a -> a -> a -> a -> a
integerRoot2' a
lo a
r ((a
r a -> a -> a
forall a. Num a => a -> a -> a
+ a
lo) a -> a -> a
forall a. Integral a => a -> a -> a
`div` 2) a
y
                | a
c a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
y     = a
r
                | Bool
otherwise  =
                    if (a
ra -> a -> a
forall a. Num a => a -> a -> a
+1)a -> Integer -> a
forall a b. (Num a, Integral b) => a -> b -> a
^2 a -> a -> Bool
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 a -> a -> a
forall a. Num a => a -> a -> a
+ a
hi) a -> a -> a
forall a. Integral a => a -> a -> a
`div` 2) a
y
                    where c :: a
c = a
ra -> Integer -> a
forall a b. (Num a, Integral b) => a -> b -> a
^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 eps :: 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) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ (Rational -> CF -> Rational
approxCF Rational
eps ([Rational] -> Rational -> CF
forall a. Fractional a => [a] -> a -> [(a, a)]
fromTaylorToCF [Rational]
s Rational
x))
        where
            x :: Rational
x = 1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%(640320Integer -> Integer -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^3)::Rational
            s :: [Rational]
s = [((-1)Integer -> Integer -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
kInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer -> Integer
fac (6Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
k))Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%((Integer -> Integer
fac Integer
k)Integer -> Integer -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^3Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer -> Integer
fac (3Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
k))))Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
*((Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
kInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
b)Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
c) | Integer
k<-[0..]]
            a :: Integer
a = 545140134
            b :: Integer
b = 13591409
            c :: Integer
c = 426880
            d :: Rational
d = 10005

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

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

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

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

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


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


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

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

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

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

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


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

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

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

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

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

asinh :: Rational -> Rational -> Rational
asinh :: Rational -> Rational -> Rational
asinh eps :: Rational
eps x :: 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 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ (Rational -> Rational -> Rational
sqrt Rational
eps (Rational
xRational -> Integer -> Rational
forall a b. (Num a, Integral b) => a -> b -> a
^2 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ 1)))

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

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

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