TRANSFORMATIONS OF INTEGER SEQUENCES N. J. A. Sloane This is plain text file giving a number of Maple procedures for performing transformations on sequences and numbers. This is a subpage of the On-Line Encyclopedia of Integer Sequences (see http://www.research.att.com/~njas/sequences/Seis.html) which makes extensive use of these transformations. The Maple Procedures Instructions for downloading these procedures: Download this file, strip off everything above "#### 1 ####"" (about 8 lines below), store the result in a file called transforms. It can then be read directly into Maple by saying (in Maple) read transforms; #### 1 #### # transforms # Jul 01 2004: I had to make many changes to get this to run # correctly under the new Maple (9.5). I really love how # Maple makes changes so that programs written in old # versions never run under the new version. # added LAH and LAHi, Oct 05 2003 # added RECORDS Sep 17 2002 # added SHADOW Aug 02 2002 # added CONTINUANTi Jul 18, 2002 # added transforms from Antti Karttunen, Jul 12 2002 # added HANKEL, lHANKEL, CONTINUANT, Jun 10, 2002 # added dirsort Jan 19 2002 # Dec 24 2001: signs corrected in REVERT and REVEGF transforms. # Modified Jun 17 2001 # A collection of Maple procedures for transforming sequences # These mostly work on lists. Example: # a:=[1,1,2,3,5,8,13,21]; # # [1, 1, 2, 3, 5, 8, 13, 21] # # BINOMIAL(a); # # [1, 2, 5, 13, 34, 89, 233, 610] # In case of trouble they return the empty list [], not an error message. # Requires that you have already done: with(numtheory, mobius): # with(combinat, stirling1, stirling2 ), and loaded gfun: with(numtheory, mobius): with(combinat, stirling1, stirling2 ): # for maple6 i changed the following line on 4/11/00 from #with(share): with(gfun): # to with(gfun): # (actually the gfun package isn't used much, and # most of these transforms will work without it) # Reference : M. Bernstein and NJAS, Some Canonical Sequences of Integers, # Linear Algebra and its Applications 1995 volume 226-228 pages 57-72 # Available via my home page # http://www.research.att.com/~njas/doc/eigen.ps or eigen.pdf # Summary: # ALLDIFFS all differences a(j)-a(i), j>i # ANDCONV AND-convolution, SUM a(k).AND.a(n-k) # BINOMIAL from 1st diag of diff table to top row # BINOMIALi inverse: from sequence to leading diag of diff table # BIN1 a variant of BINOMIAL used by Zagier # BISECT(a,0), BISECT(a,1) bisect a sequence # BOUS boustrophedon transform # BOUS2 boustrophedon transform (official boust. transform) # BOUS2i inv. boust transform (official boust. transform) # CHAR characteristic function of sequence # COMP complement of sequence # COMPl complement of sequence (long version) # COMPOSE compose two sequences # COMPSQRT functional square root # CONTINUANT continuant transform, cf. Concrete Math. p. 301 # CONTINUANTi inverse continuant transform (not always integral) # CONV simple convolution, expand A(x)*B(x) # CONVi square root of convolution, not always integral # DECIMATE(a,k,0), DECIMATE(a,k,1), ... decimate a sequence # DIFF first difference # DIGREV reverse digits (use digrev for single numbers) # DIGSUM sum of digits (use digsum for single numbers) # DIRICHLET Dirichlet convolution of two sequences # EULER Euler Xfm # EULERi inverse Euler Xfm # EXP egf of b = exp (egf of a) # EXPCONV exponential convolution, expand E1(x)*E2(x) # EXTRACT extract subsequence # GCDCONV GCD-convolution, SUM a(k).gcd.a(n-k) # RECORDS extracts records from a list # HANKEL Hankel transform # HEAP Given a, start with heap of 1, add largest entry in a <= heap to heap; b gives successive sizes of heaps # INVERT a's from b's in 1+SUM a_i x^i = 1/(1-SUM b_i x^i) # INVERTi inverse, b's from a's # LAH egf of b = egf of a at x/(1-x) # LAHi egf of b = egf of a at x/(1+x) # LCMCONV LCM-convolution, SUM a(k).lcm.a(n-k) # LEFT shift left # lHANKEL little Hankel transform # LISTTOLISTDIV divides nth term by n! # LISTTOLISTMULT multiplies nth term by n! # LOG egf of b = log (egf of a) # MASKCONV mask-convolution # MASKTRANS mask-convolution with all-1's sequence # MASKTRANSi inverse # MOBIUS Mobius (or Lambert) transform # MOBIUSi inverse Mobius (or sum of divisors) transform # MONO sort, discard duplicates # MONO2 sort, discard duplicates, but only if difft from orig. # M2 multiply all except 1st term by 2 # M2i divide all except 1st term by 2 # NEGATE negate all except 1st term # ORCONV OR-convolution, SUM a(k).OR.a(n-k) # PARTITION partition Xfm (without repetition) # PARTITIONi inverse partition Xfm (without repetition) # POLY get polynomial for seq from difference table # POWERS extracts powers of x in f through degree n # PRODS sequence of partial products # PSUM partial sums # PSUMSIGN alternating sign partial sums # REVERT reversion # REVEGF reversion of e.g.f. # RIGHT shift right, prepending a 1 # SERIESTOLISTDIV divides nth term by n! # SERIESTOLISTMULT multiplies nth term by n! # SERIESTOSERIESDIV divides nth term by n! # SERIESTOSERIESMULT multiplies nth term by n! # SERIES2 series expansion of function of 2 variables # SERIES2TOLIST converts output of SERIES2 to a list, order 00,10,01,20,11,02,... # SERIES2TOLISTMULT ditto, multiplying coefft of x^i*y^j by i!*j! # SERIES2HTOLIST converts output of SERIES2 to a list, order 00,10,11,20,21,22,... # SERIES2HTOLISTMULT ditto, multiplying coefft of x^i*y^j by i!*j! # SHADOW see Lorenz Halbeisen and Norbert Hungerbuehler reference in sequence database under A000522 # STIRLING Stirling Xfm (of 2nd kind) # STIRLINGi inverse Stirling Xfm (equivalently, Stirling transform of 1st kind) # STIRB Stirling-Bernoulli transform # SUPPORT positions where list is nonzero # TRISECT(a,0), .. TRISECT(a,2) trisect a sequence # WEIGH b from a in 1+SUM b_n x^n=PI (1+x^a_n) # XORCONV XOR-convolution, SUM a(k).XOR.a(n-k) # others: # delta delta(f,n,k) = kth term of nth difference of f # did did it divide? # dids did it divide (with sign)? # digrev reverse digits # digsort sort digits # digsum sum of digits # dim did it mask? # Etrans Euler Xfm (again) # mex minimum excluded number # nimsum nim sum # ptrans partition Xfm (with repetition) # traptrans inverse partition Xfm (with repetition) # weighout b from a in 1+SUM b_n x^n=PI (1+x^n)^a_n # weighouti a from b in 1+SUM b_n x^n=PI (1+x^n)^a_n # weighini a from b in 1+SUM b_n x^n=PI (1+x^a_n) # weigh2out b from a in 1+SUM b_n x^n=PI (x^-n+1+x^n)^a_n # weigh2outi a from b in 1+SUM b_n x^n=PI (x^-n+1+x^n)^a_n # weigh2in b from a in 1+SUM b_n x^n=PI (x^-a_n+1+x^a_n) # weigh2ini a from b in 1+SUM b_n x^n=PI (x^-a_n+1+x^a_n) # eultrans2 2nd Euler trans from paper by Donaghey and Shapiro # but this is same as BINOMIALi # pairtrans b(n)=a(n)+a(n-1) # pairtransi a(n)=b(n)-b(n-1)+b(n-2)-... # wt Compute weight or number of 1's in binary expansion of n: # ANDnos logical AND of two numbers using their binary expansions # ORnos logical OR of two numbers using their binary expansions # XORnos logical XOR of two numbers using their binary expansions # deriv derivative of a numbers using its binary expansion delta:=proc(f,n,k) local i; add((-1)^(n-i)*binomial(n,i)*f(k+i),i=0..n); end: did:=proc(m,n) if irem(m,n) = 0 then 1 else 0: fi: end: dids:=proc(m,n) if irem(m,n) = 0 then (-1)^(m/n) else 0: fi: end: # Does it Mask (j over i)? dim := proc(j,i) if(ANDnos(j,i) = i) then 1 else 0; fi; end: mob:=proc(m,n) if irem(m,n) = 0 then mobius(m/n) else 0: fi: end: # Reference: Marc LeBrun's (mlb@well.com) message # "half-baked generalized convolution sequence transforms" # posted 09 Jun 2001 on SeqFan mailing list (seqfan@ext.jussieu.fr) # URL: http://www.ccr.jussieu.fr/gmpib/seqfan/seqfan.html # Also http://www.megabaud.fi/~karttu/matikka/findnext.txt MASKCONV :=proc(a,b) local c,i,j,n; if whattype(a) <> list then RETURN([]); fi; if whattype(b) <> list then RETURN([]); fi; n := min(nops(a),nops(b)); c := []; for j from 0 to n-1 do c := [ op(c), add( dim(j,i)*a[i+1]*b[(j-i)+1], i=0..j)]; od; RETURN(c); end: # MASKTRANS(a) is equivalent to MASKCONV(a,A000012), where A000012 is # the all-1's sequence. MASKTRANS :=proc(a) local c,i,j,n; if whattype(a) <> list then RETURN([]); fi; n := nops(a); c := []; for j from 0 to n-1 do c := [ op(c), add( dim(j,i)*a[i+1], i=0..j)]; od; RETURN(c); end: MASKTRANSi := a -> MASKCONV(a,[seq((-1)^(wt(j) mod 2),j=0..nops(a)-1)]): # HEAP transform. Assumes a is monotonic and begins with 0 or 1 HEAP:=proc(a) local pttemp,pt,la,lb,MAXLEN,heap,b,i,k: if whattype(a) <> list then RETURN([]); fi: if a[1] < 0 or a[1] > 1 then RETURN([]); fi: la:=nops(a); lb:=1; heap:=1; b:=[heap]: pt:=1; MAXLEN:=min(50,la);; for k from 1 to MAXLEN do # Can we raise pt? pttemp:=pt; for i from pt+1 to la do if a[i] <= heap then pttemp:=i; else break; fi; od; pt:=pttemp; # Extend b? if lb < MAXLEN then heap:=heap+a[pt]; b:=[op(b),heap]; lb:=lb+1; else RETURN(b); fi; od; RETURN(b); end: # SUPPORT positions where list is nonzero, starting count at 1 SUPPORT:=proc(a) local b,i; if whattype(a) <> list then RETURN([]); fi: b:=[]; for i from 1 to nops(a) do if a[i] <> 0 then b:=[op(b),i]; fi; od: RETURN(b); end: # MONO: sort abs values, discard duplicates MONO:=proc(a): if whattype(a) <> list then RETURN([]); fi: sort(convert(convert(map(abs,a),set),list)); end: # MONO2: sort abs values, discard dupes, return null unless seq. is changed MONO2:=proc(a) local b,bu: if whattype(a) <> list then RETURN([]); fi: b:=sort(convert(convert(map(abs,a),set),list)); if a = b then RETURN([]); fi; bu:=b[nops(b)]; if abs( bu - nops(b) ) < 10 then RETURN([]); else RETURN(b); fi; end: # COMP: compute complement of a sequence # i.e. vector [a[0]...a[i]...a[maxn]] which are numbers not in sequence # stopping at N if N < maxn is max |term| in sequence COMP:=proc(a) local maxn,notmember,b,c,n1,n,m: if whattype(a) <> list then RETURN([]); fi: maxn:=80: # get monotonic and unique version b b:=MONO(a); n1:=nops(b); if n1 = 0 then RETURN([]); fi: n:=b[n1]; m:=min(maxn,n); notmember:=subs(_WZ=b,proc(x) not member(x,_WZ) end); c:=select(notmember,[$1..m]); # don't return anything if complement is too long or too short #if nops(c) < 10 or nops(c) >50 then RETURN([]); else RETURN(c); fi; end: # COMPl: compute complement of a sequence (long version) # i.e. vector [a[0]...a[i]...a[maxn]] which are numbers not in sequence # USE WITH CARE!!! COMPl:=proc(a) local notmember,b,n1,m: if whattype(a) <> list then RETURN([]); fi: # get monotonic and unique version b b:=MONO(a); n1:=nops(b); if n1 = 0 then RETURN([]); fi: m:=min(b[n1],1000); notmember:=subs(_WZ=b,proc(x) not member(x,_WZ) end); select(notmember,[$1..m]); end: # CHAR: compute characteristic function of a sequence # i.e. vector [a[0]...a[i]...a[100]] which is 1 if |i| in sequence, 0 if not # stopping at N if N < 100 is max |term| in sequence CHAR:=proc(a) local maxn,ap,b,i,k,n,m: if whattype(a) <> list then RETURN([]); fi: maxn:=100: for n from 0 to maxn do b[n]:=0; od; m:=0: for n from 1 to nops(a) do ap:=abs(a[n]): if ap <= maxn then b[ap]:=1; fi; if ap > m then m:=ap; fi; od; [seq(b[n],n=0.. min(maxn,m) )]; end: # return subsequence of sequence[id] starting at /start/ with period /period/ EXTRACT:=proc(id,period,start) local i,a,b; a := nops(sequence[id])-start; b:=floor(a/period); [seq(sequence[id][period*i+start], i=0..b) ]; end: DECIMATE:=proc(a,k,j) local b,i,l: if whattype(a) <> list then RETURN([]); fi: l:=trunc( (nops(a)+k-1-j)/k ): if l < 1 then RETURN([]); fi: b:=[]: for i to l do b:=[op(b), a[k*(i-1)+j+1] ]: od: RETURN(b); end: BISECT:=proc(a,j) local b,i,l: if whattype(a) <> list then RETURN([]); fi: if j < 0 or j > 1 then RETURN([]); fi: l:=trunc( (nops(a)+1-j)/2 ): if l < 1 then RETURN([]); fi: b:=[]: for i to l do b:=[op(b), a[2*(i-1)+j+1] ]: od: RETURN(b); end: TRISECT:=proc(a,j) local b,i,l: if whattype(a) <> list then RETURN([]); fi: if j < 0 or j > 2 then RETURN([]); fi: l:=trunc( (nops(a)+2-j)/3 ): if l < 1 then RETURN([]); fi: b:=[]: for i to l do b:=[op(b), a[3*(i-1)+j+1] ]: od: RETURN(b); end: # Calculate first differences of a sequence DIFF:=proc(a) local b,i: if whattype(a) <> list then RETURN([]); fi: if nops(a) <= 1 then RETURN([]); fi: b:=[]: for i from 2 to nops(a) do b:=[op(b), a[i]-a[i-1] ]: od: RETURN(b); end: # Calculate all differences of a sequence. # Warning: unless the original sequence increases rapidly, # there could be small differences later in the sequence # that are missed! ALLDIFFS:=proc(a) local b,i,j: if whattype(a) <> list then RETURN([]); fi: if nops(a) <= 1 then RETURN([]); fi: b:={}: for i from 1 to nops(a)-1 do for j from i+1 to nops(a) do b:={op(b), a[j]-a[i] }: od: od: RETURN(sort(convert(b,list))); end: # PRODS Form sequence of partial products b[n] = a[1]*..*a[n] PRODS:=proc(a) local n,b,i,t1; if whattype(a) <> list then RETURN([]); fi: n:=nops(a); b:=[]; if n > 0 then t1:=1; for i from 1 to n do b:=[ op(b), t1 ]; t1:=t1*a[i]; od; fi; RETURN(b); end: PSUM:=proc(a) local b,i,s: if whattype(a) <> list then RETURN([]); fi: if nops(a) <= 0 then RETURN([]); fi: b:=[]: s:=0: for i from 1 to nops(a) do s:=s+a[i]: b:=[op(b), s ]: od: RETURN(b); end: PSUMSIGN:=proc(a) local b,i,s: if whattype(a) <> list then RETURN([]); fi: if nops(a) <= 0 then RETURN([]); fi: b:=[]: s:=0: for i from 1 to nops(a) do s:=a[i]-s: b:=[op(b), s ]: od: RETURN(b); end: BINOMIAL:=proc(a) local b,i,k: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i to nops(a) do b:=[op(b), add( binomial(i-1,k)*a[k+1], k=0..i-1)]: od: RETURN(b); end: BINOMIALi:=proc(a) local b,i,k: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i to nops(a) do b:=[op(b), add( (-1)^(i-1-k)*binomial(i-1,k)*a[k+1], k=0..i-1)]: od: RETURN(b); end: # BIN1 was introduced by Don Zagier (see M. Kaneko, # "A recurrence formula for the Bernoulli numbers", # Proc. Japan Acad., 71 A (1995), 192-193). It is an involution # on the class of sequences a = [a_0, a_1, a_2, ...], # sending a to b where b_n = (-1)^n Sum_{i=0..n} binomial(n+1,i+1) a_i. BIN1:=proc(a) local b,i,j,k: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i to nops(a) do j:=i-1; b:=[op(b), (-1)^j*add( binomial(j+1,k+1)*a[k+1] , k=0..j)]: od: RETURN(b); end: # EULER takes [a_1,a_2,a_3,...] and produces [b_1,b_2,b_3,...] with # 1 + Sum_{n=1..inf} b_n x^n = 1 / Product_{n=1..inf} (1-x^n)^a_n. EULER:=proc(a) local b,c,i,d: if whattype(a) <> list then RETURN([]); fi: c:=[]: for i to nops(a) do c:=[op(c), add( d*did(i,d)*a[d] , d=1..i)]: od: b:=[]: for i to nops(a) do b:=[op(b),(1/i)*(c[i]+ add( c[d]*b[i-d], d=1..i-1))]: od: RETURN(b); end: EULERi:=proc(b) local a,c,i,d: if whattype(b) <> list then RETURN([]); fi: c:=[]: for i to nops(b) do c:=[op(c),i*b[i]-add(c[d]*b[i-d], d=1..i-1)]: od: a:=[]: for i to nops(b) do a:=[op(a),(1/i)*add(mob(i,d)*c[d] , d=1..i)]: od: RETURN(a); end: MOBIUS:=proc(a) local b,i,d: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i to nops(a) do b:=[op(b), add( mob(i,d)*a[d], d=1..i)]: od: RETURN(b); end: MOBIUSi:=proc(a) local b,i,d: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i to nops(a) do b:=[op(b), add( did(i,d)*a[d] , d=1..i)]: od: RETURN(b); end: POLY:=proc(a) local t1,b,i,k: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i to nops(a) do b:=[op(b), add( (-1)^(i-1-k)*binomial(i-1,k)*a[k+1] , k=0..i-1)]: od: t1:= add ( b[i]*binomial(_n,i-1) , i=1..nops(a) ); sort(expand(expand(t1))): end: # POWERS extracts powers of x in f thru degree n POWERS:=proc(f,x,n) local i,g,lis; lis:=[]: g:=collect(expand(f),x); for i from 0 to n do if coeff(g,x,i) <> 0 then lis:=[op(lis),i]; fi; od; RETURN(lis); end: STIRLING:=proc(a) local b,i,k,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[]: for i from 1 to n do b:=[op(b), add( combinat[stirling2](i,k)*a[k], k=1..i)]: od: RETURN(b); end: STIRLINGi:=proc(a) local b,i,k,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[]: for i from 1 to n do b:=[op(b), add( combinat[stirling1](i,k)*a[k], k=1..i)]: od: RETURN(b); end: # The "Stirling-Bernoulli transform" maps a sequence b_0, b_1, b_2, ... # to a sequence c_0, c_1, c_2, ..., where if B has o.g.f. B(x), c has # e.g.f. exp(x)*B(1-exp(x)). More explicitly, # c_n = Sum_{m=0..n} (-1)^m*m!*Stirling2(n+1,m+1)*b_m. # Thanks to Masanobu Kaneko for telling me about this. STIRB:=proc(a) local b,j,k,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[]: for j from 0 to n-1 do b:=[op(b), add( (-1)^k*k!*combinat[stirling2](j+1,k+1)*a[k+1], k=0..j)]: od: RETURN(b); end: # Cameron's A and inverse: 1+SUM a_n x^n = 1/(1 - SUM b_n x^n); n=1..inf # i.e. 1 + A(x) = 1/(1 - B(x) ) # a's from b's: INVERT:=proc(a) local t1,t2,x,b,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a)+1: b:=listtoseries(a,x,'ogf'): t1:=series(1/(1-x*b),x,n): t2:=subsop(1=NULL, seriestolist(t1,'ogf')): if type(t2,list(integer)) then RETURN(t2) else RETURN([]): fi: end: # inverse, b's from a's: INVERTi:=proc(a) local t1,x,b,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a)+1: b:=listtoseries(a,x,'ogf'): t1:=series(-1/(1+x*b),x,n): RETURN(subsop(1=NULL, seriestolist(t1,'ogf'))): end: REVERT:=proc(a) local t1,t2,x,b,n: if whattype(a) <> list then RETURN([]); fi: if a[1] <> 1 then RETURN([]); fi: n:=nops(a)+1: b:=listtoseries(a,x,'ogf'): t1:=seriestoseries(b,'revogf'); t2:=subsop(1=NULL, seriestolist(t1,'ogf')): if type(t2,list(integer)) then RETURN(t2) else RETURN([]): fi: end: REVEGF:=proc(a) local t0,t3,i,t1,t2,x,b,n: if whattype(a) <> list then RETURN([]); fi: if a[1] <> 1 then RETURN([]); fi: n:=nops(a)+1: b:=listtoseries(a,x,'egf'): t0:=seriestoseries(b,'revogf'): t1:=seriestolist(t0,'ogf'); t3:=[seq(t1[i]*(i-1)!,i=1..nops(t1))]; t2:=subsop(1=NULL,t3); if type(t2,list(integer)) then RETURN(t2) else RETURN([]): fi: end: COMPOSE:=proc(a,b) local t1,n,f,g,T,U,x,i,j; if whattype(a) <> list then RETURN([]); fi: if whattype(b) <> list then RETURN([]); fi: n:=min( nops(a), nops(b) ): f:=unapply(convert( [seq( op(i,a)*T^i, i=1..n )], `+`), T); g:=unapply(convert( [seq( op(j,b)*U^j, j=1..n )], `+`), U); t1:=seriestolist(series( f(g(x)) ,x,n+1)); RETURN(subsop(1=NULL,t1)); end: COMPSQRT:=proc(b,M) # Compositional square root of b(x) = Sum b[i]*x^i, i=1..M # Finds a(x) = Sum b[i]*x^i, i=1..M, such that a(a(x)) = b(x). # Example: COMPSQRT(sin(x),10); should produce # x-1/12*x^3-1/160*x^5-53/40320*x^7-23/71680*x^9+... local i,j,a1,t1,t2,t3,a,e,s,A,B; # set up the equations B:=series(b,x,M+1); a1:=abs(sqrt(coeff(B,x,1))); A:=a1*x + add(a[i]*x^i, i=2..M); t1:=unapply(A,x); t2:=t1(t1(x)); t3:=series(t2-B,x,M+1); for i from 2 to M do e[i]:=coeff(t3,x,i); od: # solve them for i from 2 to M do s[i]:=solve(e[i],a[i]); for j from 2 to M do e[j]:=eval(subs(a[i]=s[i],e[j])); od: od: # compute answer series(a1*x + add(s[i]*x^i, i=2..M), x, M+1); end: CONV:=proc(a,b) local c,i,k,n: if whattype(a) <> list then RETURN([]); fi: if whattype(b) <> list then RETURN([]); fi: n:=min( nops(a), nops(b) ): c:=[]: for i from 0 to n-1 do c:=[ op(c), add( a[k+1]*b[i-k+1] , k=0..i)]: od; RETURN(c); end: CONVi:=proc(b) local n,B,a,A,i,t1,t2,t3: if whattype(b) <> list then RETURN([]); fi: if b[1]=0 then RETURN([]); fi: n:=nops(b): B:=listtoseries(b,x,'ogf'): a:=[]: for i from 0 to n-1 do a:=[op(a),a||i]: od: A:=listtoseries(a,x,'ogf'): t1:=series(A^2-B,x,n): for i from 0 to n-1 do t2:= solve( coeff( t1,x,i ), a||i ): if whattype(t2) = exprseq then t3:=t2[1]: else t3:=t2: fi: t1:=subs(a||i=t3, t1): A:=subs(a||i=t3, A): od: RETURN(seriestolist(A,'ogf')): end: LCMCONV:=proc(a) local b,i,k,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[]: for i from 0 to n-1 do b:=[op(b), add( lcm( a[k+1], a[i-k+1] ), k=0..i)]: od: RETURN(b); end: GCDCONV:=proc(a) local b,i,k,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[]: for i from 0 to n-1 do b:=[op(b), add( gcd( a[k+1], a[i-k+1] ), k=0..i)]: od: RETURN(b); end: ANDCONV:=proc(a) local b,i,k,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[]: for i from 0 to n-1 do b:=[op(b), add( ANDnos( a[k+1], a[i-k+1] ), k=0..i)]: od: RETURN(b); end: ORCONV:=proc(a) local b,i,k,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[]: for i from 0 to n-1 do b:=[op(b), add( ORnos( a[k+1], a[i-k+1] ), k=0..i)]: od: RETURN(b); end: XORCONV:=proc(a) local b,i,k,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[]: for i from 0 to n-1 do b:=[op(b), add( XORnos( a[k+1], a[i-k+1] ), k=0..i)]: od: RETURN(b); end: # form AND, OR, or XOR of two integers using their binary expansions ANDnos:=proc(n,m) local ans,n1,m1,k1,k2,q1,q2,sc; if n < 0 or m < 0 then ERROR(`negative argument`): fi: ans:=0: n1:=n: m1:=m: sc:=1: while n1 <> 0 or m1 <>0 do k1:=irem(n1,2,'q1'): k2:=irem(m1,2,'q2'): ans:=ans+sc*k1*k2; sc:=2*sc; n1:=q1: m1:=q2: od; RETURN(ans); end: ORnos:=proc(n,m) local k,ans,n1,m1,k1,k2,q1,q2,sc; if n < 0 or m < 0 then ERROR(`negative argument`): fi: ans:=0: n1:=n: m1:=m: sc:=1: while n1 <> 0 or m1 <>0 do k1:=irem(n1,2,'q1'): k2:=irem(m1,2,'q2'): if k1 = 1 or k2 = 1 then k:= 1; else k:= 0; fi; ans:=ans+sc*k; sc:=2*sc; n1:=q1: m1:=q2: od; RETURN(ans); end: XORnos:=proc(n,m) local k,ans,n1,m1,k1,k2,q1,q2,sc; if n < 0 or m < 0 then ERROR(`negative argument`): fi: ans:=0: n1:=n: m1:=m: sc:=1: while n1 <> 0 or m1 <>0 do k1:=irem(n1,2,'q1'): k2:=irem(m1,2,'q2'): if k1+k2 = 1 then k:= 1; else k:= 0; fi; ans:=ans+sc*k; sc:=2*sc; n1:=q1: m1:=q2: od; RETURN(ans); end: # Derivative of binary string of length n defined # by the binary expansion of u # Example: u=5 regarded as a binary string of length 4 is # 0101 with derivative (or difference) 111, so deriv(5,4) # returns 7. deriv:=proc(u,n) local t0,t1,t2; t0:=2^n+u; t1:=floor(t0/2); t2:=XORnos(t0,t1); ANDnos(t2,2^(n-1)-1); end: # LISTTOLISTDIV converts list to list by dividing nth term by n! LISTTOLISTDIV:=proc(a) local n,b,i; if whattype(a) <> list then RETURN([]); fi: n:=nops(a); b:=[]; for i from 1 to n do b:=[ op(b), a[i]/(i-1)! ]; od; end: # LISTTOLISTMULT converts list to list by multiplying nth term by n! LISTTOLISTMULT:=proc(a) local n,b,i; if whattype(a) <> list then RETURN([]); fi: n:=nops(a); b:=[]; for i from 1 to n do b:=[ op(b), a[i]*(i-1)! ]; od; end: # SERIESTOLISTDIV converts series to list by dividing nth term by n! SERIESTOLISTDIV:=proc(a) local t1,n,b,i; if whattype(a) <> series then RETURN([]); fi: t1:=seriestolist(a); n:=nops(t1); b:=[]; for i from 1 to n do b:=[ op(b), t1[i]/(i-1)! ]; od; end: # SERIESTOLISTMULT converts series to list by multiplying nth term by n! SERIESTOLISTMULT:=proc(a) local t1,n,b,i; if whattype(a) <> series then RETURN([]); fi: t1:=seriestolist(a); n:=nops(t1); b:=[]; for i from 1 to n do b:=[ op(b), t1[i]*(i-1)! ]; od; end: # SERIESTOSERIESDIV converts series to series by dividing nth term by n! SERIESTOSERIESDIV:=proc(a) local x,t1,n,b,i; if whattype(a) <> series then RETURN(0); fi: x:=op(0,a); t1:=seriestolist(a); n:=nops(t1); b:=0; for i from 1 to n do b:=b+t1[i]*x^(i-1)/(i-1)!; od; series(b,x,n+1); end: # SERIESTOSERIESMULT converts series to series by multiplying nth term by n! SERIESTOSERIESMULT:=proc(a) local x,t1,n,b,i; if whattype(a) <> series then RETURN(0); fi: x:=op(0,a); t1:=seriestolist(a); n:=nops(t1); b:=0; for i from 1 to n do b:=b+t1[i]*x^(i-1)*(i-1)!; od; series(b,x,n+1); end: # Procedure PARTITION takes a list l, # and determines the number P(i) of ways of # partitioning an integer i between 1 and max(list) into numbers from the list # (NOT counting repetitions). The output is the list [P(1), P(2),...)]. # l should start with 1, otherwise the 0's in the output will be skipped. # A second argument can be used to cut down the length of the output, which can # otherwise quickly grow unmanageable under iterations: the argument is a # number >=0, which specifies how much longer the output list should be than # the input list. PARTITION := proc() local f, g, i,translist,n,lp,l: l:=args[1]: if whattype(l) <> list then RETURN([]); fi: l:=convert(l,set); lp:=convert(l,list); if nargs=1 then n:=max(op(lp)) else n:=args[2] fi: f := 1: for i to nops(lp) while lp[i]<=n do f:=f/(1-t^op(i,lp)) od: g:=taylor(f, t=0, n+1): translist := []: for i from 2 to nops(g)/2-1 do translist:=[op(translist),op(2*i-1,g)]: od: translist; end: PARTITIONi:=proc(l) local t,i,j,l1: if whattype(l) <> list then RETURN([]); fi: t:=EULERi(l); l1:=[]: for i to nops(t) do for j from 1 to t[i] do l1:=[op(l1),i] od: od: l1; end: WEIGH := proc() local a,x,f,n,g,i: a:=args[1]: if whattype(a) <> list then RETURN([]); fi: if nargs=1 then n:=max(op(a)) else n:=args[2] fi: f := 1: for i to nops(a) do f:=f*(1+x^a[i]) od: g:=taylor(f, x=0, n+1): RETURN( subsop( 1=NULL, seriestolist(g,'ogf'))): # RETURN( seriestolist(g)): end: # EXP converts [a_1, a_2, ...] to [b_1, b_2,...] where # 1 + EGF_B (x) = exp EGF_A (x) (Ref. Eigen Seq paper page 61) EXP:=proc(a) local i,t0,t1,t2,x,b,n: if whattype(a) <> list then RETURN([]); fi: b:=[0,op(a)]; n:=nops(b)+1: t0:=series(exp(listtoseries(b,x,'egf')),x,n): t1:=seriestolist(t0,'ogf'); t2:=[seq(t1[i]*(i-1)!,i=1..nops(t1))]; RETURN(subsop(1=NULL,t2)); end: EXPCONV:=proc(a,b) local i,t0,t1,x,n: if whattype(a) <> list then RETURN([]); fi: if whattype(b) <> list then RETURN([]); fi: n:=min(nops(a),nops(b)): t0:=series(listtoseries(a,x,'egf')*listtoseries(b,x,'egf'),x,n); t1:=seriestolist(t0,'ogf'); [seq(t1[i]*(i-1)!,i=1..nops(t1))]; end: # LOG converts [a_1, a_2, ...] to [b_1, b_2,...] where # 1 + EGF_A (x) = exp EGF_B (x) (Ref. Eigen Seq paper page 61) # i.e. EGF_A (x) = log( 1 + EGF_B (x) ). LOG:=proc(a) local t2,t0,i,t1,x,b,n: if whattype(a) <> list then RETURN([]); fi: b:=[1,op(a)]; n:=nops(b)+1: t0:=series(log(listtoseries(b,x,'egf')),x,n): t1:=seriestolist(t0,'ogf'); t2:=[seq(t1[i]*(i-1)!,i=1..nops(t1))]; RETURN(subsop(1=NULL,t2)); end: # shift left LEFT:=proc(a) RETURN(subsop(1=NULL, a)): end: # shift right RIGHT:=proc(a) RETURN([1,op(a)]): end: M2:=proc(a) local i,b,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[a[1]]: if n=1 then RETURN(b); fi: for i from 2 to n do b:=[op(b), 2*a[i] ]; od: RETURN(b); end: M2i:=proc(a) local i,b,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[a[1]]: if n=1 then RETURN(b); fi: for i from 2 to n do b:=[op(b), a[i]/2 ]; od: RETURN(b); end: NEGATE:=proc(a) local i,b,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[a[1]]: if n=1 then RETURN(b); fi: for i from 2 to n do b:=[op(b), -a[i] ]; od: RETURN(b); end: # other transforms: # weigh2out: b's from a's in # 1+SIGMA b_n x^n = nonnegative part of PI (x^-n + 1 + x^n)^a_n # a_n ---> b_n weigh2out := proc(a) local t0,n,x,f,i: if whattype(a) <> list then RETURN([]); fi: n:=nops(a)+1: f:= 1: for i to nops(a) do f:=f*(x^(-i)+1+x^i)^a[i]: od: t0:=seriestoseries( series( expand(f),x,n ),'ogf'): RETURN( seriestolist(t0,'ogf')): end: # eultrans2 eultrans2:=proc(a) local t0,x,n: if whattype(a) <> list then RETURN([]); fi: n:=nops(a)+1: t0:=series( (1/(1+x))*subs( x=x/(1+x), listtoseries(a,x,'ogf') ) , x, n): RETURN( seriestolist(t0,'ogf')): end: # weighout: b's from a's in # 1+SIGMA b_n x^n = PI (1+x^n)^a_n # (was Etrans2): a_n ---> b_n weighout := proc(l) local x,f, g, i : if whattype(l) <> list then RETURN([]); fi: f := 1: for i to nops(l) do f:=f*(1+x^i)^l[i] od: g:=taylor(f, x=0, nops(l)+1): RETURN( subsop( 1=NULL, seriestolist(g,'ogf'))): end: # weighouti: a's from b's in # 1+SIGMA b_n x^n = PI (1+x^n)^a_n # b_n ---> n_n weighouti := proc(b) local c,i,d,a: if whattype(b) <> list then RETURN([]); fi: c:=[]: for i to nops(b) do c:=[op(c),i*b[i]-add(c[d]*b[i-d], d=1..i-1)]: od: a:=[c[1]]; for i from 2 to nops(b) do a:=[op(a), c[i]/i- (1/i)*add(-dids(i,d)*a[d]*d , d=1..i-1)]: od: RETURN(a): end: # EULER (again) # Procedure Etrans takes a list l and determines the Euler transform # of the list. The output list starts with the *1st* coefficient in # the Taylor series. Etrans := proc(l) local f, g, i, j, translist: if whattype(l) <> list then RETURN([]); fi: f := 1: for i to nops(l) do f:=f/((1-t^i)^op(i,l)) od: g:=taylor(f, t=0, nops(l)+1): translist := []: for i from 2 to nops(g)/2-1 do for j from op(2*i-2,g)+1 to op(2*i,g)-1 do translist:=[op(translist),0]: od: translist:=[op(translist),op(2*i-1,g)]: od: translist; end: # Procedure ptrans takes a list l, # and determines the number P(i) of ways of partitioning an integer i between 1 # and max(list) into numbers from the list (WITH REPETITIONS). The output is # the list [P(1), P(2),...,P(m)]. l should start with 1, otherwise the 0's # in the output will be skipped. # A second argument can be used to cut down the length of the output, which can # otherwise quickly grow unmanageable under iterations: the argument is a # number >=0, which specifies how much longer the output list should be than # the input list. ptrans := proc() local f, g, i,translist,n,l: l:=args[1]: if whattype(l) <> list then RETURN([]); fi: if nargs=1 then n:=max(op(l)) else n:=args[2]+nops(l) fi: f := 1: for i to nops(l) do f:=f/(1-t^op(i,l)) od: g:=taylor(f, t=0, n+1): translist := []: for i from 2 to nops(g)/2-1 do translist:=[op(translist),op(2*i-1,g)]: od: translist; end: # pairtrans pairtrans:=proc(a) local i,b,n; if whattype(a) <> list then RETURN([]); fi: n:=nops(a): b:=[a[1]]: if n=1 then RETURN(b); fi: for i from 2 to n do b:=[op(b), a[i]+a[i-1] ]; od: RETURN(b); end: # I clarified these definitions: # The BOUS transform takes [a1,a2,a3,...] to [b0=1,b1,b2,...], # where the triangle is # b0=1 # a1 ---> b1 # b2 * <--- a2 # a3 ---> * * b3 # etc. # # The BOUS2 transform takes [a0,a1,a2,a3,...] to [b0,b1,b2,...], # where the triangle is # a0=b0 # a1 ---> b1 # b2 * <--- a2 # a3 ---> * * b3 # etc. BOUS:=proc(a) local c,i,j,n: if whattype(a) <> list then RETURN([]); fi: n:= min( nops(a), 60); c[0,0]:=1; for i to n do c[i,0]:=a[i]; od; for i to n do for j to i do c[i,j] := c[i,j-1] + c[i-1,i-j]; od; od; RETURN([seq(c[i,i],i=0..n)]); end: BOUS2:=proc(a) local c,i,j,n: if whattype(a) <> list then RETURN([]); fi: n:= min( nops(a), 60); for i from 0 to n-1 do c[i,0]:=a[i+1]; od; for i to n-1 do for j to i do c[i,j] := c[i,j-1] + c[i-1,i-j]; od; od; RETURN([seq(c[i,i],i=0..n-1)]); end: BOUS2i:=proc(b) local c,i,j,n: if whattype(b) <> list then RETURN([]); fi: n:= min( nops(b), 60); for i from 0 to n-1 do c[i,0]:=b[i+1]; od; for i to n-1 do for j to i do c[i,j] := c[i,j-1] - c[i-1,i-j]; od; od; RETURN([seq(c[i,i],i=0..n-1)]); end: CONTINUANT:=proc(a) option remember; local i,n,t0; # The n-th term of the transformed sequence is K_n(a_1, a_2, ..., a_n) # - cf. Graham, Knuth and Patashnik, Concrete Math., 2nd. ed., p. 302. # See A072347 for definition of continuant. if whattype(a) <> list then RETURN([]); fi: n:=nops(a); if n = 0 then RETURN([]); fi; if n = 1 then RETURN(a); fi; t0:=[1,a[1]]; for i from 2 to n do t0:=[op(t0), a[i]*t0[i]+t0[i-1] ]; od; expand(subsop(1=NULL,t0)); end: CONTINUANTi:=proc(a) option remember; local i,n,t0,t1; # Inverse CONTINUANT transform # Warning: does not in general produce integers if whattype(a) <> list then RETURN([]); fi: n:=nops(a); if n = 0 then RETURN([]); fi; if n = 1 then RETURN(a); fi; t0:=[a[1]]; t1:=[1,op(a)]; for i from 2 to n do t0:=[op(t0), (t1[i+1]-t1[i-1])/t1[i] ]; od; expand(t0); end: digsum:=proc(n) local i,m,n1; m:=0; n1:=abs(n); while n1<>0 do i:= n1 mod 10; m:=m+i; n1:=(n1-i)/10; od; RETURN(m); end: digrev:=proc(n) local b,i,j,n1,m; m:=0; b:=[]; n1:=abs(n); while n1<>0 do i:= n1 mod 10; b:=[op(b), i]; n1:=(n1-i)/10; od; for j from 1 to nops(b) do m:=10*m+b[j]; od; RETURN(m); end: digsort:=proc(n) local b,i,j,n1,m; m:=0; b:=[]; n1:=abs(n); while n1<>0 do i:= n1 mod 10; b:=[op(b), i]; n1:=(n1-i)/10; od; b:=sort(b); for j from 1 to nops(b) do m:=10*m+b[j]; od; RETURN(m); end: DIGSUM:=proc(a) local b,i: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i to nops(a) do b:=[op(b), digsum(a[i])]: od: RETURN(b); end: DIRICHLET:=proc(a,b) local c,i,d,s: # Dirichlet convoltion of [a_1, ...] and [b_1, ...] # Ref. Apostol, Intro. to Analytic Number Theory, Section 11.4, page 228. if whattype(a) <> list then RETURN([]); fi: if whattype(b) <> list then RETURN([]); fi: c:=[]: for i to min(nops(a),nops(b)) do s:=0; for d from 1 to i do if did(i,d) = 1 then s:=s+a[d]*b[i/d]; fi; od: c:=[op(c),s]: od: RETURN(c); end: DIGREV:=proc(a) local b,i: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i to nops(a) do b:=[op(b), digrev(a[i])]: od: RETURN(b); end: # mex = minimal excluded member of a set mex := proc(s) local n; for n from 0 while member(n,s) do od; n end: # nimsum of two numbers # write as sums of powers of 2, add them cancelling equal pairs nimsum:=proc(n1,n2) local v1,v2,l1,l2,i,ans,t1,t2; ans:=0; v1:=convert(n1,base,2); v2:=convert(n2,base,2); l1:=nops(v1); l2:=nops(v2); if l1 > l2 then t1:=v1; v1:=v2; v2:=t1; t2:=l1; l1:=l2; l2:=t2; fi; # now l1 is the shorter if l1 > 0 then for i from 1 to l1 do if v1[i] <> v2[i] then ans:=ans+2^(i-1); fi; od; fi; if l2 > l1 then for i from l1+1 to l2 do if v2[i] <> 0 then ans:=ans+2^(i-1); fi; od; fi; ans; end: # Compute weight or number of 1's in binary expansion of n: wt:=proc(n) local w,m,i; w:=0; m:=n; while m > 0 do i := m mod 2; w:=w+i; m:=(m-i)/2; od; w; end: # little Hankel lHANKEL:=proc(a) local b,i,k: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i from 2 to nops(a)-1 do b:=[op(b), a[i]^2 - a[i+1]*a[i-1] ]; od: RETURN(b); end: # Hankel HANKEL:=proc(a) local h,M,b,i,j,n,m0: if whattype(a) <> list then RETURN([]); fi: m0:=floor((nops(a))/2); b:=[]: h:=(i,j)->a[i+j-1]; for n from 1 to m0 do M:=linalg[matrix](n,n,h); b:=[op(b),det(M)]; od: RETURN(b); end: # Shadow - see Lorenz Halbeisen and Norbert Hungerbuehler reference in sequence database under A000522 # Given a sequence a(0), a(1), ... the shadow is the sequence s(0), s(1), ... # where s(n) = number of i with 0 <= i < n such that n divides a(i). SHADOW:=proc(a) local n,b,t,j,t1: if whattype(a) <> list then RETURN([]); fi: n:=nops(a); if n = 0 then RETURN([]); fi: if n = 1 then RETURN([0]); fi: b:=[0]: for t from 1 to n-1 do t1:=0; for j from 0 to t-1 do if a[j+1] mod t = 0 then t1:=t1+1; fi; od: b:=[op(b), t1]: od: RETURN(b); end: # SERIES2(f,x,y,n) returns the sum of all monomials x^i*y^j # in f with i < n and j < n. SERIES2:=proc(f,x,y,n) local i,j,t0,t1,t2,t3,t4; t0:=0; t1:=series(f,x,n); for i from 0 to n-1 do t2:=coeff(t1,x,i); t3:=series(t2,y,n); for j from 0 to n-1 do t4:=coeff(t3,y,j); if i+j < n then t0:=t0+t4*x^i*y^j; fi: od; od; sort(t0,[x,y], tdeg); end: # SERIES2TOLIST(f,x,y,n) takes output of SERIES2(f,x,y,n) # and produces a list obtained by reading the monomials in f # in the order 00, 10, 01, 20, 11, 02, 30, 21, 12, 03, ... # where total degree is < n SERIES2TOLIST:=proc(f,x,y,n) local i,j,k,t0,t1,t2; t0:=[]; for k from 0 to n-1 do for j from 0 to k do i:=k-j; t1:=coeff(f,x,i); t2:=coeff(t1,y,j); t0:=[op(t0),t2]; od; od; t0; end: # SERIES2TOLISTMULT(f,x,y,n) takes output of SERIES2(f,x,y,n) # and produces a list obtained by reading the monomials in f # in the order 00, 10, 01, 20, 11, 02, 30, 21, 12, 03, ... # where total degree is < n # and multipying the coefficient of x^i*y^j by i!*j!. SERIES2TOLISTMULT:=proc(f,x,y,n) local i,j,k,t0,t1,t2; t0:=[]; for k from 0 to n-1 do for j from 0 to k do i:=k-j; t1:=coeff(f,x,i); t2:=i!*j!*coeff(t1,y,j); t0:=[op(t0),t2]; od; od; t0; end: # SERIES2HTOLIST(f,x,y,n) takes output of SERIES2(f,x,y,n) # and produces a list obtained by reading the monomials in f # in the order 00, 10, 11, 20, 21, 22, 30, 31, 32, 33, ... # where both exponents are < n # It will complain if this would omit any terms - if there is a term x^i*y^j # with j>i for i 0 then lprint("Error in SERIES2HTOLIST: missed term in row", i, t2*x^i*y^j); error "Perhaps try SERIES2TOLIST?"; fi; od: od; t0; end: # SERIES2HTOLISTMULT(f,x,y,n) takes output of SERIES2(f,x,y,n) # and produces a list obtained by reading the monomials in f # in the order 00, 10, 11, 20, 21, 22, 30, 31, 32, 33, ... # and multipying the coefficient of x^i*y^j by i!*j! # where both exponents are < n # It will complain if this would omit any terms - if there is a term x^i*y^j # with j>i for i 0 then lprint("Error in SERIES2HTOLISTMULT: missed term in row", i, t2*x^i*y^j); error "Perhaps try SERIES2TOLISTMULT?"; fi; od: od; t0; end: # Transforms ADDONE, RIGHT0, SUMTABL, RAST and RASTxx contributed by # Antti Karttunen 12. July 2002. add1 := n -> n+1: ADDONE := a -> map(add1,a): # Increment each term by one. RIGHT0 := a -> [0,op(a)]: # Shift right, prepending 0. # SUMTABL(A001477,A001477) gives A003056. # SUMTABL(A000027,A000027) gives A003057. # and SUMTABL(A000027,A001477) & SUMTABL(A001477,A000027) give A002024. A025581 := n -> binomial(1+floor((1/2)+sqrt(2*(1+n))),2) - (n+1): # The X-projection & A002262 := n -> n - binomial(floor((1/2)+sqrt(2*(1+n))),2): # Y-projection (of A001477) SUMTABL := proc(a,b) local c,i,u; u := binomial(min(nops(a),nops(b))+1,2); c := []; for i from 0 to u-1 do c := [ op(c), a[A025581(i)+1]+b[A002262(i)+1] ]; od; RETURN(c); end: # Theorem: The set of the sequences where each n >= 0 occurs A000108(n) times # (i.e. the permutations of A072643) is closed under these two transformations. # The fixed point of RASTxx is A071673. RAST := (a,b) -> RIGHT0(ADDONE(SUMTABL(a,b))): RASTxx := a -> RIGHT0(ADDONE(SUMTABL(a,a))): # End of contributions from Antti Karttunen 12. July 2002. # RECORDS: Extract records and their positions from a list # Returns pair b, c where b = records, c = where they occurred (starting at 1) # Example: with(numtheory): t1:=[seq(sigma(n),n=1..200)]; #gives A000203; RECORDS(t1); # gives [A034885, A002093]; RECORDS := proc(a) local i,n,b,c,rec; if whattype(a) <> list then RETURN([[],[]]); fi: n:=nops(a); if n=0 then RETURN([[],[]]); fi: if n=1 then RETURN([[a[1]],[1]]); fi; rec:=a[1]; b:=[rec]; c:=[1]; for i from 2 to n do if a[i] > rec then rec:=a[i]; b:=[op(b),rec]; c:=[op(c),i]; fi; od: RETURN([b,c]); end: # LAH and LAHi were suggested by Vladeta Jovovic (vladeta(AT)Eunet.yu), Oct 03, 2003 # The "LAH transform" maps a sequence a_0, a_1, a_2, ... # to a sequence b_0, b_1, b_2, ..., where if A has e.g.f. A(x), b has # e.g.f. A(x/(1-x)). More explicitly, # b_n = Sum_{m=0..n} (n!/m!)*binomial(n-1,m-1)*a_m. # Examples are sequences A084357 (from A000110) and A084358 (from A000670). LAH:=proc(a) local b,i,m: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i to nops(a) do b:=[op(b), add( ((i-1)!/m!)*binomial(i-2,m-1)*a[m+1] , m=0..i-1)]: od: RETURN(b); end: # LAH and LAHi were suggested by Vladeta Jovovic (vladeta(AT)Eunet.yu), Oct 03, 2003 # The "LAHi transform" maps a sequence a_0, a_1, a_2, ... # to a sequence b_0, b_1, b_2, ..., where if A has e.g.f. A(x), b has # e.g.f. A(x/(1+x)). More explicitly, # b_n = Sum_{m=0..n} (-1)^(n-m)*(n!/m!)*binomial(n-1,m-1)*a_m. # Examples are sequences A000110 (from A084357) and A000670 (from A084358). LAHi:=proc(a) local b,i,m: if whattype(a) <> list then RETURN([]); fi: b:=[]: for i to nops(a) do b:=[op(b), add( (-1)^(i-1-m)*((i-1)!/m!)*binomial(i-2,m-1)*a[m+1] , m=0..i-1)]: od: RETURN(b); end: