-- |
-- Module : LinRec
-- Copyright : (c) 2006-2007 Bertram Felgenhauer
-- License : BSD3
--
-- Maintainer : Bertram Felgenhauer
-- Stability : experimental
-- Portability : portable
--
-- Module for evaluating linear recurrences in arbitrary rings.
--
-- The linear recurrence a_{k+n} = b_0a_k + b_1a_{k+1} ... + b_{n-1}a_{k+n-1}
-- is represented by the list [b_0,...,b_{n-1}].
--
module LinRec (
generate,
step,
get,
) where
import Data.List (tails)
-- | return scalar product of two vectors (given as lists)
(.*.) :: Num a => [a] -> [a] -> a
a .*. b = sum $ zipWith (*) a b
-- | apply the recurrence exactly once to a prefix, yielding the next prefix
--
-- @
-- step [1,1] [2,3]
-- @
--
-- results in [3,5]
step :: Num a => [a] -> [a] -> [a]
step rec pre = tail $ pre ++ [rec .*. pre]
-- | generate the sequence with prefix pre.
--
-- @
-- generate [1,1] [0,1]
-- @
--
-- results in the Fibonacci sequence [0,1,1,2,...]
generate :: Num a => [a] -> [a] -> [a]
generate rec pre = let res = pre ++ map (rec .*.) (tails res) in res
-- | Get the k-th element of the sequence with prefix pre.
-- @get@ obeys the law
--
-- @
-- get rec pre k == generate rec pre !! k
-- @
--
-- but takes O(n^2 * log(k)) @Num@ operations, where n is the order of
-- the recurrence relation, while @generate@ takes O(n*k) operations.
--
-- @
-- get [1,1] [0,1] 200000
-- @
--
-- calculates the 200,000th Fibonacci number reasonably quickly.
get :: Num a => [a] -> [a] -> Int -> a
get rec pre k = u0inv [] (u k) .*. pre where
-- Let
-- / 0 1 ... 0 \
-- T = | . . . |
-- | 0 0 ... 1 |
-- \ b_0 b_1 ... b_{n-1} /
--
-- The code below is largely concerned with calculating T^k efficiently.
--
-- Let
-- u_0 = ... = u_{k-2} = 0, u_{k-1} = 1,
-- u_n+k =
--
-- and define
-- / u_k u_{k+1} ... u_{k+n-1} \
-- U_k = | . . . |
-- \ u_{k+n-1} u_{k+n} ... u_{k+2n-2} /
--
-- Then, U_k = T^k U_0, and T^k = U_k U_0^{-1}.
--
-- Each U_k can be extrapolated easily from the first row by applying
-- the linear recurrence, and U_0 was picked such that it's unitarian,
-- and thus can be inverted without using divisions.
n = length rec
-- u_0 is the first row of U_0, u_n the first row of U_n.
u_0 = iterate (0:) [1] !! (n-1)
u_n = take n $ drop n $ generate rec u_0
-- u k calculates the first row of U_k,
-- using U_{2k} = U_k U_0^{-1} U_k
u 0 = u_0
u k | odd k = step rec $ u (k-1)
| otherwise = let u' = u (k `div` 2)
u'' = u0inv [] u'
extrapolate xs = take n $ tails $ generate rec xs
in map (u'' .*.) $ extrapolate u'
-- u0inv [] x calculates the first row of X U_0^{-1}, where X is the
-- matrix obtained by extrapolating the first row, x, by using the linear
-- recurrence. It's a streamlined version of Gaussian elimination on U_0.
u0inv acc [] = acc
u0inv acc (x:xs) = u0inv (x:acc) $ zipWith (-) xs $ map (x*) u_n