# Maple procedures:
# The procedures calculate generating functions for linear recurrences
# given the coefficients for some polynomial terms of low order plus
# the coefficients of the linear recurrence of the higher order terms.
# Richard J. Mathar, 2007-11-14, http://www.strw.leidenuniv.nl/~mathar
# Decompose a polynomial in n into a list of coefficients
# [N0,N1,N2,...] such that N0+N1*n+N2*n(n-1)+N3*n(n-1)(n-2)+... equals the polynomial.
# arguments:
# polyn: the polynomial with n as the variable
# returns:
# the list of the coefficients [N0,N1,N2,...]
#
# Example:
# If the argument polyn is n^2+7*n*(n-1)+3 on input, the function returns
# [3,1,8] because 3+1*n+8*n(n-1)= n^2+7*n*(n-1)+3.
#
# Note: this is only an auxiliary procedure called from other functions
# below.
FallFact := proc(polyn)
local i,o,c,t ;
c := [0] ;
# Work around what seems to be a bug in Maple 9.0 which returns
# a polynomial degree of -infnity if the polynomial equals 0.
if polyn = 0 then
;
else
c := [seq(0,i=0..degree(polyn,n))] ;
c := subsop(1=coeftayl(polyn,n=0,0),c) ;
for o from 1 to degree(polyn,n) do
for t from 0 to o do
op(t+1,c)+coeftayl(polyn,n=0,o)*combinat[stirling2](o,t) ;
c := subsop(t+1=%,c) ;
end do:
end do:
end if ;
RETURN(c) ;
end proc:
# From a list coefficients [[c0,d0],[c1,d1],...] construct a polynomial
# sum_i c_i*n(n-1)(n-2)....*(n-di+1)
# arguments:
# poch: the list with the duples [ci,di] of coefficients and depth of
# the falling factorial.
# returns:
# the polynomial
#
# Example:
# If the argument polyn is [[2,5],[7,2]], the polynomial
# 2n(n-1)(n-2)(n-3)(n-4)+7n(n-1) is returned.
#
FallFactInv := proc(poch)
local p,t,d ;
p := 0 ;
for t in poch do
p := p+op(1,t)*mul(n-d,d=0..op(2,t)-1) ;
end do:
RETURN(p) ;
end proc:
# Print the taylor coefficients of the function
# @param genF the function of x
# @param ord the order of the expansion
# @since 2008-04-16
GenFVerify := proc(genF,ord)
local j,g,T ;
print(genF) ;
if type(genF, ratpoly) = true then
g := convert(genF,parfrac,x) ;
print(g) ;
else
g := genF ;
end if ;
T := taylor(g,x=0,ord+1) ;
# The following step is introduced to work around a Maple 9.5
# bug which would otherwise create a coeftayl(T,x=0,3)=0 for
# T=x+2*x^3.
T := convert(T,polynom) ;
for j from 0 to ord do
printf("%a,", coeftayl(T,x=0,j) ) ;
end do ;
printf("\n") ;
end proc:
# Return a generating function with x as the independent variable
# which is derived from the generating function genF (again x as the
# main variable) with each Taylor coefficient a(n) multiplied by
# the polynomial polyn, a polynomial with independent variable n.
#
# Example:
# GenFpolyMul(n^3+n^2,1/(1-2*x)) returns 4x(1+4x)/(2x-1)^4, where
# the original expansion is 1/(1-2x)=1+2x+4x^2+8x^3+16x^4+...
# and the expansion returned is 4x+48x^2+288x^3+1280x^4+..,
# where a(3)=288=8*(3^3+3^2) has been created from the coefficient
# 8 of the input's Taylor expansion by multiplication with n^3+n^2
# at n=3.
#
# @since 2007-12-01
GenFpolyMul := proc(polyn,genF,verb)
local gf,stir,j ;
stir := FallFact( polyn ) ;
gf :=0 ;
# accumulate expansion coefficients of falling factorials
for j from 0 to nops(stir)-1 do
if j > 0 then
gf := gf+ op(j+1,stir)*x^j* diff(genF,x$j) ;
else
gf := gf+ op(j+1,stir)* genF ;
end if ;
end do ;
if verb then
GenFVerify(gf,18) ;
GenFpurrs(gf) ;
fi ;
RETURN(factor(gf)) ;
end proc:
# Compute the generating function of constants raised to the nth power.
# arguments:
# cbn: a list of pairs of the form [[c1,b1],[c2,b2],[c3,b3],...]
# which symbolize a sum of powers of the form c1*b1^n+c2*b2^n+c3*b3^n+...
# The procedure returns a generating function which matches this sum
# in the coefficient of the term for x^n, n >= offs. The ci may be
# polynomials in n, the bi must be independent of n.
# offs: an integer offset value. This is the minimum n (minimum power in
# the x^n Taylor expansion of the generating function) at which the sum
# should be matched. If set to a value larger than 0, the constant, linear
# etc powers will get a zero coefficient from the generating function.
# verb: a boolean value. If set to true, the first few Taylor coefficients
# of the generating function are displayed as a check.
#
# Example 1:
# GenPown([[1/3,-7],[6,2]],3,false) returns -x^3(1694x-199)/3/(1+7x)/(2x-1)
# = -199/3*x^3+2689/3*x^4+... where (1/3)*(-7)^3+6*2^3= -199/3 is
# the Taylor coefficient of the generating function at n=3, where
# (1/3)*(-7)^4+6*2^4= 2689/3 is the Taylor coefficient of the
# generating function at n=4, and the coefficients of x^0, x^1 and x^2
# are zero as demanded by the offset of 3.
#
# Example 2:
# GenPown([[1/3,-7],[6,2]],0,true) returns -(124x+19)/3/(1+7x)/(2x-1)
# = 19/3+29/3*x... where (1/3)*(-7)^0+6*2^0= 19/3 is
# the Taylor coefficient of the generating function at n=0, where
# (1/3)*(-7)^1+6*2^1= 29/3 is the Taylor coefficient of the
# generating function at n=1.
#
# Note: This is only an auxiliary routine for GenPolyn() which has a more
# generic interface and returns generating functions of a more general
# input specification.
GenPown := proc(cbn,offs,verb)
local gf,i,j,geos,stir;
gf :=0 ;
for i in cbn do
geos := 1/(1-op(2,i)*x) - add((op(2,i)*x)^j,j=0..offs-1) ;
gf := gf + GenFpolyMul( op(1,i), geos, false ) ;
end do:
# If the verb parameter has been set to true, display the first
# coefficients of the Taylor series as a check.
if verb then
interface(prettyprint=0) ;
print(taylor(gf,x=0,23)) ;
end if ;
RETURN(factor(gf)) ;
end proc:
# Compute the generating function for the non-recurrent case
# of the n-th coefficient of the Taylor expansion being a polynomial in
# n plus a sum of powers with exponents n.
# arguments:
# polyn: a polynomial with n the main parameter
# offs: an integer offset
# verb: a parameter to display the first few terms of the series
# pown: a specification of a sum of n-th powers as described
# for the parameter cbn in the procedure GenPown
# returns:
# The generating function in the independent variable x which
# has coefficient zero for the components from x^0 up to x^(offs-1),
# and has coeffients equal to the sum polyn+pown for the powers
# x^n where n>= offs. Here, pown is encoded on input as described
# for the first parameter of the procedure GenPown.
#
# Example 1:
# GenPolyn(1,0,20,[]) returns 1/(1-x), where the first argument
# represents the polynomial 1 (the degenerate case where only a
# constant term is present), where the offset has been set to zero,
# and where the contributions of n-th powers is absent. The
# generating function 1/(1-x)=1+x+1*x^2+1*x^3+1*x^4+.. has each
# coefficient, starting with the power x^offs=x^0, set to the polynomial.
#
# Example 2:
# GenPolyn(n^3-1,0,20,[[2,5]]) returns
# (1+x-7x^2-11x^3-8x^4)/(1-x)^4/(1-5x), and displays the list of
# Taylor coeficients 1, 10, 57, 276 etc, where the first argument
# represents the polynomial n^3-1, where the offset has been set to zero,
# and where the contributions of n-th powers reads 2*5^n. The Taylor
# coefficients of the generating functions are to be interpreted as
# 1=0^3-1+2*5^0, 10=1^3-1+2*5^1, 57=2^3-1+2*5^2 etc, just the sum of
# the polynomial and the sum (here: single term) of the n-th powers.
#
# Example 3:
# GenPolyn(n*(n-1)*(n-2)*(n-3),0,20,[]) returns the generating function
# for the specification a(n)=n(n-1)(n-2)(n-3), 24x^4/(1-x)^5,
# see A052762 in the OEIS.
#
# Example 4:
# GenPolyn(0,0,20,[[n^2,3],[3-n,5]]) computes the generating
# function for coefficients n^2*3^n+(3-n)*5^n, that is
# -(-44*x+240*x^2-636*x^3+765*x^4+3)/(-1+3*x)^3/(-1+5*x)^2
# = 3+13x+61x^2+243x^3+....
#
# Note: a more general interface is provided by GenFLinRec(), which
# also handles recurrences.
GenPolyn := proc(polyn,offs,verb,pown)
local fn,gf,i ;
gf := GenPown([[polyn,1]],offs,false) ;
gf := gf+ GenPown(pown,offs,false) ;
# To beautify the output, we try to factorize the result.
gf := factor(gf) ;
gf := sort(gf,"<") ;
# If verbosity is demanded, the first few Taylor coefficients
# of the generating function are displayed as a check.
if verb >0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Compute generating function for linear recurrences with constant coefficients
# and optional inhomogeneous elements in the recurrences of some polynomial
# format.
# arguments:
# initL: a list of the first initial terms of the Taylor expansion of
# the generating function, [a0,a1,....], where the first entry specifies
# the coefficient of x^off, the second the coefficient of x^(off+1)
# and so on. The list is at least as long as the list of elements in
# the argument recurL. If it is longer, the recurrence specified by
# recurL is "delayed" to be valid only at higher powers of x.
# recurL: a list of constant coefficients [d1,d2,d3,..] to specify the
# linear relation a(n)=d1*a(n-1)+d2*a(n-2)+d3*a(n-3)+...
# Ignoring some of the earlier coefficients is simply achieved by
# specifying a 0 at their positions. Note that the leftmost element
# in the list recurL is attached to Taylor coefficients just before
# the element a(n), the next one to the element even one index earlier
# etc.
# off: an integer offset; if larger than zero, the first off elements
# of the Taylor expansion of the generating function (starting at the
# constant term x^0) become zero.
# polyn: optional additional inhomogeneity at the right hand side of the
# recurrence, a polynomial in the variable n. To be set to zero if not used.
# a(n)=d1*a(n-1)+d2*a(n-2)+...+polyn+pown
# The format is the same as for the polyn argument of GenPolyn().
# verb: if set to a value larger than zero, the number of expansion coefficients
# printed.
# pown: optional additional inhomogeneity at the right hand side of the
# recurrence, a sum of n-th powers of constants.
# The format is the same as for the pown argument of GenPolyn().
#
# Example 1:
# GenFLinRec( [1,5], [0], 0, 0, 10,[[4,2]]) returns the generating function
# specified by a(0)=1, a(1) =5, a(n)=4*2^n, see OEIS A132479.
#
# Example 2:
# GenFLinRec( [1,7], [12,-1], 1, -5, 10,[]) returns the generating function
# specified by a(1)=1, a(2)=7, a(n)=12a(n-1)-a(n-2)-5, see OEIS A133272
#
# Example 3:
# GenFLinRec( [1,2,2], [0,0,1], 0, 0, 10,[]) returns the generating
# function specified by a(0)=1, a(1)=a(2)=2, a(n)=a(n-3) (so it is a periodic
# sequence of the a(n) in the taylor series a(n)*x^n with period 3).
#
# Example 4:
# GenFLinRec( [0,1], [3,k], 0, 0, 10,[]) returns the generating function
# for a(0)=0, a(1)=1, a(n)=3a(n-1)+k*a(n-2), see A083857 in the OEIS.
#
# Example 5:
# GenFLinRec( [1,10], [2,1], 0, 7*n+1, 10,[]) returns the generating
# function for the case a(0)=1, a(1)=10, a(n)=2a(n-1)+a(n-2)+7n+1,
# see A051959 in the OEIS.
#
GenFLinRec := proc(initL, recurL, off, polyn, verb, pown)
local n, # numerator
d, # denominator
i,j,g ,
polynOff;
n := 0 ;
for i from 1 to nops(initL) do
n := n+op(i,initL)*x^(i-1+off) ;
end do:
for i from 1 to nops(recurL) do
for j from 1 to nops(initL)-i do
n := n-op(i,recurL)*op(j,initL) * x^(i+j-1+off) ;
end do:
end do:
n := n+ GenPolyn(polyn,nops(initL)+off,0,pown) ;
d := 1 ;
for i from 1 to nops(recurL) do
d := d-op(i,recurL)*x^i ;
end do:
g := factor(n/d) ;
if verb > 0 then
GenFVerify(convert(g,parfrac,x),verb) ;
try:
GenFpurrs(g) ;
catch:
end try:
fi ;
RETURN(g) ;
end proc:
# Compute generating function for linear recurrences with constant coefficients
# and optional inhomogeneous elements in the recurrences of some gen. func.
# format.
# arguments:
# initL: a list of the first initial terms of the Taylor expansion of
# the generating function, [a0,a1,....], where the first entry specifies
# the coefficient of x^off, the second the coefficient of x^(off+1)
# and so on. The list is at least as long as the list of elements in
# the argument recurL. If it is longer, the recurrence specified by
# recurL is "delayed" to be valid only at higher powers of x.
# recurL: a list of constant coefficients [d1,d2,d3,..] to specify the
# linear relation a(n)=d1*a(n-1)+d2*a(n-2)+d3*a(n-3)+...
# Ignoring some of the earlier coefficients is simply achieved by
# specifying a 0 at their positions. Note that the leftmost element
# in the list recurL is attached to Taylor coefficients just before
# the element a(n), the next one to the element even one index earlier
# etc.
# off: an integer offset; if larger than zero, the first off elements
# of the Taylor expansion of the generating function (starting at the
# constant term x^0) become zero.
# ogf: optional additional inhomogeneity at the right hand side of the
# recurrence, a generating function in the variable x.
# verb: if set to a value larger than zero, the number of expansion coefficients
# printed.
# pown: optional additional inhomogeneity at the right hand side of the
# recurrence, a sum of n-th powers of constants.
# The format is the same as for the pown argument of GenPolyn().
#
# Example 1:
# GenFLinRecG( [1], [4], 0, (1+x+x^2+x^3)/(1-x)/(1-x^3), 10,[]) returns the generating function
# specified by a(0)=1, a(n)=4*a(n-1)+b(n), see OEIS A141844 and A042968.
#
# @since 2008-09-06
GenFLinRecG := proc(initL, recurL, off, ogf, verb, pown)
local n, # numerator
d, # denominator
i,j,g ;
n := 0 ;
for i from 1 to nops(initL) do
n := n+(op(i,initL)-coeftayl(ogf,x=0,i-1))*x^(i-1+off) ;
end do:
for i from 1 to nops(recurL) do
for j from 1 to nops(initL)-i do
n := n-op(i,recurL)*op(j,initL) * x^(i+j-1+off) ;
end do:
end do:
n := n+ ogf ;
d := 1 ;
for i from 1 to nops(recurL) do
d := d-op(i,recurL)*x^i ;
end do:
g := factor(n/d) ;
if verb > 0 then
GenFVerify(convert(g,parfrac,x),verb) ;
GenFpurrs(g) ;
end if ;
RETURN(g) ;
end proc:
# Compute a generating function with expansion coefficients
# obtained by shuffling the expansion coefficients of two or more
# other generating functions.
#
# arguments: genFList is a list of k generating functions with main
# variable x with
# expansion gf(1)=a(1,0)+a(1,1)*x+a(1,2)*x^2+a(1,3)*x^3+.. for the
# first, gf(n)=a(n,0)+a(n,1)*x+a(n,2)*x^2+... for the nth, 1<=n<=k.
#
# returns: the generating function equivalent to the expansion
# a(1,0)+a(1,1)*x^k+a(1,2)*x^(2k)+a(1,3)*x^(3k)
# +...+a(n,0)*x^(n-1)+a(n,1)*x^(n-1+k)+a(n,3)*x^(n-1+2k)+...
# equivalent to fanning out the first sequence by creating
# space for intermediate powers, and to create the output by
# insertion of the g(2) with places in the coefficients right after
# those of g(1), g(3) with places right after those of g(2) etc.
#
# Example 1:
# GenFShuff([2*exp(x),1/(1-x)]) returns a generating function
# (-2*exp(x^2)+2*exp(x^2)*x^2-x)/(x-1)/(x+1), where the even
# indices in the Taylor expansion are from the coefficients of
# 2*exp(x)=2+2x+x^2+x^3/3+.., the odd indices from 1/(1-x)
# =1+x+x^2+x^3+x^4+.. So the function returned has the
# expansion 2+1*x+2*x^2+1*x^3+x^4+x^5+x^6/3+x^7+...
#
# Example 2:
# GenFShuff([1/(1+x),x/(1-x-2*x^2)],true) returns the generating
# function for A133730, mixing the generating function of A033999
# with the generating function of A001045.
#
# @since 2007-12-01
GenFShuff := proc(genFList,verb)
local gf,i,k ;
gf := 0 ;
k := nops(genFList) ;
for i from 1 to k do
gf := gf+x^(i-1)*subs(x=x^k,op(i,genFList)) ;
end do:
gf := factor(gf) ;
if verb >0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Return two generating functions obtained by bi-section of the
# coefficients of the generating function provided.
# argument:
# genFmerg is the original function a0+a1*x+a2*x^2+...
#
# returns:
# two generating functions, the first representing the series
# a0+a2*x+a4*x^2+a6*x^3+..., the second a1+a3*x+a5*x^2+a7*x^3+...
#
# Example 1:
# First create a rational polynomial 2+x+3x^2+4x^3+7x^4+11*x^5+..
# luc := (2-x)/(1-x-x^2) ; # Lucas numbers, OEIS A000032
# Then assign (1+x)/(1-3x+x^2)= 1+4*x+11*x^2+... to biluc:
# biluc := op(2,GenFDis(luc)) ;
# taylor(biluc,x=0,14) ; # essentially OEIS A093960
#
# Example 2:
# First dissect the generating function of A024495
# a024495 := GenFDis( x^2/((1-x)^3-x^3) ) ;
# Then create the generating function of A135574
# GenFShuff([op(2,a024495),op(1,a024495)],true) ;
#
# @since 2007-12-01
GenFDis := proc(genFmerg)
local gminx,geven,godd ;
# seperate the input function into the even and odd part
# with respect to the variable x.
gminx := subs(x=-x,genFmerg) ;
geven := simplify((genFmerg+gminx)/2) ;
godd := simplify((genFmerg-gminx)/2) ;
geven := subs(x=sqrt(x),expand(numer(geven)))/subs(x=sqrt(x),expand(denom(geven))) ;
godd := factor(godd)/x ;
godd := subs(x=sqrt(x),expand(numer(godd)))/subs(x=sqrt(x),expand(denom(godd))) ;
geven := numer(geven)/expand(denom(geven)) ;
godd := numer(godd)/expand(denom(godd)) ;
RETURN([factor(simplify(geven)),factor(simplify(godd))]) ;
end proc:
# Return the generating function which has the terms from a(0)
# over a(1)*x up to a(off-1)*x^(off-1) removed from the input.
# argument:
# genf the input generating function with typically some
# nonzero Taylor expasion coefficients for powers 0 to off-1.
#
# @since 2007-12-01
GenFDelToOffs := proc(genf,off)
local p,gf;
gf := genf ;
for p from 0 to off-1 do
gf := gf- coeftayl(gf,x=0,p)*x^p
end do:
RETURN( factor(gf)) ;
end proc:
# Return the generating functions obtained by using only each
# m-th coefficient of list from the generating function provided.
# argument:
# genf is the original function a(0)+a(1)*x+a(2)*x^2+...
# with x being the independent variable.
#
# returns:
# the function equivalent to a(0+c)+a(m+c)*x+a(2m+c)*x^2+a(3m+c)*x^3+..
# again with x the independent variable, a multisection of the
# original series.
#
# The method is described in J. Riordan's book "Combinatorial Identities"
# (John Wiley, 1968) Section 4.3
#
# Example 1:
# First create a rational polynomial 2+x+3x^2+4x^3+7x^4+11*x^5+18*x^6+..
# luc := (2-x)/(1-x-x^2) ; # Lucas numbers, OEIS A000032
# Then return (2-4x)/(1-4x-x^2) = 2+4x+18x^2+..., keeping each 3rd
# of the coefficients above, essentially creating OEIS A014448
# GenFm(luc,3,0) ;
#
# Example 2:
# Compute the generating function of A132397 given the generating
# function of A024494:
# GenFm(x*(1-x)/(1-2*x)/(1-x+x^2),3,2) ;
#
# @since 2007-12-01
GenFm := proc(genf,m,c)
local diss,k;
if c < 0 then
error "negative third argument in GenFm"
elif c > 0 then
diss := GenFDelToOffs(genf,c)/x^c ;
else
diss := genf ;
end if ;
if m = 1 then
RETURN(diss) ;
elif m = 2 then
RETURN(op(1,GenFDis(diss))) ;
elif m <= 0 then
error "zero or negative second argument in GenFm"
elif 2^ilog2(m) = m then
diss := op(1,GenFDis(diss)) ;
RETURN( GenFm(diss,m/2) ) ;
else
diss := add( subs(x=exp(2*Pi*I*k/m)*x^(1/m),diss),k=0..m-1)/m ;
diss := factor(diss) ;
diss := numer(diss)/expand(denom(diss)) ;
diss := factor(diss) ;
end if:
end proc:
# Generating function for a m-st order difference.
# arguments:
# genf is the original generating function, a(0)+a(1)*x+a(2)*x^2+..
# with x the independent variable.
# m is the non-negative order of the difference scheme
#
# returns:
# the generating function which has the m-st order differences (essentially
# a linear recurrence with alternatingly signed binomial coefficients) as
# Taylor expansion coefficients. For the 1st order differences this
# is just a multiplication of genf by 1-x on return.
#
# Example:
# For m=1, the generating function returned is a(1)-a(0)+[a(2)-a(1)]*x+..
#
# @since 2007-12-02
GenFmOrdDiff := proc(genf,m,verb)
local gf,i;
gf := genf ;
for i from 1 to m do
gf := (1-x)*gf ;
gf := GenFDelToOffs(gf,1) ;
gf := factor(gf/x) ;
end do:
if verb > 0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Generating function given a m-st order difference
# This is the inverse functionality of GenFmOrdDiff().
# arguments:
# genf is the generating function for the m-th order difference
#
# a0L is the list of leading terms of the differences. a0L[1] is the
# a(0) coefficient of the (m-1)st order generating function, a0L[2]
# the a(0) coefficient of the (m-2)nd order generating function. The
# length of this list implies the order of the differences.
#
# verb turns on an explicit display of the first terms of the results
#
# off is an integer which (if non-zero) multiplies by powers of x
# to translate the indices
#
# returns:
# the generating function which has the m-st order differences as provided.
# For m=1, this computes the generating function of the partial
# sums of the initial genf, just a division of genf through 1-x.
#
# Example:
# First generate the generating function for the periodic 1,0,0,1,
# which yields per=(x^2-x+1)/(1-x)/(x^2+1).
# per := GenFLinRec([1,0,0,1],[0,0,0,1],0 ,0, false, []) ;
# Then compute the generating function which has per as the first
# differences and the leading coefficient 3:
# GenFMOrdDiffInv(per,[3],0,true) ;
# which is 3+4x+4x^2+4x^3+5x^4+6x^5+6x^6+6x^7+7x^8+...
#
# @since 2007-12-05
GenFmOrdDiffInv := proc(genf,a0L,off,verb)
local gf,i;
gf := genf ;
for i from 1 to nops(a0L) do
gf := x*gf+op(i,a0L) ;
gf := factor(gf/(1-x)) ;
end do:
gf := factor(x^off*gf) ;
if verb >0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Generating function of the Motzkin Transform.
# arguments:
# genf is the generating function for the original sequence
#
# verb turns on an explicit display of the first terms of the results
# if larger than zero.
#
# returns:
# the generating function of the Motzkin transform of the
# original series.
#
# @since 2008-11-08
GenFMotzkinT := proc(genf,verb)
local gf,i;
gf := subs(x= (1-x-sqrt(1-2*x-3*x^2))/2/x,genf) ;
gf := simplify(gf) ;
if verb >0 then
GenFVerify(gf,verb) ;
end if ;
RETURN(gf) ;
end proc:
# Generating function of the Inverse Motzkin Transform.
# arguments:
# genf is the generating function for the transformed sequence
#
# verb turns on an explicit display of the first terms of the results
# if larger than zero.
#
# returns:
# the generating function of the original series.
#
# @since 2008-11-08
GenFMotzkinTInv := proc(genf,verb)
local gf,i;
gf := subs(x= x/(1+x+x^2),genf) ;
gf := factor(gf) ;
if verb >0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Generating function of the L-Schroeder Transform.
# arguments:
# genf is the generating function for the original sequence
#
# verb turns on an explicit display of the first terms of the results
# if larger than zero.
#
# returns:
# the generating function of the Large-Schroeder transform of the
# original series.
#
# @since 2008-11-13
GenFSchroederT := proc(genf,verb)
local gf,i;
gf := subs(x= (1-x-sqrt(1-6*x+x^2))/2,genf) ;
gf := factor(expand(gf)) ;
if verb >0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Generating function of the inverse L-Schroeder Transform.
# arguments:
# genf is the generating function for the original sequence
#
# verb turns on an explicit display of the first terms of the results
# if larger than zero.
#
# returns:
# the generating function of the inverse Large-Schroeder transform of the
# original series.
#
# @since 2008-11-13
GenFSchroederTInv := proc(genf,verb)
local gf,i;
gf := subs(x= x*(1-x)/(1+x),genf) ;
gf := factor(expand(gf)) ;
if verb >0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Generating function of the Catalan Transform.
# arguments:
# genf is the generating function for the original sequence
#
# verb turns on an explicit display of the first terms of the results
# if larger than zero.
#
# returns:
# the generating function of the Catalan transform of the
# original series.
#
# @since 2008-11-06
# @see P. Barry, J. Int. Seq. vol 8 (2005) #05.4.4
GenFCatalanT := proc(genf,verb)
local gf,i;
gf := subs(x= (1-sqrt(1-4*x))/2,genf) ;
gf := factor(expand(gf)) ;
if verb >0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Generating function of the Inverse Catalan Transform.
# arguments:
# genf is the generating function for the original sequence
#
# verb turns on an explicit display of the first terms of the results
# if larger than zero.
#
# returns:
# the generating function of the Inverse Catalan transform of the
# original series.
#
# @since 2008-11-06
# @see P. Barry, J. Int. Seq. vol 8 (2005) #05.4.4
GenFCatalanTInv := proc(genf,verb)
local gf,i;
gf := subs(x= x*(1-x),genf) ;
gf := factor(gf) ;
if verb >0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Generating function of Inverse Binomial Transform.
# arguments:
# genf is the generating function for the original sequence
#
# verb turns on an explicit display of the first terms of the results
# if larger than zero.
#
# returns:
# the generating function of the inverse binomial transform of the
# original series.
#
#
# @since 2008-04-02
GenFBinomialTInv := proc(genf,verb)
local gf,i;
gf := subs(x= x/(1+x),genf)/(1+x) ;
gf := factor(gf) ;
if verb >0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Generating function of Binomial Transform.
# arguments:
# genf is the generating function for the original sequence
#
# verb turns on an explicit display of the first terms of the results of larger than zero
#
# returns:
# the generating function of the binomial transform of the
# original series.
#
# Example:
# First compute the g.f. of the periodic sequence bar(2,3,1) with period 3,
# which is -(x+2)(x+1)/(x-1)/(1+x+x^2):
# gf := GenFLinRec( [2,3,1],[0,0,1],0,0,false,[]) ;
# Then compute the generating function of its inverse binomial transform, as done
# in A130752:
# GenFBinomialTInv(gf,true) ;
#
# @since 2008-04-02
GenFBinomialT := proc(genf,verb)
local gf,i;
gf := subs(x= x/(1-x),genf)/(1-x) ;
gf := factor(gf) ;
if verb > 0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Generating function of Euler Transform.
# arguments:
# genf is the generating function for the original sequence
#
# verb turns on an explicit display of the first terms of the results of larger than zero
#
# returns:
# the generating function of the Euler transform of the
# original series.
#
# \begin{latexonly}
# \cite{FlajoletTCS144}
# \end{latexonly}
# @since 2009-03-07
GenFEulerT := proc(genf,verb)
local gf,i;
gf := subs(x= -x/(1-x),genf)/(1-x) ;
gf := factor(gf) ;
if verb > 0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Print a power representation of the coefficients of the generating function.
# @param genf the generating function as a function of x
# @since 2008-04-16
GenFpurrs1Term := proc(genf)
local g,res,den,nu,c,denPr,xroo,tt,ctt,dentt,nutt ;
res := 0 ;
if degree(denom(genf),x) = 0 then
res := res+ genf*Dirac(n);
elif degree( denom(genf),x) = 1 then
den := denom(genf) ; nu := numer(genf) ;
c := [coeftayl(den,x=0,0),-coeftayl(den,x=0,1)] ;
if op(1,c) = 0 then
RETURN() ;
end if ;
res := res+ nu*(op(2,c)/op(1,c))^n/op(1,c) ;
elif degree( denom(genf),x) = 2 then
den := denom(genf) ; nu := numer(genf) ;
denPr := coeftayl(den,x=0,2) ;
xroo := solve(den,x) ;
if xroo[1] <> xroo[2] then
xroo := convert(nu/denPr/(x-xroo[1])/(x-xroo[2]),parfrac,x) ;
for tt in xroo do
dentt := denom(tt) ; nutt := numer(tt) ;
ctt := [coeftayl(dentt,x=0,0),-coeftayl(dentt,x=0,1)] ;
res := res+ nutt*(op(2,ctt)/op(1,ctt))^n/op(1,ctt) ;
end do:
else
# term format is nu/[denPr*(x-xroo[1])^2]
res := res+ nu/denPr*(n+1)*(1/xroo[1])^(n+2) ;
fi ;
else
RETURN() ;
end if ;
RETURN(res) ;
end proc:
# Print a power representation of the coefficients of the generating function.
# @param genf the generating function as a function of x
#
# Example 1 (A108981):
# GenFpurrs( -(1+2*x)/(1+x)/(4*x-1)) ;
# prints -(-1)^n/5+6*4^n/5
#
# Example 2 (A003948):
# GenFpurrs( (1+x)/(1-5*x) ) ;
# prints -(-1)^n/5+6*4^n/5
# @since 2008-04-16
GenFpurrs := proc(genf)
local g,res,t,tt ;
# maximum polynomial order of x in the denominators of the fractions
if type(genF, ratpoly) = true then
g := convert(genf,parfrac,x) ;
res := 0 ;
if type(g,`+`) then
for t in g do
tt := GenFpurrs1Term(t) ;
if whattype(tt) = exprseq then
RETURN() ;
else
res := res+tt ;
end if ;
end do;
else
tt := GenFpurrs1Term(g) ;
if whattype(tt) = exprseq then
RETURN() ;
else
res := tt ;
end if ;
end if ;
print("test purrs") ;
res := simplify(res,symbolic) ;
print(factor(res)) ;
print(res) ;
print(taylor(genf,x=0,15)) ;
for t from 0 to 15 do
printf("%a,", simplify(subs(n=t,res)) ) ;
end do:
printf("\n") ;
for t from 0 to 10 do
printf("%a,", evalf(simplify(subs(n=t,res))) ) ;
end do:
printf("\n") ;
end if ;
return ;
end proc:
# r-th Witt transform
#
# Example 1 (A000012 -> A110654):
# GenFWitt( 1/(1-x),2,20) ;
# Example 2 (A000027 -> A006584):
# GenFWitt( x/(1-x)^2,2,20) ;
# Example 3 (A000045 -> A089089):
# GenFWitt( x/(1-x-x^2),2,20) ;
# Example 4 (A000217 -> A032092):
# GenFWitt( 1/(1-x)^3,2,20) ;
# Example 5 (A000292 -> A032094):
# GenFWitt( 1/(1-x)^4,2,20) ;
# Example 6 (A000012 -> A006918):
# GenFWitt( 1/(1-x),4,20) ;
# Example 7 (A000027 -> 2*A031164):
# GenFWitt( 1/(1-x)^2,4,20) ;
# Example 8 (A000045 -> A089116):
# GenFWitt( 1/(1-x-x^2),3,20) ;
#
# @param genf the generating function as a function of x
# @param r the order of the Witt transform
# @param verb highest number of the taylor series coefficients printed for the
# transformed sequence.
# @since 2008-07-27
# @see Pieter Moree, The formal series Witt transform, Discr Math 295 (2005) 143-160
GenFWitt := proc(genf,r,verb)
local gf,d;
gf := 0 ;
for d in numtheory[divisors](r) do
gf := gf + numtheory[mobius](d)*(subs(x= x^d,genf))^(r/d) ;
end do:
gf := factor(gf/r) ;
if verb > 0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# r-th Akiyama-Tanigawa transform
# @param genf the generating function of the original function
# @param r the order of the transform
# @param verb the number of taylor series coefficients to be printed
# @see D. Merlini, R. Sprugnoli, M. Cecilia Verri, INTEGERS: El. J. Combinat. Numb.
# Theory vol 5 (2005), #A05
# @since 2008-07-30
GenFAkiyama := proc(genf,r,verb)
local gf,k;
gf := genf ;
for k from 0 to r do
gf := 1-(1-x)*diff(gf,x) ;
end do:
if verb > 0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Invert transform.
# Supposed the input is of the form genf=sum_{i=0} b_i x^i with b_0=1,
# the output is 1/(1-sum_{i=1} b_ix^i).
# Note that generally sequences need to be shifted before using it as input
# to be compatible with the transforms.txt implementation of the OEIS, which
# assumes that the first entry of the list is the term in front of x^1 of the
# generating function.
#
# @param genf the generating function as a function of x
# @param verb highest number of the taylor series coefficients printed for the
# transformed sequence.
# @since 2008-07-28
#
# Example: (A004227 to A098182)
# g:= (1+x^2)/(1-x)^2 ; # as in A004227
# GenFInvrt(1+g*x,20) ; # similar to A098182
#
GenFInvrt := proc(genf,verb)
local gf;
gf := 0 ;
if coeftayl(genf,x=0,0) = 1 then
gf := factor((1/(2-genf)-1)/x) ;
end if:
if verb > 0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
end if ;
RETURN(gf) ;
end proc:
# Guess a linear recurrence.
# Example A086515:
# GenFGuess([1,3,7,9,13,15,19,21,25,27,31,33,37,39,43,45,49,51],3,0,20) ;
#
# @param L the list of the first terms of the sequence
# @param o the order of the recurrence
# @param offs offset assigned to the first term
# @param verb the number of terms printed
# @return a generating function for the sequence if all terms follow a common
# linear recurrence with the same constant coefficients
# @since 2008-09-06
# \beginlatexonly
# \cite{MansourAJC30}
# \endlatexonly
GenFGuess := proc(L,o,offs,verb)
local f,l,A,r,c,B,s,sold,consi,g;
l := nops(L) ;
A := Matrix(o,o) ;
B := Vector(o) ;
consi := true ;
for f from 1 to l-2*o+1 do
for r from 1 to o do
for c from 1 to o do
A[r,c] := op(f+r+c-2,L) ;
end do:
B[r] := op(f+r+o-1,L) ;
end do:
s := LinearAlgebra[LinearSolve](A,B) ;
if verb > 0 then
printf("a(n)=") ;
for c from 1 to o do
if s[o+1-c] = 1 then
printf(" +a(n-%d)",c) ;
elif s[o+1-c] = -1 then
printf(" -a(n-%d)",c) ;
elif s[o+1-c] < 0 then
printf(" %a*a(n-%d)",s[o+1-c],c) ;
elif s[o+1-c] > 0 then
printf(" +%a*a(n-%d)",s[o+1-c],c) ;
fi;
end do:
printf("\n") ;
end if;
if f >1 and not LinearAlgebra[Equal](s,sold) then
consi := false ;
end if;
sold := s;
end do:
if consi then
g := GenFLinRec( [op(1..o,L)], [seq(s[o+1-i],i=1..o)], offs, 0, verb,[]) ;
if verb > 0 then
printf("G.f. ( %a ) / ( %a ). - _R. J. Mathar_, %s\n", factor(numer(g)),factor(denom(g)),StringTools[FormatTime]("%b %d %Y")) ;
end if;
return g ;
else
return () ;
end if;
end proc:
# Generate a polynomial representation from terms of the form 1/(x-1)^k+1/(x+1)^k2
# @param L a list of the form [[n1,k1],[n2,k2],....] where
# each component represents a term nj/(x-1)^kj in a generating function.
# @param Lplus a list of the form [[n1,k1],[n2,k2],....] where
# each component represents a term nj/(x+1)^kj in a generating function.
# @return a polynomial in n representing the sum of these two generating functions
# @since 2008-09-07
GenF2pol := proc(L,Lplus,verb)
local p,t,k,nu,j;
p := 0 ;
for t in L do
k := op(2,t) ;
nu := op(1,t) ;
if k = 1 then
p := p-nu ;
elif k = 2 then
p := p+nu*(n+1) ;
elif k = 3 then
p := p-nu*(n+1)*(n+2)/2 ;
else
p := p+(-1)^k*nu/(k-1)!*mul(n+j,j=1..k-1) ;
end if;
end do:
for t in Lplus do
k := op(2,t) ;
nu := op(1,t) ;
if k = 1 then
p := p+(-1)^n*nu ;
else
p := p+(-1)^n*nu/(k-1)!*mul(n+j,j=1..k-1) ;
end if;
end do:
p := expand(p) ;
if verb > 0 then
for k from 0 to verb do
printf("%a,",subs(n=k,p)) ;
end do;
printf("\n") ;
end if;
RETURN(p) ;
end proc:
# Generate a polynomial representation from a list of equally-spaced abscissa points.
# @param L a list of the form [f[o],f[o+1],..] with the sampled points
# @param o the offset, declaring the index of L[1].
# @param po the order of the polynomial
# @return a polynomial in n
# @since 2011-02-10
L2pol := proc(L,o,po,verb)
# Lagrange interpolation formula
local p,i,lnumer,j;
p := 0 ;
if po+1 > nops(L) then
error "L2pol: insufficient number of terms for polynomial order "
end if;
for i from 0 to po do
lnumer := 1 ;
for j from 0 to po do
if j <> i then
lnumer := lnumer*(n-(j+o))/(i-j) ;
end if;
end do:
p := p+lnumer*op(i+1,L) ;
p := expand(p) ;
end do;
if verb > 0 then
for i from 0 to verb do
printf("%a,",subs(n=i,p)) ;
end do;
printf("\n") ;
end if;
return factor(p);
end proc:
# Generating function of an integer polynomial modulo a constant.
# @param poln the polynomial with n as the independent variable
# @param per the period of the modulo-evaluation
# @param ceil if true, return the ceiling function
# @param verb the number of terms to be evaluated for verification
# @return the generating function for the sequence of that mod function
# p(n) mod per, n=0,1,2,3...
#
# Example:
# GenFPmodk(n,10,false,8) ;
# returns the generating function of n mod 10.
# @since 2008-09-09
GenFPmodk := proc(poln,per,ceil,verb)
local gf,resi,i,perio,t,lcD,pol ;
# work with polynomials with integer coefficients
lcD := 1 ;
for i from 0 to degree(poln) do
lcD := lcm(lcD, denom(coeff(poln,x,i)) ) ;
od:
pol := poln*lcD ;
# add the generating function of the per-periodic remainder
if ceil then
resi := [] ;
for i from 0 to per*lcD-1 do
t := subs(n=i,pol) mod (per*lcD);
if t <> 0 then
t := lcd*per-t ;
end if;
resi := [op(resi),t] ;
end do:
else
resi := [seq(subs(n=i,pol) mod (per*lcD),i=0..lcD*per-1)] ;
end if ;
perio := [seq(0,i=1..lcD*per-1),1] ;
gf := GenFLinRec( resi, perio, 0, 0, 0,[]) ;
gf := factor(gf/lcD) ;
if verb > 0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
GenFgf2Rec(gf) ;
end if ;
return gf ;
end proc:
# Convert a rational polynomial generating function in x into a linear recurrence.
# This prints the denominator of the polynomial in the style of a(n) = sum_i ... a(n-i).
# @param gf the rational polynomial in x
# @since 2010-04-15
GenFgf2Rec := proc(gf)
local g,a0,d ;
g := denom(factor(gf)) ;
# ensure that this is a monic polynomial
a0 := simplify(subs(x=0,g)) ;
g := expand(g/a0) ;
# print(g) ; # debugging
if degree(g,x) < 11 then
# printf("Index to sequences with linear recurrences with constant coefficients, signature (") ;
printf("Index to sequences with linear recurrences with constant coefficients, signature (",degree(g,x)) ;
for d from 1 to degree(g,x) do
a0 := -coeff(g,x,d) ;
printf("%d",a0) ;
if d <> degree(g,x) then
printf(",") ;
end if;
end do ;
printf(").\n") ;
end if;
printf("a(n) =") ;
for d from 1 to degree(g,x) do
a0 := -coeff(g,x,d) ;
if a0 = 1 then
printf(" +a(n-%d)",d) ;
elif a0 = -1 then
printf(" -a(n-%d)",d) ;
elif a0 < 0 then
printf(" %d*a(n-%d)",a0,d) ;
elif a0 > 0 then
printf(" +%d*a(n-%d)",a0,d) ;
end if;
end do ;
printf(".\n") ;
return ;
end proc:
# Generating function of the floor an integer polynomial divided by a constant.
# @param poln the polynomial with n as the independent variable
# @param den the denominator
# @param ceil if true, work with ceil() instead of floor().
# @param offs offset to remove initial coefficients
# @param verb the number of terms to be evaluated for verification
# @return the generating function for the sequence of that floor function
# floor(p(n)/den)
#
# Example 1:
# GenFFloorP(3*n,4,false,0,10) ;
# returns the generating function of floor(3n/4), A057353
# Example 2:
# GenFFloorP(FallFactInv([[1,5]]),5!*5,false,0,10) ;
# returns the generating function of floor(binomial(n,5)/5), A011851
# @since 2008-09-09
GenFFloorP := proc(poln,den,ceil,offs,verb)
local gf,resi,i,per,pol ;
# since poln may contain binomial() expressions,
# we can cover these by expanding first.
pol := expand(poln) ;
# first the generating function of the polynomial
gf := GenPolyn(pol,0,0,[]) ;
if ceil then
# decompose floor(p(n)) = p(n) - p(n) mod den
gf := gf+ GenFPmodk(pol,den,ceil,0) ;
else
# decompose floor(p(n)) = p(n) - p(n) mod den
gf := gf- GenFPmodk(pol,den,ceil,0) ;
end if;
gf := GenFDelToOffs(gf,offs) ;
gf := factor(gf/den) ;
if verb > 0 then
GenFVerify(gf,verb) ;
GenFpurrs(gf) ;
GenFgf2Rec(gf) ;
end if ;
return gf ;
end proc: