module HOPS.GF.Series
( module HOPS.GF.Rat
, Series (..)
, polynomial
, series
, xpow
, nil
, infty
, precision
, coeffVector
, coeffList
, constant
, leadingCoeff
, rationalPrefix
, integerPrefix
, intPrefix
, eval
, (.*)
, (./)
, (!^!)
, (^!)
, (?)
, o
, blackDiamond
, derivative
, integral
, revert
, sec
, fac
, rseq
, rseq'
) where
import GHC.TypeLits
import Data.Proxy
import Data.Ratio
import Data.Maybe
import Data.Vector (Vector, (!), (//))
import qualified Data.Vector as V
import qualified Data.ByteString.Char8 as B
import HOPS.GF.Rat
import HOPS.Pretty
newtype Series (n :: Nat) = Series (Vector Rat)
instance Show (Series n) where
showsPrec p f =
let s = "series (Proxy :: Proxy " ++ show (precision f) ++ ") "
in showParen (p > 10) $ showString s . shows (coeffList f)
instance Eq (Series n) where
f == g = coeffVector f == coeffVector g
instance KnownNat n => Num (Series n) where
(Series u) + (Series v) = Series $ V.zipWith (+) u v
(Series u) (Series v) = Series $ V.zipWith () u v
(Series u) * (Series v) = Series $ u `mul` v
fromInteger x = polynomial (Proxy :: Proxy n) [fromIntegral x]
abs f = signum f * f
signum f = polynomial (Proxy :: Proxy n) [signum (leadingCoeff f)]
instance KnownNat n => Fractional (Series n) where
fromRational r = polynomial (Proxy :: Proxy n) [fromRational r]
(Series u) / (Series v) = Series $ u `divide` v
instance KnownNat n => Floating (Series n) where
pi = polynomial (Proxy :: Proxy n) [pi]
exp = exp'
log f = c0 log f + integral (derivative f / restrict f)
f ** g = if isConstant g then f ^! constant g else exp (g * log f)
sin f = c0 sin f + integral (derivative f * cos (restrict f))
cos f = c0 cos f integral (derivative f * sin (restrict f))
asin f = c0 asin f + integral (derivative f / sqrt (1 restrict f^(2::Int)))
acos f = c0 acos f integral (derivative f / sqrt (1 restrict f^(2::Int)))
atan f = c0 atan f + integral (derivative f / (1 + restrict f^(2::Int)))
sinh f = c0 sinh f + integral (derivative f * cosh (restrict f))
cosh f = c0 cosh f + integral (derivative f * sinh (restrict f))
asinh f = c0 asinh f + integral (derivative f / sqrt (restrict f^(2::Int) + 1))
acosh f = c0 acosh f + integral (derivative f / sqrt (restrict f^(2::Int) 1))
atanh f = c0 atanh f + integral (derivative f / (1 restrict f^(2::Int)))
sqrt f = f ^! (1/2)
instance Pretty (Series n) where
pretty f = B.concat ["{", B.intercalate "," (map pretty (coeffList f)), "}"]
coeffVector :: Series n -> Vector Rat
coeffVector (Series v) = v
coeffList :: Series n -> [Rat]
coeffList = V.toList . coeffVector
leadingCoeff :: Series n -> Rat
leadingCoeff f =
case dropWhile (==0) (coeffList f) of
[] -> 0
(c:_) -> c
rationalPrefix :: Series n -> [Rational]
rationalPrefix = mapMaybe maybeRational . takeWhile isRational . coeffList
integerPrefix :: Series n -> [Integer]
integerPrefix = mapMaybe maybeInteger . takeWhile isInteger . coeffList
intPrefix :: Series n -> [Int]
intPrefix = mapMaybe maybeInt . takeWhile isInt . coeffList
precision :: Series n -> Int
precision (Series v) = V.length v
polynomial :: KnownNat n => Proxy n -> [Rat] -> Series n
polynomial n xs = Series $ V.fromListN (fromInteger (natVal n)) (xs ++ repeat 0)
series :: KnownNat n => Proxy n -> [Rat] -> Series n
series n xs = Series $ V.fromListN (fromInteger (natVal n)) (xs ++ repeat Indet)
nil :: KnownNat n => Series n
nil = series (Proxy :: Proxy n) []
infty :: KnownNat n => Series n
infty = series (Proxy :: Proxy n) [DZ]
xpow :: KnownNat n => Int -> Series n
xpow k = polynomial (Proxy :: Proxy n) $ replicate k 0 ++ [1]
(.*) :: Series n -> Series n -> Series n
(.*) (Series u) (Series v) = Series $ V.zipWith (*) u v
(./) :: Series n -> Series n -> Series n
(./) (Series u) (Series v) = Series $ V.zipWith (/) u v
mul :: Vector Rat -> Vector Rat -> Vector Rat
mul u v =
let n = V.length u
in if n < 17
then convolution n u v
else let mu = V.mapM maybeInt u
mv = V.mapM maybeInt v
in case (mu, mv) of
(Just u', Just v') ->
V.fromList $ map fromIntegral (pmult n (V.toList u') (V.toList v'))
_ -> convolution n u v
convolution :: Int -> Vector Rat -> Vector Rat -> Vector Rat
convolution n u v = V.generate n $ \j -> sum [u!i * v!(ji) | i <- [0..j]]
peval :: [Int] -> Integer -> Integer
peval p x = foldr (\c s -> x * s + fromIntegral c) 0 p
maxNorm :: [Int] -> Integer
maxNorm = maximum . (0:) . map (abs . toInteger)
bound :: Int -> [Int] -> [Int] -> Integer
bound prec p q =
let b = max 1 (toInteger prec * maxNorm p * maxNorm q)
in 2^((2 :: Integer) + ceiling (logBase 2 (fromIntegral b) :: Double))
pmult :: Int -> [Int] -> [Int] -> [Integer]
pmult prec p q =
take prec $ unpack (peval p a * peval q a) ++ repeat 0
where
a = bound prec p q
unpack i = if i < 0 then map (*(1)) $ unp (i) else unp i
unp 0 = []
unp i = let (c,r) = quotRem i a
in if 2*r > a
then (r a) : unp (c+1)
else r : unp c
divide :: Vector Rat -> Vector Rat -> Vector Rat
divide u v | V.null v || V.null u = V.empty
divide u v =
let u' = V.tail u
v' = V.tail v
in case (V.head u, V.head v) of
(0, 0) -> V.snoc (divide u' v') Indet
(c, d) -> let q = c/d
qv' = V.map (q*) v'
in V.cons q $ divide (V.zipWith () u' qv') (V.init v)
derivative :: Series n -> Series n
derivative (Series v) = Series $ V.snoc u Indet
where
u = V.imap (\i a -> a * fromIntegral (i+1)) (V.tail v)
integral :: Series n -> Series n
integral (Series v) = Series (V.cons 0 u)
where
u = V.imap (\i a -> a / fromIntegral (i+1)) (V.init v)
revert :: Series n -> Series n
revert (Series u) = Series rev
where
rev | V.null u = V.empty
| V.head u /= 0 = V.fromListN n (repeat Indet)
| otherwise = iterate f ind !! n
where
n = V.length u
ind = V.fromListN n (repeat Indet)
one = V.fromListN (n1) (1 : repeat 0)
u' = V.tail u
f w = V.cons 0 $ one `divide` (u' `comp` V.init w)
eval :: Series n -> Rat -> Rat
eval (Series v) x = V.foldr (\c s -> x * s + c) 0 v
comp :: Vector Rat -> Vector Rat -> Vector Rat
comp _ v | V.null v = V.empty
comp u v =
case V.head v of
0 -> V.cons c (v' `mul` (u' `comp` w))
_ -> V.imap addc (v `mul` V.snoc (u' `comp` w) 0)
where
addc 0 a = a + c
addc _ a = a
u' = V.tail u
v' = V.tail v
c = V.head u
w = V.init v
dropTrailing :: Eq a => a -> [a] -> [a]
dropTrailing x = reverse . dropWhile (== x) . reverse
rseq :: Series n -> Series n
rseq f = Series $ V.replicate (precision f) 0 // zip cs (repeat 1)
where
cs = dropTrailing 0 [ x | x <- intPrefix f, x >= 0, x < precision f ]
rseq' :: Series n -> Series n
rseq' f = Series $ V.map (1) $ coeffVector (rseq f)
infixr 7 ?
(?) :: KnownNat n => Series n -> Series n -> Series n
(?) (Series v) c | c == 0 = series (Proxy :: Proxy n) [v ! 0]
(?) (Series v) c =
series (Proxy :: Proxy n) $ map (v !) $
dropTrailing 0 [ x | x <- intPrefix c, x >= 0, x < precision c ]
infixr 7 `o`
o :: Series n -> Series n -> Series n
o (Series u) (Series v) = Series $ u `comp` v
c0 :: KnownNat n => (Rat -> Rat) -> Series n -> Series n
c0 f (Series v) = fromRat $ f (V.head v)
constant :: Series n -> Rat
constant (Series v) = V.head v
fromRat :: KnownNat n => Rat -> Series n
fromRat c = polynomial (Proxy :: Proxy n) [c]
kRestrict :: Int -> Series n -> Series n
kRestrict k (Series v) = Series (v // [(k, Indet)])
restrict :: Series n -> Series n
restrict f = kRestrict (precision f 1) f
exp' :: KnownNat n => Series n -> Series n
exp' f = expAux (precision f 1) f
where
expAux 0 _ = 1
expAux d g =
let h = expAux (d1) (kRestrict d g)
in c0 exp g + integral (derivative g * h)
(!^!) :: Rat -> Rat -> Rat
(!^!) c r = let Series v = polynomial (Proxy :: Proxy 2) [c] ^! r in V.head v
(^!) :: KnownNat n => Series n -> Rat -> Series n
(^!) _ DZ = polynomial (Proxy :: Proxy n) [DZ]
(^!) _ Indet = polynomial (Proxy :: Proxy n) [Indet]
(^!) f s@(Val r) =
case (numerator r, denominator r) of
(0, _) -> polynomial (Proxy :: Proxy n) [1]
(n, _) | n < 0 -> 1 / (f ^! (s))
(n, 1) -> f ^^ n
(n, k) ->
case (d `mod` fromInteger k, k) of
(0, 2) -> y * squareRoot (f/xpow d) ^^ n
(0, _) -> y * exp' (fromRational r * log (f/xpow d))
(_, _) -> polynomial (Proxy :: Proxy n) [Indet]
where
d = fromInteger $ leadingExponent f :: Int
y = xpow ((d `div` fromInteger k) * fromInteger n)
leadingExponent :: Series n -> Integer
leadingExponent f =
case span (==0) (rationalPrefix f) of
(_ ,[]) -> 0
(xs,_ ) -> fromIntegral (length xs)
isConstant :: Series n -> Bool
isConstant (Series u) = V.all (==0) (V.tail u)
squareRoot :: KnownNat n => Series n -> Series n
squareRoot f = squareRoot' (precision f 1) f
where
squareRoot' 0 _ = 1
squareRoot' d g =
let h = 2 * squareRoot' (d1) (kRestrict d g)
in c0 sqrt g + integral (derivative g / h)
blackDiamond :: Series n -> Series n -> Series n
blackDiamond (Series u) (Series v) =
Series $ V.generate (V.length u) $ \n ->
sum [ u!(a+c) * v!(b+c) * multinomial [a,b,c]
| [a,b,c] <- compositions 3 n
]
compositions :: Int -> Int -> [[Int]]
compositions 0 0 = [[]]
compositions 0 _ = []
compositions 1 n = [[n]]
compositions 2 n = [[ni,i] | i <- [0..n]]
compositions k 0 = [ replicate k 0 ]
compositions k n = [0..n] >>= \i -> map (i:) (compositions (k1) (ni))
sec :: KnownNat n => Series n -> Series n
sec f = 1 / cos f
fac :: KnownNat n => Series n -> Series n
fac f | isConstant f = polynomial (Proxy :: Proxy n) [factorial (constant f)]
| otherwise = polynomial (Proxy :: Proxy n) [Indet]