TRANSFORMATIONS OF INTEGER SEQUENCES
N. J. A. Sloane
This is a 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://oeis.org)
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 ####
# This is the official version of transforms
# May 10 2020: Added MP, COMPl2, binary
# May 09 2020: corrected POLY, WEIGH
# May 26 2017: added new operations, fixed REVERT etc
# Oct 12 2013: added digprod, digprod0, improved digsum
# Jan 05 2013: Corrections and test code (at end of file) from Michael Somos
# Nov 26 2011: Added FREQ
# Jul 16 2011: Added LAGRANGE
# Oct 22 2010: added two Van Eck transforms
# Oct 16 2010: Richard Mathar sent AERATE, DIRICHLETi, LODUMO, digcat2, digcatL
# added RUNS Nov 29 2009
# added stoint Jun 27 2009
# added ConvOffs and ConvOffsStoT from Richard Mathar, Feb 01 2008
# added INVERSE Oct 20 2005
# added INVERSE2 Feb 16 2017
# 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
# Memo to myself: to create the file transforms.txt do
# make transforms.txt
# 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://neilsloane.com/doc/eigen.pdf
# Summary:
# AERATE aerate by insertion of zeros into a sequence
# 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)
# COMPl2 compute complement of a sequence (long version 2)
# 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
# ConvOffs see code for description
# ConvOffsStoT see code for description
# 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
# DIRICHLETi Dirichlet inverse of an arithmetic sequence
# ECKb Backwards Van Eck transform
# ECKf Forwards Van Eck transform
# 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
# FREQ b[n] = number of i <= n such that a[i]=a[n]
# GCDCONV GCD-convolution, SUM a(k).gcd.a(n-k)
# 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
# INVERSE inverse of permutation
# INVERSE2 smallest inverse of n in a, or 0 if missing value
# 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
# LAGRANGE How many terms of a are needed to sum to n
# 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!
# LODUMO Deleham's Lodumo_m transform
# 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 different from orig.
# MP Maple Print (prints a sequence on one long line)
# myhead first L (default 100) terms of a sequence
# mytail last L (default 100) terms of a sequence
# 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 sequence from leading diagonal of 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
# RECORDS extracts records from a list
# REVERT reversion
# REVEGF reversion of e.g.f.
# RIGHT shift right, prepending a 1
# RUNS lengths of successive runs of identical terms
# SEQPI produces the pi function for a monotonic sequence a[]
# 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
# UNBASE Does opposite of convert(n,base,b)
# WEIGH b from a in 1+SUM b_n x^n=PROD(1+x^a_n)
# XORCONV XOR-convolution, SUM a(k).XOR.a(n-k)
# others:
# ANDnos logical AND of two numbers using their binary expansions
# binary binary representation of n, in human order
# delta delta(f,n,k) = kth term of nth difference of f
# deriv derivative of a numbers using its binary expansion
# did did it divide?
# dids did it divide (with sign)?
# digcat2 concatenate two numbers
# digcatL concatenate a list of numbers
# digrev reverse digits
# digprod product of digits
# digprod0 product of nonzero digits
# digsort sort digits into increasing order
# digsortrev sort digits into decreasing order
# digsum sum of digits
# dim did it mask?
# Etrans Euler Xfm (again)
# eultrans2 2nd Euler trans from paper by Donaghey and Shapiro
# but this is same as BINOMIALi
# mex minimum excluded number
# myconvert myconvert(n,b) = convert(n,base,b) in right order
# nimsum nim sum
# ORnos logical OR of two numbers using their binary expansions
# pairtrans b(n)=a(n)+a(n-1)
# pairtransi a(n)=b(n)-b(n-1)+b(n-2)-...
# ptrans partition Xfm (with repetition)
# stoint convert string to integer
# 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)
# wt Compute weight or number of 1's in binary expansion of n:
# XORnos logical XOR of two numbers using their binary expansions
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 ([email protected]) message
# "half-baked generalized convolution sequence transforms"
# posted 09 Jun 2001 on SeqFan mailing list ([email protected])
# 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([]); 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 S starting at /start/ with period /period/
EXTRACT:=proc(S,period,start) local i,a,b;
if whattype(S) <> list then RETURN([]); fi:
a := nops(S)-start; b:=floor(a/period);
[seq(S[period*i+start], i=0..b) ]; end:
# b[n] = number of i <= n such that a[i]=a[n].
# Cf. A200779, A200780
FREQ:=proc(a) local b,i,c,n,t1:
if whattype(a) <> list then RETURN([]); fi:
b:=[];
for n from 1 to nops(a) do
c:=0; t1:=a[n];
for i from 1 to n do if a[i]=t1 then c:=c+1; fi; od:
b:=[op(b),c];
od;
RETURN(b);
end:
# return subsequence of terms == j mod k
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
t1:=t1*a[i];
b:=[ op(b), t1 ];
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 get polynomial for sequence from leading diagonal of difference table
# Corrected by njas May 09 2020
POLY:=proc(a) local t1,b,i,k:
if whattype(a) <> list then RETURN([]); fi:
sort(expand(add(binomial(_n,k)*a[k+1], k=0..nops(a)-1)));
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 Reversion of a sequence.
# Corrected (so that it works with the latest version of Maple), May 26 2017
# To test this, set a:=[1,3,5,7,9,11,13,15,17]; then REVERT(a) should produce
# [1, -3, 13, -67, 381, -2307, 14589, -95235, 636925]
REVERT:=proc(a) local ashift,t1,t2,x,b,n:
if whattype(a) <> list then RETURN([]); fi:
if a[1] <> 1 then RETURN([]); fi:
ashift:=[0,op(a)];
b:=listtoseries(ashift,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 Exponential reversion of a sequence.
# Corrected (so that it works with the latest version of Maple), May 26 2017
# To test this, set a:=[1,3,5,7,9,11,13,15,17]; then REVEGF(a) should produce
# [1, -3, 22, -262, 4336, -91984, 2381408, ...]
REVEGF:=proc(a) local ashift,t0,t3,i,t1,t2,x,b,n:
if whattype(a) <> list then RETURN([]); fi:
if a[1] <> 1 then RETURN([]); fi:
ashift:=[0,op(a)];
b:=listtoseries(ashift,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 transform: b from a in 1+SUM b_n x^n = PROD(1+x^a_n)
# Usage: WEIGH([1,2,4,8]); will produce [1, 1, 1, 1, 1, 1, 1, 1]
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:
# Calculate sequence of lengths of runs of identical terms
RUNS:=proc(a) local b,i,prev,ct:
if whattype(a) <> list then RETURN([]); fi:
if nops(a) = 1 then RETURN([1]); fi;
b:=[]:
for i from 1 to nops(a) do
if i=1 then prev:=a[1]; ct:=1;
else
if a[i] = prev then ct:=ct+1;
else b := [op(b), ct]; ct:=1; prev:=a[i];
fi;
fi;
if i = nops(a) then b:=[op(b),ct]; fi;
od;
RETURN(b);
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) add(d, d=convert(n, base, 10)) ; end proc: # A007953
digprod := proc(n) mul(d, d=convert(n, base, 10)) ; end proc: # A007954
digprod0 := proc(n) local d, j: d:=convert(n, base, 10): return mul(`if`(d[j]=0, 1, d[j]), j=1..nops(d)): end: # A051801
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:
digsortrev:=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);
b:=ListTools[Reverse](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),linalg[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