{-# LANGUAGE ParallelListComp #-}
{-# LANGUAGE CPP #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
{-# LANGUAGE Safe #-}
#endif
module Math.ContinuedFraction
( CF
, cf, gcf
, asCF, asGCF
, truncateCF
, equiv
, setNumerators
, setDenominators
, partitionCF
, evenCF
, oddCF
, convergents
, steed
, lentz, lentzWith
, modifiedLentz, modifiedLentzWith
, sumPartialProducts
) where
import Control.Arrow ((***))
import Data.List (tails, mapAccumL)
data CF a
= CF a [a]
| GCF a [(a,a)]
cf :: a -> [a] -> CF a
cf :: a -> [a] -> CF a
cf = a -> [a] -> CF a
forall a. a -> [a] -> CF a
CF
gcf :: a -> [(a,a)] -> CF a
gcf :: a -> [(a, a)] -> CF a
gcf = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
GCF
instance Show a => Show (CF a) where
showsPrec :: Int -> CF a -> ShowS
showsPrec Int
p (CF a
b0 [a]
ab) = Bool -> ShowS -> ShowS
showParen (Int
pInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
10)
( String -> ShowS
showString String
"cf "
ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> a -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 a
b0
ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar Char
' '
ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [a] -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 [a]
ab
)
showsPrec Int
p (GCF a
b0 [(a, a)]
ab) = Bool -> ShowS -> ShowS
showParen (Int
pInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
10)
( String -> ShowS
showString String
"gcf "
ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> a -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 a
b0
ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar Char
' '
ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [(a, a)] -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 [(a, a)]
ab
)
instance Functor CF where
fmap :: (a -> b) -> CF a -> CF b
fmap a -> b
f (CF a
b0 [a]
cf) = b -> [b] -> CF b
forall a. a -> [a] -> CF a
CF (a -> b
f a
b0) ((a -> b) -> [a] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map a -> b
f [a]
cf)
fmap a -> b
f (GCF a
b0 [(a, a)]
gcf) = b -> [(b, b)] -> CF b
forall a. a -> [(a, a)] -> CF a
GCF (a -> b
f a
b0) (((a, a) -> (b, b)) -> [(a, a)] -> [(b, b)]
forall a b. (a -> b) -> [a] -> [b]
map (a -> b
f (a -> b) -> (a -> b) -> (a, a) -> (b, b)
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** a -> b
f) [(a, a)]
gcf)
asCF :: Fractional a => CF a -> (a, [a])
asCF :: CF a -> (a, [a])
asCF (CF a
b0 [a]
cf) = (a
b0, [a]
cf)
asCF (GCF a
b0 []) = (a
b0, [])
asCF (GCF a
b0 [(a, a)]
cf) = (a
b0, (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
bs [a]
cs)
where
(a
a:[a]
as, [a]
bs) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
cf
cs :: [a]
cs = a -> a
forall a. Fractional a => a -> a
recip a
a a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a -> a
forall a. Fractional a => a -> a
recip (a
aa -> a -> a
forall a. Num a => a -> a -> a
*a
c) | a
c <- [a]
cs | a
a <- [a]
as]
asGCF :: (Num a, Eq a) => CF a -> (a,[(a,a)])
asGCF :: CF a -> (a, [(a, a)])
asGCF (CF a
b0 [a]
cf) = (a
b0, [(a
1, a
b) | a
b <- [a]
cf])
asGCF (GCF a
b0 [(a, a)]
gcf) = (a
b0, ((a, a) -> Bool) -> [(a, a)] -> [(a, a)]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile ((a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=a
0)(a -> Bool) -> ((a, a) -> a) -> (a, a) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(a, a) -> a
forall a b. (a, b) -> a
fst) [(a, a)]
gcf)
truncateCF :: Int -> CF a -> CF a
truncateCF :: Int -> CF a -> CF a
truncateCF Int
n (CF a
b0 [a]
ab) = a -> [a] -> CF a
forall a. a -> [a] -> CF a
CF a
b0 (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n [a]
ab)
truncateCF Int
n (GCF a
b0 [(a, a)]
ab) = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
GCF a
b0 (Int -> [(a, a)] -> [(a, a)]
forall a. Int -> [a] -> [a]
take Int
n [(a, a)]
ab)
equiv :: (Num a, Eq a) => [a] -> CF a -> CF a
equiv :: [a] -> CF a -> CF a
equiv [a]
cs CF a
orig
= a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf a
b0 ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
as' [a]
bs')
where
(a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
([a]
as,[a]
bs) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
terms
as' :: [a]
as' = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
cs' (a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
cs')) [a]
as
bs' :: [a]
bs' = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
cs' [a]
bs
cs' :: [a]
cs' = [a]
cs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
1
setDenominators :: (Fractional a, Eq a) => [a] -> CF a -> CF a
setDenominators :: [a] -> CF a -> CF a
setDenominators [a]
denoms CF a
orig
= a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf a
b0 ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
as' [a]
bs')
where
(a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
([a]
as,[a]
bs) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
terms
as' :: [a]
as' = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
as ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
cs (a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
cs))
bs' :: [a]
bs' = ((a -> a) -> a -> a) -> [a -> a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
($) ((a -> a -> a) -> [a] -> [a -> a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a -> a
forall a b. a -> b -> a
const [a]
denoms [a -> a] -> [a -> a] -> [a -> a]
forall a. [a] -> [a] -> [a]
++ (a -> a) -> [a -> a]
forall a. a -> [a]
repeat a -> a
forall a. a -> a
id) [a]
bs
cs :: [a]
cs = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
bs' [a]
bs
setNumerators :: (Fractional a, Eq a) => [a] -> CF a -> CF a
setNumerators :: [a] -> CF a -> CF a
setNumerators [a]
numers CF a
orig
= a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf a
b0 ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
as' [a]
bs')
where
(a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
([a]
as,[a]
bs) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
terms
as' :: [a]
as' = ((a -> a) -> a -> a) -> [a -> a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
($) ((a -> a -> a) -> [a] -> [a -> a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a -> a
forall a b. a -> b -> a
const [a]
numers [a -> a] -> [a -> a] -> [a -> a]
forall a. [a] -> [a] -> [a]
++ (a -> a) -> [a -> a]
forall a. a -> [a]
repeat a -> a
forall a. a -> a
id) [a]
as
bs' :: [a]
bs' = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
bs [a]
cs
cs :: [a]
cs = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
as' ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
as (a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
cs))
partitionCF :: (Fractional a, Eq a) => CF a -> (CF a, CF a)
partitionCF :: CF a -> (CF a, CF a)
partitionCF CF a
orig = case [(a, a)]
terms of
[] -> (CF a
orig, CF a
orig)
[(a
a1,a
b1)] ->
let final :: CF a
final = a -> [a] -> CF a
forall a. a -> [a] -> CF a
cf (a
b0 a -> a -> a
forall a. Num a => a -> a -> a
+ a
a1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
b1) []
in (CF a
final, CF a
final)
[(a, a)]
_ -> (CF a
evenPart, CF a
oddPart)
where
(a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
([a]
as, [a]
bs) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
terms
pairs :: [b] -> [(b, b)]
pairs (b
a:b
b:[b]
rest) = (b
a,b
b) (b, b) -> [(b, b)] -> [(b, b)]
forall a. a -> [a] -> [a]
: [b] -> [(b, b)]
pairs [b]
rest
pairs [b
a] = [(b
a,b
0)]
pairs [b]
_ = []
alphas :: [a]
alphas@(a
alpha1:a
alpha2:[a]
_) = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
as ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
bs (a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
bs))
evenPart :: CF a
evenPart = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf a
b0 ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
cs [a]
ds)
where
cs :: [a]
cs = a
alpha1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [(-a
aOdd) a -> a -> a
forall a. Num a => a -> a -> a
* a
aEven | (a
aEven, a
aOdd) <- [a] -> [(a, a)]
forall b. Num b => [b] -> [(b, b)]
pairs ([a] -> [a]
forall a. [a] -> [a]
tail [a]
alphas)]
ds :: [a]
ds = a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
alpha2 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
aOdd a -> a -> a
forall a. Num a => a -> a -> a
+ a
aEven | (a
aOdd, a
aEven) <- [(a, a)] -> [(a, a)]
forall a. [a] -> [a]
tail ([a] -> [(a, a)]
forall b. Num b => [b] -> [(b, b)]
pairs [a]
alphas)]
oddPart :: CF a
oddPart = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf (a
b0a -> a -> a
forall a. Num a => a -> a -> a
+a
alpha1) ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
cs [a]
ds)
where
cs :: [a]
cs = [(-a
aOdd) a -> a -> a
forall a. Num a => a -> a -> a
* a
aEven | (a
aOdd, a
aEven) <- [a] -> [(a, a)]
forall b. Num b => [b] -> [(b, b)]
pairs [a]
alphas]
ds :: [a]
ds = [a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
aOdd a -> a -> a
forall a. Num a => a -> a -> a
+ a
aEven | (a
aEven, a
aOdd) <- [a] -> [(a, a)]
forall b. Num b => [b] -> [(b, b)]
pairs ([a] -> [a]
forall a. [a] -> [a]
tail [a]
alphas)]
evenCF :: (Fractional a, Eq a) => CF a -> CF a
evenCF :: CF a -> CF a
evenCF = (CF a, CF a) -> CF a
forall a b. (a, b) -> a
fst ((CF a, CF a) -> CF a) -> (CF a -> (CF a, CF a)) -> CF a -> CF a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CF a -> (CF a, CF a)
forall a. (Fractional a, Eq a) => CF a -> (CF a, CF a)
partitionCF
oddCF :: (Fractional a, Eq a) => CF a -> CF a
oddCF :: CF a -> CF a
oddCF = (CF a, CF a) -> CF a
forall a b. (a, b) -> b
snd ((CF a, CF a) -> CF a) -> (CF a -> (CF a, CF a)) -> CF a -> CF a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CF a -> (CF a, CF a)
forall a. (Fractional a, Eq a) => CF a -> (CF a, CF a)
partitionCF
convergents :: (Fractional a, Eq a) => CF a -> [a]
convergents :: CF a -> [a]
convergents CF a
orig = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1 ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
nums [a]
denoms)
where
(a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
nums :: [a]
nums = a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:a
b0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
x1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
x0 | a
x0:a
x1:[a]
_ <- [a] -> [[a]]
forall a. [a] -> [[a]]
tails [a]
nums | (a
a,a
b) <- [(a, a)]
terms]
denoms :: [a]
denoms = a
0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:a
1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
x1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
x0 | a
x0:a
x1:[a]
_ <- [a] -> [[a]]
forall a. [a] -> [[a]]
tails [a]
denoms | (a
a,a
b) <- [(a, a)]
terms]
steed :: (Fractional a, Eq a) => CF a -> [a]
steed :: CF a -> [a]
steed (CF a
b0 []) = [a
b0]
steed (GCF a
b0 []) = [a
b0]
steed (CF a
0 ( a
a :[a]
rest)) = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/) (CF a -> [a]
forall a. (Fractional a, Eq a) => CF a -> [a]
steed (a -> [a] -> CF a
forall a. a -> [a] -> CF a
CF a
a [a]
rest))
steed (GCF a
0 ((a
a,a
b):[(a, a)]
rest)) = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
a a -> a -> a
forall a. Fractional a => a -> a -> a
/) (CF a -> [a]
forall a. (Fractional a, Eq a) => CF a -> [a]
steed (a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
GCF a
b [(a, a)]
rest))
steed CF a
orig
= (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl a -> a -> a
forall a. Num a => a -> a -> a
(+) a
b0 [a]
dxs
where
(a
b0, (a
a1,a
b1):[(a, a)]
gcf) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
dxs :: [a]
dxs = a
a1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
b1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [(a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
1) a -> a -> a
forall a. Num a => a -> a -> a
* a
dx | (a
a,a
b) <- [(a, a)]
gcf | a
d <- [a]
ds | a
dx <- [a]
dxs]
ds :: [a]
ds = a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
b1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a -> a
forall a. Fractional a => a -> a
recip (a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
d) | (a
a,a
b) <- [(a, a)]
gcf | a
d <- [a]
ds]
lentz :: (Fractional a, Eq a) => CF a -> [a]
lentz :: CF a -> [a]
lentz = (a -> a) -> (a -> a -> a) -> (a -> a) -> CF a -> [a]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
lentzWith a -> a
forall a. a -> a
id a -> a -> a
forall a. Num a => a -> a -> a
(*) a -> a
forall a. Fractional a => a -> a
recip
{-# INLINE lentzWith #-}
lentzWith :: (Fractional a, Eq a) => (a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
lentzWith :: (a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
lentzWith a -> b
f b -> b -> b
op b -> b
inv (CF a
0 ( a
a :[a]
rest)) = (b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map b -> b
inv ((a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
lentzWith a -> b
f b -> b -> b
op b -> b
inv (a -> [a] -> CF a
forall a. a -> [a] -> CF a
CF a
a [a]
rest))
lentzWith a -> b
f b -> b -> b
op b -> b
inv (GCF a
0 ((a
a,a
b):[(a, a)]
rest)) = (b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (b -> b -> b
op (a -> b
f a
a) (b -> b) -> (b -> b) -> b -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. b -> b
inv) ((a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
lentzWith a -> b
f b -> b -> b
op b -> b
inv (a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
GCF a
b [(a, a)]
rest))
lentzWith a -> b
f b -> b -> b
op b -> b
inv CF a
c = (b -> a -> b) -> b -> [a] -> [b]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl b -> a -> b
opF (a -> b
f a
b0) ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
cs [a]
ds)
where
opF :: b -> a -> b
opF b
x a
y = b -> b -> b
op b
x (a -> b
f a
y)
(a
b0, [a]
cs, [a]
ds) = CF a -> (a, [a], [a])
forall a. (Fractional a, Eq a) => CF a -> (a, [a], [a])
lentzRecurrence CF a
c
lentzRecurrence :: (Fractional a, Eq a) => CF a -> (a,[a],[a])
lentzRecurrence :: CF a -> (a, [a], [a])
lentzRecurrence CF a
orig
| [(a, a)] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(a, a)]
terms = (a
b0,[],[])
| Bool
otherwise = (a
b0, [a]
cs, [a]
ds)
where
(a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
cs :: [a]
cs = [ a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
aa -> a -> a
forall a. Fractional a => a -> a -> a
/a
c | (a
a,a
b) <- [(a, a)]
terms | a
c <- a
b0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
cs]
ds :: [a]
ds = [a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/(a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
aa -> a -> a
forall a. Num a => a -> a -> a
*a
d) | (a
a,a
b) <- [(a, a)]
terms | a
d <- a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
ds]
modifiedLentz :: (Fractional a, Eq a) => a -> CF a -> [[a]]
modifiedLentz :: a -> CF a -> [[a]]
modifiedLentz = (a -> a) -> (a -> a -> a) -> (a -> a) -> a -> CF a -> [[a]]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
modifiedLentzWith a -> a
forall a. a -> a
id a -> a -> a
forall a. Num a => a -> a -> a
(*) a -> a
forall a. Fractional a => a -> a
recip
{-# INLINE modifiedLentzWith #-}
modifiedLentzWith :: (Fractional a, Eq a) => (a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
modifiedLentzWith :: (a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
modifiedLentzWith a -> b
f b -> b -> b
op b -> b
inv a
z (CF a
0 ( a
a :[a]
rest)) = ([b] -> [b]) -> [[b]] -> [[b]]
forall a b. (a -> b) -> [a] -> [b]
map ((b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map b -> b
inv ) ((a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
modifiedLentzWith a -> b
f b -> b -> b
op b -> b
inv a
z (a -> [a] -> CF a
forall a. a -> [a] -> CF a
CF a
a [a]
rest))
modifiedLentzWith a -> b
f b -> b -> b
op b -> b
inv a
z (GCF a
0 ((a
a,a
b):[(a, a)]
rest)) = ([b] -> [b]) -> [[b]] -> [[b]]
forall a b. (a -> b) -> [a] -> [b]
map ((b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (b -> b -> b
op (a -> b
f a
a) (b -> b) -> (b -> b) -> b -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. b -> b
inv)) ((a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
modifiedLentzWith a -> b
f b -> b -> b
op b -> b
inv a
z (a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
GCF a
b [(a, a)]
rest))
modifiedLentzWith a -> b
f b -> b -> b
op b -> b
inv a
z CF a
orig = [(Bool, b)] -> [[b]]
forall b. [(Bool, b)] -> [[b]]
separate (((Bool, b) -> (Bool, a) -> (Bool, b))
-> (Bool, b) -> [(Bool, a)] -> [(Bool, b)]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (Bool, b) -> (Bool, a) -> (Bool, b)
opF (Bool
False, a -> b
f a
b0) [(Bool, a)]
cds)
where
(a
b0, [(Bool, a)]
cs, [(Bool, a)]
ds) = a -> CF a -> (a, [(Bool, a)], [(Bool, a)])
forall a.
(Fractional a, Eq a) =>
a -> CF a -> (a, [(Bool, a)], [(Bool, a)])
modifiedLentzRecurrence a
z CF a
orig
cds :: [(Bool, a)]
cds = ((Bool, a) -> (Bool, a) -> (Bool, a))
-> [(Bool, a)] -> [(Bool, a)] -> [(Bool, a)]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Bool, a) -> (Bool, a) -> (Bool, a)
forall b. Num b => (Bool, b) -> (Bool, b) -> (Bool, b)
mult [(Bool, a)]
cs [(Bool, a)]
ds
mult :: (Bool, b) -> (Bool, b) -> (Bool, b)
mult (Bool
xa,b
xb) (Bool
ya,b
yb) = (Bool
xa Bool -> Bool -> Bool
|| Bool
ya, b
xb b -> b -> b
forall a. Num a => a -> a -> a
* b
yb)
opF :: (Bool, b) -> (Bool, a) -> (Bool, b)
opF (Bool
xa,b
xb) (Bool
ya,a
yb) = (Bool
xa Bool -> Bool -> Bool
|| Bool
ya, b -> b -> b
op b
xb (a -> b
f a
yb))
separate :: [(Bool, b)] -> [[b]]
separate [] = []
separate ((Bool
_,b
x):[(Bool, b)]
xs) = case ((Bool, b) -> Bool) -> [(Bool, b)] -> ([(Bool, b)], [(Bool, b)])
forall a. (a -> Bool) -> [a] -> ([a], [a])
break (Bool, b) -> Bool
forall a b. (a, b) -> a
fst [(Bool, b)]
xs of
([(Bool, b)]
xs, [(Bool, b)]
ys) -> (b
xb -> [b] -> [b]
forall a. a -> [a] -> [a]
:((Bool, b) -> b) -> [(Bool, b)] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (Bool, b) -> b
forall a b. (a, b) -> b
snd [(Bool, b)]
xs) [b] -> [[b]] -> [[b]]
forall a. a -> [a] -> [a]
: [(Bool, b)] -> [[b]]
separate [(Bool, b)]
ys
modifiedLentzRecurrence :: (Fractional a, Eq a) => a -> CF a -> (a,[(Bool, a)],[(Bool, a)])
modifiedLentzRecurrence :: a -> CF a -> (a, [(Bool, a)], [(Bool, a)])
modifiedLentzRecurrence a
z CF a
orig
| [(a, a)] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(a, a)]
terms = (a
b0, [], [])
| Bool
otherwise = (a
b0, [(Bool, a)]
cs, [(Bool, a)]
ds)
where
(a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
cs :: [(Bool, a)]
cs = [a -> (a -> a) -> (Bool, a)
forall b. a -> (a -> b) -> (Bool, b)
reset (a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
aa -> a -> a
forall a. Fractional a => a -> a -> a
/a
c) a -> a
forall a. a -> a
id | (a
a,a
b) <- [(a, a)]
terms | a
c <- a
b0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: ((Bool, a) -> a) -> [(Bool, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (Bool, a) -> a
forall a b. (a, b) -> b
snd [(Bool, a)]
cs]
ds :: [(Bool, a)]
ds = [a -> (a -> a) -> (Bool, a)
forall b. a -> (a -> b) -> (Bool, b)
reset (a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
aa -> a -> a
forall a. Num a => a -> a -> a
*a
d) a -> a
forall a. Fractional a => a -> a
recip | (a
a,a
b) <- [(a, a)]
terms | a
d <- a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: ((Bool, a) -> a) -> [(Bool, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (Bool, a) -> a
forall a b. (a, b) -> b
snd [(Bool, a)]
ds]
reset :: a -> (a -> b) -> (Bool, b)
reset a
x a -> b
f
| a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = (Bool
True, a -> b
f a
z)
| Bool
otherwise = (Bool
False, a -> b
f a
x)
sumPartialProducts :: Num a => [a] -> CF a
sumPartialProducts :: [a] -> CF a
sumPartialProducts [] = a -> [a] -> CF a
forall a. a -> [a] -> CF a
cf a
0 []
sumPartialProducts (a
x:[a]
xs) = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf a
0 ((a
x,a
1)(a, a) -> [(a, a)] -> [(a, a)]
forall a. a -> [a] -> [a]
:[(a -> a
forall a. Num a => a -> a
negate a
x, a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
x) | a
x <- [a]
xs])