------------------------------------------------------------------------------ |-- Module : Math.LinearRecursive.Monad-- Copyright : (c) Bin Jin 2011-- License : BSD3---- Maintainer : bjin1990+haskell@gmail.com-- Stability : experimental-- Portability : portable---- A monad to calculate linear recursive sequence efficiently. Matrix-- multiplication and fast exponentiation algorithm are used to speed-- up calculating the number with particular index in the sequence. This-- library also provides a monadic DSL to describe the sequence.---- As an example, here is the fibonacci sequence-- -- >fib = do-- > [f0, f1] <- newVariables [1, 1]-- > f0 <:- f0 <+> f1-- > return f1-- -- >>> map (runLinearRecursive fib) [0..10]-- [1,1,2,3,5,8,13,21,34,55,89]-- ----------------------------------------------------------------------------moduleMath.LinearRecursive.Monad(-- * Vector-- ** vector typesVectorLike(..),LinearCombination,Variable-- ** vector operators and constant,(<+>),(<->),(<*),(*>),zeroVector-- * Polynomial,Polynomial,P.x-- * Monad,LinearRecursive,newVariable,newVariables,(<:-),(<+-),runLinearRecursive,simulateLinearRecursive-- * Utility,getConstant,getPartialSum,getStep,getPowerOf,getPolynomial,getPartialSumWith)whereimportControl.Monad(zipWithM_)importControl.Applicative((<$>))importData.IntMap(IntMap)importqualifiedData.IntMapasIntMapimportMath.LinearRecursive.Internal.VectorimportMath.LinearRecursive.Internal.MatriximportMath.LinearRecursive.Internal.Polynomialhiding(fromList,toList,x)importqualifiedMath.LinearRecursive.Internal.PolynomialasP-- | A vector represents linear combination of several variables.typeLinearCombination=Vector-- | An unit vector represents dependence on a particular variable.typeVariable=Vector1dataLRVariablea=LRV{initialValue::a,dependency::LinearCombinationa}dmap::Numa=>(LinearCombinationa->LinearCombinationa)->LRVariablea->LRVariableadmapf(LRVvaldep)=LRVval(fdep)typeLRVariablesa=IntMap(LRVariablea)-- | The monad to specify the calculation of next number of a linear recursive sequence.---- All linear recursive sequences can be generated by iteration, the next number can-- be represented by linear combination of some previous numbers. This can be regarded-- as linear transformation between states, and it's actually multiply a transform matrix.---- In order to formalize and simply this procedure, this monad use mutable-like variables to -- denote the states, and mutable-like assignment to denote the transform matrix.---- To evaluate this sequence, the monad will be simulated step by step, after each step, all-- variable will be updated. Besides, if the monad returns a 'LinearCombination', a number-- will be generated each step. (well, actual calculation uses fast exponentiation algorithm -- to speed up this calculation)--dataLinearRecursiveab=LR{unLR::Int->(b,Int,LRVariablesa->LRVariablesa)}-- unLR prevDeclaredVars = (return value, newDeclaredVars, changes to variables)instanceNuma=>Functor(LinearRecursivea)wherefmapfm=m>>=return.finstanceNuma=>Monad(LinearRecursivea)wherereturna=LR(const(a,0,id))a>>=b=LR$\v->let(ra,nva,ma)=unLRav(rb,nvb,mb)=unLR(bra)(v+nva)in(rb,nva+nvb,mb.ma)-- | Declare a new variable, with its initial value (the value before step 0).---- >test = do-- > v <- newVariable 1-- > v <:- v <+> v-- > return v---- >>> map (runLinearRecursive test) [0..10]-- [1,2,4,8,16,32,64,128,256,512,1024]newVariable::(Eqa,Numa)=>a->LinearRecursivea(Variablea)newVariableval0=LR$\v->(vector1v,1,IntMap.insertvvariable)wherevariable=LRV{initialValue=val0,dependency=zeroVector}-- | Declare a new sequence, with their initial value.---- After each step, each variable except the first one will be assigned with-- the value of its predecessor variable before this turn. ---- It's not encouraged to assign any value to the variables other -- than the first one.---- >test = do-- > [v1, v2, v3] <- newVariables [1,2,3]-- > v1 <:- v3-- > return v3---- >>> map (runLinearRecursive test) [0..10]-- [3,2,1,3,2,1,3,2,1,3,2]newVariables::(Eqa,Numa)=>[a]->LinearRecursivea[Variablea]newVariablesvals=doret<-mapMnewVariablevalszipWithM_(<:-)(tailret)retreturnret-- | return a constent number. Use one extra variable.---- >>> map (runLinearRecursive (getConstant 3)) [0..10]-- [3,3,3,3,3,3,3,3,3,3,3]getConstant::(Eqa,Numa)=>a->LinearRecursivea(LinearCombinationa)getConstantval=doone<-newVariable1one<:-onereturn(toVectorone*>val)-- | return sum of a linear combination in steps before current one. Use one extra variable.---- >>> map (runLinearRecursive (getConstant 3 >>= getPartialSum)) [0..10]-- [0,3,6,9,12,15,18,21,24,27,30]getPartialSum::(Eqa,Numa)=>LinearCombinationa->LinearRecursivea(LinearCombinationa)getPartialSumval=dos<-newVariable0s<:-s<+>valreturn(toVectors)-- | return the current step number. Use two extra variables.---- >>> map (runLinearRecursive getStep) [0..10]-- [0,1,2,3,4,5,6,7,8,9,10]getStep::(Eqa,Numa)=>LinearRecursivea(LinearCombinationa)getStep=getConstant1>>=getPartialSum-- | @getPowerOf a@ return power of @a@ with order equal to current step number. -- Use one extra variable.---- >>> map (runLinearRecursive (getPowerOf 3)) [0..10]-- [1,3,9,27,81,243,729,2187,6561,19683,59049]getPowerOf::(Eqa,Numa)=>a->LinearRecursivea(LinearCombinationa)getPowerOfa=doprod<-newVariable1prod<:-prod*>areturn(toVectorprod)-- | given n polynomials, the i-th (0 indexed) polynomial's degree is i and with first -- coeff equal to one. find the linear combination for x^i for each i in [0, n)inverseTrans::(Eqa,Numa)=>[Polynomiala]->MatrixainverseTranspolys=inverseMatrixDiag1mawheren=lengthpolysma=matrix[[vcomponent(unPolypolyi)j|j<-[0..n-1]]|polyi<-polys]getPartialSumWith::(Eqa,Numa,VectorLikev)=>Polynomiala->va->LinearRecursivea(LinearCombinationa)getPartialSumWithpolyv|n<0=returnzeroVector|otherwise=dobasisValue<-go(toVectorv)0letvars=map(foldl(<+>)zeroVector.zipWith(*>)basisValue)transreturn$foldl(<+>)zeroVector[powi*>coeffi|(i,powi)<-zip[0..]vars,letcoeffi=vcomponentveci]wheren=degreepolybasisPoly=scanl(*)1[P.x+fromIntegrali|i<-[1..n]]goprevpos|pos>n=return[]|otherwise=donext<-(*>fromIntegral(pos`max`1)).(<+>prev)<$>getPartialSumprev(next:)<$>gonext(pos+1)trans=unMatrix(inverseTransbasisPoly)vec=unPolypoly-- | @getPolynomial poly@ evaluate polynomial @poly@ with variable @x@ replaced by current step number. -- Use @n@ extra variables, where @n@ is the degree of @poly@---- >>> map (runLinearRecursive (getPolynomial ((x+1)^2))) [0..10]-- [1,4,9,16,25,36,49,64,81,100,121]getPolynomial::(Eqa,Numa)=>Polynomiala->LinearRecursivea(LinearCombinationa)getPolynomialpoly=newVariable1>>=getPartialSumWithpoly-- | Variable accumulated assignment. @v \<+- a@ replace variable @v@ with @v \<+\> a@.---- Be aware that @v@ will be zero before any assignment.(<+-)::(Eqa,Numa,VectorLikev)=>Variablea->va->LinearRecursivea()(<+-)vardep=LR(const((),0,IntMap.adjust(dmap(<+>toVectordep))(unVector1var)))-- | Variable assignment. @v \<:- a@ replace variable @v@ with @a@ after this step. -- If there are multiple assignments to one variable, only the last one counts.(<:-)::(Eqa,Numa,VectorLikev)=>Variablea->va->LinearRecursivea()(<:-)vardep=LR(const((),0,IntMap.adjust(dmap(const(toVectordep)))(unVector1var)))infix1<:-,<+-buildMatrix::(Eqa,Numa)=>LRVariablesa->(Matrixa,Matrixa)buildMatrixmapping=(matrixtrans,matrix$map(:[])initValues)whereinitValues=mapinitialValue(IntMap.elemsmapping)rawDep=map(unVector'.dependency)(IntMap.elemsmapping)varCount=lengthinitValuestrans=map(\m->[IntMap.findWithDefault0im|i<-[0..varCount-1]])rawDep-- | /O(v^3 * log n)/, where /v/ is the number of variables, and /n/ is steps to simulate.---- @runLinearRecursive m n@ simulate the monad by @n@ steps, and return the actual value denoted-- by returned 'LinearCombination'.---- n must be non-negative.runLinearRecursive::(Eqa,Numa,Integralb,VectorLikev)=>LinearRecursivea(va)->b->arunLinearRecursive_steps|steps<0=error"runLinearRecursive: steps must be non-negative"runLinearRecursivemsteps=sum[head(res!!i)*ai|(i,ai)<-IntMap.assocs(unVector'target)]where(target,_,g)=unLRm0dep=gIntMap.empty(trans,initCol)=buildMatrixdepres=unMatrix'(trans^steps*initCol)-- | /O(v^2 * n)/. similar to @runLinearRecursive@, but return an infinite list instead of a particular index.simulateLinearRecursive::(Eqa,Numa,VectorLikev)=>LinearRecursivea(va)->[a]simulateLinearRecursivem=map(\res->sum[head(res!!i)*ai|(i,ai)<-IntMap.assocs(unVector'target)])colswhere(target,_,g)=unLRm0dep=gIntMap.empty(trans,initCol)=buildMatrixdepcols=mapunMatrix'$scanl(flip(*))initCol(repeattrans)