{-# LANGUAGE NoImplicitPrelude
#-}
-----------------------------------------------------------------------------
-- |
-- Module : Math.Combinatorics.Species.CycleIndex
-- Copyright : (c) Brent Yorgey 2010
-- License : BSD-style (see LICENSE)
-- Maintainer : byorgey@cis.upenn.edu
-- Stability : experimental
--
-- The Newton-Raphson iterative method for computing with recursive
-- species. Any species @T@ which can be written in the form @T =
-- X*R(T)@ (the species of "@R@-enriched rooted trees") may be
-- computed by a quadratically converging iterative process. In fact
-- we may also compute species of the form @T = N + X*R(T)@ for any
-- integer species @N@, by iteratively computing @T' = X*R(T' + N)@
-- and then adding @N@.
--
-----------------------------------------------------------------------------
module Math.Combinatorics.Species.NewtonRaphson
(
newtonRaphsonIter
, newtonRaphson
, newtonRaphsonRec
, solveForR
) where
import NumericPrelude
import PreludeBase
import Math.Combinatorics.Species.Class
import Math.Combinatorics.Species.AST
import Math.Combinatorics.Species.AST.Instances (reflect)
import Math.Combinatorics.Species.Simplify
import Data.Typeable
import Control.Monad (guard)
import Data.List (delete)
-- | A single iteration of the Newton-Raphson method.
-- @newtonRaphsonIter r k a@ assumes that @a@ is a species having
-- contact of order @k@ with species @t = x '*' (r ``o`` t)@ (that
-- is, @a@ and @t@ agree on all label sets of size up to and
-- including @k@), and returns a new species with contact of order
-- @2k+2@ with @t@.
--
-- See BLL section 3.3.
newtonRaphsonIter :: Species s => s -> Integer -> s -> s
newtonRaphsonIter r k a = a + sum as
where p = x * (r `o` a)
q = x * (oneHole r `o` a)
ps = map (p `ofSizeExactly`) [k+1..2*k+2]
qs = map (q `ofSizeExactly`) [1..k+1]
as = zipWith (+) ps
(map (sum . zipWith (*) qs) $ map reverse (inits' as))
-- | Lazier version of inits.
inits' :: [a] -> [[a]]
inits' xs = [] : inits'' xs
where inits'' [] = []
inits'' (x:xs) = map (x:) (inits' xs)
-- | Given a species @r@ and a desired accuracy @k@, @'newtonRaphson'
-- r k@ computes a species which has contact at least @k@ with the
-- species @t = x '*' (r ``o`` t)@.
newtonRaphson :: Species s => s -> Integer -> s
newtonRaphson r n = newtonRaphson' 0 0
where newtonRaphson' a k
| k >= n = a
| otherwise = newtonRaphson' (newtonRaphsonIter r k a) (2*k + 2)
-- | @'newtonRaphsonRec' f k@ tries to compute the recursive species
-- represented by the code @f@ up to order at least @k@, using
-- Newton-Raphson iteration. Returns 'Nothing' if @f@ cannot be
-- written in the form @f = X*R(f)@ for some species @R@.
newtonRaphsonRec :: (ASTFunctor f, Species s) => f -> Integer -> Maybe s
newtonRaphsonRec code k = fmap (\(n,r) -> n + newtonRaphson r k) (solveForR code)
-- | Given a code @f@ representing a recursive species, try to find an
-- integer species N and species R such that @f = N + X*R(f)@. If
-- such species can be found, return @'Just' (N,R)@; otherwise
-- return 'Nothing'.
solveForR :: (ASTFunctor f, Species s) => f -> Maybe (s, s)
solveForR code = do
let terms = sumOfProducts . erase' $ apply code (TRec code)
guard . not . null $ terms
-- If there is a constant term, it will be the first one; pull it
-- out.
let (n, terms') = case terms of
([One] : ts) -> (One, ts)
([N n] : ts) -> (N n, ts)
ts -> (Zero, ts)
-- Now we need to be able to factor an X out of the rest.
guard $ all (X `elem`) terms'
-- XXX this is wrong, what if there are still occurrences of X remaining?
-- Now replace every recursive occurrence by (n + X).
let r = foldr1 (+) $ map ( foldr1 (*)
. map (substRec code (n + x))
. delete X)
terms'
return (reflect n, reflect r)