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 (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([]); 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<n.
# The H stands for "horizontal".
SERIES2HTOLIST:=proc(f,x,y,n) local i,j,k,t0,t1,t2;
t0:=[]; 
for i from 0 to n-1 do
t1:=coeff(f,x,i);
for j from 0 to i do
t2:=coeff(t1,y,j);
t0:=[op(t0),t2];
od;
# check for possible lost terms
for j from i+1 to n-1 do
t2:=coeff(t1,y,j);
if t2 <> 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<n.
# The H stands for "horizontal".
SERIES2HTOLISTMULT:=proc(f,x,y,n) local i,j,k,t0,t1,t2;
t0:=[]; 
for i from 0 to n-1 do
t1:=i!*coeff(f,x,i);
for j from 0 to i do
t2:=j!*coeff(t1,y,j);
t0:=[op(t0),t2];
od;
# check for possible lost terms
for j from i+1 to n-1 do
t2:=j!*coeff(t1,y,j);
if t2 <> 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:

# SEQPI: produces the pi function for a monotonic sequence a[]
# Returns list b[], where b[n+1} gives number of terms in range 0 to n.
# Generates at least M terms of b (if possible) and then stops at the next convenient stopping point.
SEQPI:=proc(a,M)  local b,la,L,i,s,at,ct:
if whattype(a) <> list then RETURN([]); fi:
la:=nops(a); b:=[]: L:=min(a[la],M);
ct:=0;
   s:=0;
for at from 1 to la do
   for i from s to a[at]-1 do b:=[op(b),ct]; od:
   ct:=ct+1;
   s:=a[at];
   if s>L then break; fi;
                    od:
b:=[op(b),ct];
RETURN(b);
end:


# UNBASE: Does opposite of convert(n,base,b)
UNBASE:=proc(a)  local N,i;
if whattype(a) <> list then RETURN([]); fi:
N:=add(a[i]*10^(i-1),i=1..nops(a));
RETURN(N);
end:


INVERSE := proc(a) local b,c,i,n,hit,ulim;
if whattype(a) <> list then RETURN([]); fi:
c:=sort(a);
b:=[];
ulim:=min(1000,c[nops(c)]);
for n from c[1] to ulim do
  hit:=0;
  for i from 1 to nops(a) do
  if a[i]=n then b:=[op(b),i]; hit:=1; break; fi;
  od;
if hit = 0 then RETURN(b); fi;
od;
RETURN(b);
end:

# INVERSE2 smallest inverse of n in a, or 0 if missing value
# ff = offset of a
INVERSE2 := proc(a,ff) local nov,b,c,i,n,llim,ulim;
if whattype(a) <> list then RETURN([]); fi:
c:=sort(a);
nov:=nops(convert(c,set));
lprint("Contains",nov,"values, in the range",c[1],"to",c[-1]);
llim:=max(-1000,c[1]);
ulim:=min(1000,c[-1]);
b:=Array(llim..ulim,0);
for n from llim to ulim do
  if member(n,a) then 
    for i from 1 to nops(a) do if a[i]=n then b[n]:=i+ff-1; break; fi; od:
  fi;
                     od;
RETURN([seq(b[n],n=llim..ulim)]);
end:

# ConvOffs: from  Richard J. Mathar, Jan 10 2008
# input:
#   L a list
# output:
#   a list C with one more element than L. C[1]=L[1].
#   C[2]=L[-1]. Recursively C[i]=C[i-1]*L[1-i]/L[i-1], 3<=i.
#   If L is empty, the list [1] is returned.
#
# Example: L=[1,4,10,20] returns [1,20,50,20,1]
# Example: L=[1,2,3] returns [1,3,3,1]
# Example: L=[1] returns [1,1]
# Example: L=[] returns [1]
#
ConvOffs := proc(L)
	local C,i ;
	if nops(L) > 0 then
		C := [op(1,L),op(-1,L)] ;
		for i from 2 to nops(L) do
			C := [op(C), op(-i,L)*op(i,C)/op(i,L) ] ;
		od:
	else
		C := [1] ;
	fi ;
	RETURN(C) ;
end:

# ConvOffsStoT: from  Richard J. Mathar, Jan 10 2008
# input:
#    S a list (that is, the first members of a virtually infinite sequence)
#
#    verb verbosity parameter: if set to true, the output is printed
#     as a stream of numbers in a row-by-row fashion separated by commas.
#
# output:
#    a list of lists, where each of the lists represents a row n of
#    a triangle, n=0 up to the length of S, and is the list defined by
#    ConvOffs(sublist of first n elements of S).
#    
# Example A056939: L=[1,4,10,20] returns
#   [[1], [1, 1], [1, 4, 1], [1, 10, 10, 1], [1, 20, 50, 20, 1]]
#   where [1] is generated by ConvOffs([]), [1,1] is generated by
#   ConvOffs([1]), [1,4,1] is generated by ConvOffs([1,4]),
#   [1,10,10,1] is generated by ConvOffs([1,4,10]) etc.
#    
# tests:
# ConvOffs([1,2,3]) ;
# ConvOffsStoT([1,4,10,20],true) ;
#
ConvOffsStoT := proc(S,verb)
	local T,n, nrow,c ;
	# first row is from empty sub-list
	T := [[1]] ;
	for n from 1 to nops(S) do
		nrow := ConvOffs([op(1..n,S)]) ;
		if verb then
			for c in nrow do
				printf("%d, ",c) ;
			od:
		fi ;
		T := [op(T),nrow] ; 
	od:
	RETURN(T) ;
end:


# stoint: convert string to integer
stoint := proc(s)
local n,a,a0,i;
n:=length(s);
a:=convert(s,'bytes');
a0:=0:
for i from 1 to n do
       a0:=(a[i]-48)*10^(n-i)+a0:
end do:
a0;
end:


# Dirichlet inverse of an arithmetic function.
# The input list L represents [a1,a2,a3,a4,....], the values of the original
#   arithmetic function at integer values of 1, 2, 3 etc. The value of a1 must
#   be nonzero.
# The list returned is the inverse of the Dirichlet convolution, [b1, b2, b3,...].
# This means that the DIRICHLET(.,.) function returns the
# identity if called with two parameters equal to the input and output of
# DIRICHLETi(.).
# 
# Example:
# L := [seq(n mod 2, n=1..50)] ; # this is A000035
# DIRICHLETi(L) ; # this is A087003
# Richard J. Mathar, 2010-07-15
DIRICHLETi := proc(L)
	local Linv,n,d,finv;
	if not type(L,'list') then
		ERROR("Expected list type argument instead of ",whattype(L) ) ;
	end if;
	Linv := [1/op(1,L)] ;
	for n from 2 to nops(L) do
		finv := 0 ;
		for d in ( numtheory[divisors](n) minus {n} ) do
			op(n/d,L)*op(d,Linv) ;
			finv := finv - % * op(1,Linv) ;
		end do:
		Linv := [op(Linv),finv] ;
	end do:
	return Linv ;
end proc:

# Deleham's Lodumo_m transform of a sequence L.
# The list a(n) is returned and defined by the smallest number not yet in
# the list a(.) satisfying a(n)=L(n) (mod m).
# @param L the input list of non-negative numbers 
# @param m  the remainder of the modulo operation. This is  only useful
#           if it's a number larger than 1.
# @return the sequence of nonnegative numbers a(n)
#
# Example: generate A160016 from A159833:
#     A159833 := [seq(n^2*(n^2+15)/4,n=0..100) ];
#     LODUMO(A159833,2) ;
#
# Example: generate A160081 from A000045:
#     A000045 := [combinat[fibonacci](n),n=0..100) ];
#     LODUMO(A000045,5) ;
#
# Richard J. Mathar, 2009-04-30
LODUMO := proc(L,m)
	local a,n, an;
	if not type(L,'list') then
		ERROR("Expected list type argument instead of ",whattype(L) ) ;
	end if;
	a := [] ;
	for n from 1 to nops(L) do
		for an from op(n,L) mod m by m do
			if not an in a then
				a := [op(a),an] ;
				break;
			fi;
		od:
	end do;
	a ;
end proc:

# Integer value of the concatenation a//b of the base-10 decimal representations of two integers.
# The arguments a and b are non-negative integers.
# @param a The integer that represents the leading decimal digits of the result
# @param b The integer that represents the trailing (low significant) decimal digits of the result
# @return The integer obtained by concatenation of the base-10 digits of a and b
# @author N. J. A. Sloane
# @since 2010-10-16
# Example:
#     digcat2(3,4)  ; # returns 34
#     digcat2(3,40)  ; # returns 340
#     digcat2(3,0)  ; # returns 30
#     digcat2(37,0)  ; # returns 370
#     digcat2(0,800)  ; # returns 800
#     digcat2(-4,8)  ; # error: negative value of an argument
digcat2 := proc(a,b)
	if a < 0 or b < 0 then
		ERROR("Negative argument ",a," or ", b) ;
	end if;
	parse(cat(a,b)) ;
end proc:

# Integer value of the concatenation of the base-10 decimal representations of the integers in a list.
# @param L The list of non-negative integers [a,b,c,d,...]
# @return The integer obtained by concatenation of the base-10 digits, a//b//c//...
# @since 2010-10-16
# @author N. J. A. Sloane
# Example:
#     digcatL([3,4])  ; # returns 34
#     digcatL([3,0,4])  ; # returns 304
#     digcatL([123,456])  ; # returns 123456
#     digcatL([3])  ; # returns 3
#     digcatL([])  ; # empty list returns 0
#     digcatL(3/4)  ; # error: argument not of the list type
digcatL := proc(L)
	local ret,i ;
	if not type(L,'list') then
		ERROR("Expected list type argument instead of ",whattype(L) ) ;
	end if;
	ret := 0 ;
	for i from 1 to nops(L) do
		ret := digcat2(ret,op(i,L)) ;
	end do;
	return ret ;
end proc:

# Aerate by insertion of zeros into the original list.
# The list L(1),0,0,0,..., L(2),0,0,.. is returned where the number of
# zeros spacing the original members of the list is controlled by a parameter.
# @param L the input list
# @param n the non-negative count of zeros to be inserted.
# @return 
# @since Richard J. Mathar, 2009-06-23
# @author Richard J. Mathar
AERATE := proc(L,n)
	local a,i,j;
	if not type(L,'list') then
		ERROR("Expected list type argument instead of ",whattype(L) ) ;
	end if;
	a := [] ;
	for i from 1 to nops(L) do
		a := [op(a),op(i,L),seq(0,j=1..n)] ;
	end do;
	a ;
end proc:

# Backwards Van Eck transform. Suggested by A181391.
# If a[n] has appeared earlier in a, let a[m] be the most recent 
# occurrence, and set b[n]=n-m; otherwise b[n]=0.
ECKb:=proc(a) local b,i,m,n;
if whattype(a) <> list then RETURN([]); fi:
b:=[0];
for n from 2 to nops(a) do
# has a(n) appeared before?
m:=0;
for i from n-1 by -1 to 1 do
	if (a[i]=a[n]) then m:=n-i; break; fi
od:
b:=[op(b),m];
od:
RETURN(b);
end proc:

# Forwards Van Eck transform. Suggested by A181391.
# If a[n] appears again in a, let a[m] be the next 
# occurrence, and set b[n]=m-n; otherwise b[n]=0.
# Warning: One should provide enough terms of a! 
ECKf:=proc(a) local b,i,m,n;
if whattype(a) <> list then RETURN([]); fi:
b:=[];
for n from 1 to nops(a)-1 do
# does a(n) appear again?
m:=0;
for i from n+1 to nops(a) do
	if (a[i]=a[n]) then m:=i-n; break; fi
od:
b:=[op(b),m];
od:
b:=[op(b),0];
RETURN(b);
end proc:


# LAGRANGE transform of a sequence {a(n)}
# Suggested by Lagrange's theorem that at most 4 squares are needed to sum to n.
# Returns b(n) = minimal number of terms of {a} needed to sum to n for 1 <= n <= M.
# C = maximal number of terms of {a} to try to build n
# M = upper limit on n
# Internally, the initial terms of both a and b are taken to be 0, but since this is a number-theoretic function, the output starts at n=1
# Example. To produce A002828:
# sq:=[seq(n^2, n=1..20)];
# LAGRANGE(sq,4,120); which produces
# [1, 2, 3, 1, 2, 3, 4, 2, 1, 2, 3, 3, ...
LAGRANGE:=proc(a, C, M)
local t1, ip, i, j, a1, a2, b, c, N1, N2, Nc;
if whattype(a) <> list then RETURN([]); fi:
# sort a, remove duplicates, include 0
t1:=sort(a);
a1:=sort(convert(convert(a, set), list));
if not member(0, a1) then a1:=[0, op(a1)]; fi;
N1:=nops(a1);
b:=Array(1..M+1, -1);
for i from 1 to N1 while a1[i]<=M do b[a1[i]+1]:=1; od;
a2:=a1; N2:=N1;
for ip from 2 to C do
c:={}:
   for i from 1 to N1 while a1[i] <= M do
      for j from 1 to N2 while a1[i]+a2[j] <= M do
c:={op(c), a1[i]+a2[j]};
                                                od;
                                       od;
c:=sort(convert(c, list));
Nc:=nops(c);
   for i from 1 to Nc do
      if b[c[i]+1] = -1 then b[c[i]+1]:= ip; fi;
                      od;
a2:=c; N2:=Nc;
                   od;
RETURN([seq(b[i], i=2..M+1)]);
end: 


myhead:=proc(s, L::integer := 100 ) local s2,n;
if whattype(s) <> list then s2:=convert(s,list); else s2:=s: fi;
[seq(s2[n],n=1..min(L,nops(s2)))];
end:

mytail:=proc(s, L::integer := 100 ) local s2,n;
if whattype(s) <> list then s2:=convert(s,list); else s2:=s: fi;
[seq(s2[n],n=max(1,nops(s2)-L+1)..nops(s2))];
end:


# myconvert(n,b) = convert(n,base,b) but puts the digits in human order
myconvert:=proc(n,b) local t1,t2,len,i;
if n=0 then RETURN([0]); fi;
t1:=convert(n,base,b); len:=nops(t1);
t2:=[seq(t1[len+1-i],i=1..len)];
end:


# Maple Print (prints a sequence on one long line) Added Feb 2018
MP:=proc(s)
interface(screenwidth=1200):
lprint(seq(s[i],i=1..nops(s)));
interface(screenwidth=79):
end:


# COMPl2:  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!!!
COMPl2:=proc(a,M)  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],M);
notmember:=subs(_WZ=b,proc(x) not member(x,_WZ) end);
select(notmember,[$1..m]);
end:


# binary: binary representation of n, in human order
binary:=proc(n) local t1,L;
if n<0 then ERROR("n must be nonnegative"); fi;
if n=0 then return([0]); fi;
t1:=convert(n,base,2); L:=nops(t1);
[seq(t1[L+1-i],i=1..L)];
end:




### %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
### 
### # Added Jan 05 2013: Michael Somos has very kindly supplied the
### # following code for checking these procedures.
### # To run it, save this is a file and strip off the ### characters from the beginning of each line
### 
### # Maple test code for OEIS transforms, Michael Somos 31DEC2012
### 
### read("transforms");
### 
### chk:=proc(ex)
###   if not evalb(parse(ex,statement)) then printf("%s: false\n",ex) end if;
### end:
### 
### chk("delta(x->x^3,2,2)=18");
### chk("did(8,3)=0");
### chk("did(9,3)=1");
### chk("dids(8,3)=0");
### chk("dids(9,3)=-1");
### chk("dim(26,10)=1");
### chk("dim(10,26)=0");
### chk("mob(10,2)=-1");
### chk("mob(10,3)=0");
### # 0 indexed MASKCONV
### chk("MASKCONV([1,2,3,4],[5,6,7,8])=[5,16,22,60]");
### # 0 indexed MASKTRANS
### chk("MASKTRANS([1,2,3,4])=[1,3,4,10]");
### chk("MASKTRANS([a0,a1,a2,a3,a4])=[a0,a0+a1,a0+a2,a0+a1+a2+a3,a0+a4]");
### # 0 indexed MASKTRANSi
### chk("MASKTRANSi([1,2,3,4])=[1,1,2,0]");
### chk("MASKTRANSi([a0,a0+a1,a0+a2,a0+a1+a2+a3,a0+a4])=[a0,a1,a2,a3,a4]");
### # 1 indexed SUPPORT
### chk("SUPPORT([0,1,2,0,5,4,0,-3])=[2,3,5,6,8]");
### chk("SUPPORT([a1,a2,0,a4,0,a6,a7,0,0,a10])=[1,2,4,6,7,10]");
### # no index MONO. numeric only in.
### chk("MONO([10,3,-2,-3,-10])=[2,3,10]");
### # MONO2 has "abs(bu - nops(b))<10" test at end. don't use!   TODO
### # 1 index out COMP. numeric only in.
### chk("COMP([6,4,2,10,7])=[1,3,5,8,9]");
### # 0 indexed CHAR. numeric only in.
### chk("CHAR([6,4,2,10,7])=[0,0,1,0,1,0,1,1,0,0,1]");
### # The following is a very peculiar function with "sequence" hardcoded
### sequence:=[[a],[b,c],[d],[x1,x2,x3,x4,x5,x6]]:
### # 1 indexed EXTRACT
### chk("EXTRACT(4,3,2)=[x2,x5]");
### # 0 indexed DECIMATE
### chk("DECIMATE([0,0,2,0,0,5,0],3,2)=[2,5]");
### chk("DECIMATE([a0,a1,a2,a3,a4,a5,a6],3,2)=[a2,a5]");
### # 0 indexed BISECT
### chk("BISECT([a0,a1,a2,a3,a4,a5,a6],0)=[a0,a2,a4,a6]");
### chk("BISECT([a0,a1,a2,a3,a4,a5,a6],1)=[a1,a3,a5]");
### # 0 indexed TRISECT
### chk("TRISECT([a0,a1,a2,a3,a4,a5,a6],0)=[a0,a3,a6]");
### chk("TRISECT([a0,a1,a2,a3,a4,a5,a6],1)=[a1,a4]");
### chk("TRISECT([a0,a1,a2,a3,a4,a5,a6],2)=[a2,a5]");
### # no index DIFF
### chk("DIFF([9,4,1,0,1,4,9])=[-5,-3,-1,1,3,5]");
### chk("DIFF([a,b,c,d,e,f])=[b-a,c-b,d-c,e-d,f-e]");
### # no index PRODS.  comment is wrong about  b[n] = a[1]*..*a[n]  TODO
### chk("PRODS([1,2,3,4,5])=[1,2,6,24,120]");
### chk("PRODS([a,b,c,d,e])=[a,a*b,a*b*c,a*b*c*d,a*b*c*d*e]");
### # no index PSUM.
### chk("PSUM([1,2,3,4,5])=[1,3,6,10,15]");
### chk("PSUM([a,b,c,d,e])=[a,a+b,a+b+c,a+b+c+d,a+b+c+d+e]");
### # no index PSUMSIGN.
### chk("PSUMSIGN([1,2,3,4,5])=[1,1,2,2,3]");
### chk("PSUMSIGN([a,b,c,d,e])=[a,b-a,c-b+a,d-c+b-a,e-d+c-b+a]");
### # 0 indexed BINOMIAL
### chk("BINOMIAL([1,2,4,8,16])=[1,3,9,27,81]");
### chk("BINOMIAL([a0,a1,a2,a3])=[a0,a0+a1,a0+2*a1+a2,a0+3*a1+3*a2+a3]");
### # 0 indexed BINOMIALi
### chk("BINOMIALi([1,3,9,27,81])=[1,2,4,8,16]");
### chk("BINOMIALi([a0,a0+a1,a0+2*a1+a2,a0+3*a1+3*a2+a3])=[a0,a1,a2,a3]");
### # 1 indexed BIN1
### chk("BIN1([a1,a2,a3,a4])=[a1,-2*a1-a2,3*a1+3*a2+a3,-4*a1-6*a2-4*a3-a4]");
### chk("BIN1([2,4,8,16])=[2,-8,26,-80]");
### # 1 indexed EULER.
### chk("EULER([1,1,0,0,0,0,0])=[1,2,2,3,3,4,4]");
### chk("simplify(EULER([a1,a2,a3]))=[a1,a2+a1/2+a1^2/2,a3+a1*a2+a1/3+a1^2/2+a1^3/6]");
### # 1 indexed EULERi
### chk("EULERi([1,2,2,3,3,4,4])=[1,1,0,0,0,0,0]");
### chk("simplify(EULERi([a1,a2+a1/2+a1^2/2,a3+a1*a2+a1/3+a1^2/2+a1^3/6]))=[a1,a2,a3]");
### # 1 indexed MOBIUS
### chk("MOBIUS([1,3,4,7,6,12])=[1,2,3,4,5,6]");
### chk("MOBIUS([a1,a2,a3,a4,a5,a6])=[a1,a2-a1,a3-a1,a4-a2,a5-a1,a6-a3-a2+a1]");
### # 1 indexed MOBIUSi
### chk("MOBIUSi([1,2,3,4,5,6])=[1,3,4,7,6,12]");
### chk("MOBIUSi([a1,a2-a1,a3-a1,a4-a2,a5-a1,a6-a3-a2+a1])=[a1,a2,a3,a4,a5,a6]");
### # 1 indexed STIRLING
### chk("STIRLING([1,2,3,4,5])=[1,3,10,37,151]");
### chk("STIRLING([a1,a2,a3,a4])=[a1,a1+a2,a1+3*a2+a3,a1+7*a2+6*a3+a4]");
### # 1 indexed STIRLINGi
### chk("STIRLINGi([1,3,10,37,151])=[1,2,3,4,5]");
### chk("STIRLINGi([a1,a1+a2,a1+3*a2+a3,a1+7*a2+6*a3+a4])=[a1,a2,a3,a4]");
### # 1 indexed STIRB
### chk("STIRB([1,2,3,4,5,6])=[1,-1,1,-1,1,-1]");
### chk("STIRB([a1,a2,a3,a4])=[a1,a1-a2,a1-3*a2+2*a3,a1-7*a2+12*a3-6*a4]");
### # 1 indexed INVERT. numeric only bug.  should be INVERT:=proc(b).    TODO
### chk("INVERT([1,2,3,4,5])=[1,3,8,21,55]");
### # 1 indexed INVERTi. numeric only bug.
### chk("INVERTi([1,3,8,21,55])=[1,2,3,4,5]");
### # 1 indexed REVERT. numeric only bug.
### chk("REVERT([1,2,3,4,5])=[1,-2,5,-14,42]");
### # 1 indexed REVEGF. numeric only bug.
### chk("REVEGF([1,2,3,4,5])=[1,-4,39,-616,13505]");
### # 1 indexed COMPOSE.
### chk("COMPOSE([4,3,2,1],[1,2,3,4])=[4,11,26,59]");
### chk("COMPOSE([a1,a2,a3],[b1,b2,b3])=[a1*b1,a1*b2+a2*b1^2,a1*b3+2*a2*b1*b2+a3*b1^3]");
### # 0 indexed CONV.
### chk("CONV([4,3,2,1,0],[0,1,2,3,4])=[0,4,11,20,30]");
### chk("CONV([a0,a1,a2],[b0,b1,b2])=[a0*b0,a0*b1+a1*b0,a0*b2+a1*b1+a2*b0]");
### # 0 indexed CONVi.
### chk("CONVi([1,4,10,20,35])=[1,2,3,4,5]");
### chk("CONVi([a0*a0,2*a0*a1,2*a0*a2+a1*a1])=[a0,a1,a2]");
### # 0 indexed LISTTOLSITDIV.
### chk("LISTTOLISTDIV([1,4,9,16,25])=[1,4,9/2,8/3,25/24]");
### chk("LISTTOLISTDIV([a0,a1,a2,a3,a4])=[a0,a1,a2/2,a3/6,a4/24]");
### # 0 indexed LISTTOLSITMUL.
### chk("LISTTOLISTMULT([1,4,9,16,25])=[1,4,18,96,600]");
### chk("LISTTOLISTMULT([a0,a1,a2,a3,a4])=[a0,a1,2*a2,6*a3,24*a4]");
### # no index PARTITION. 1 indexed out. numeric input only okay.
### chk("PARTITION([1,3,5,13])=[1,1,2,2,3,4,4,5,6,7,8,9,11]");
### # 1 indexed PARITITION1. numeric input only okay.
### chk("PARTITIONi([1,1,2,2,3,4,4,5,6,7,8,9,11])=[1,3,5,13]");
### # 1 indexed EXP.
### chk("EXP([1,2,3,4])=[1,3,10,41]");
### chk("EXP([a1,a2,a3])=[a1,a2+a1*a1,a3+3*a1*a2+a1*a1*a1]");
### # 1 indexed LOG.
### chk("LOG([1,3,10,41])=[1,2,3,4]");
### chk("simplify(LOG([a1,a2+a1*a1,a3+3*a1*a2+a1*a1*a1]))=[a1,a2,a3]");
### # 0 indexed EXPCONV.
### chk("EXPCONV([1,2,3,4],[4,3,2,1])=[4,11,26,56]");
### chk("EXPCONV([a0,a1,a2],[b0,b1,b2])=[a0*b0,a0*b1+a1*b0,a0*b2+2*a1*b1+a2*b0]");
### # no index LEFT.
### chk("LEFT([4,3,2,1])=[3,2,1]");
### chk("LEFT([a,b,c,d,e])=[b,c,d,e]");
### # no index RIGHT.
### chk("RIGHT([4,3,2,1])=[1,4,3,2,1]");
### chk("RIGHT([a,b,c,d,e])=[1,a,b,c,d,e]");
### # no index RUNS.
### chk("RUNS([4,3,2,2,2,1,3,3,4])=[1,1,3,1,2,1]");
### chk("RUNS([a,b,c,c,c,d,b,b,a])=[1,1,3,1,2,1]");
### # 0 indexed M2.
### chk("M2([5,4,3,2,1])=[5,8,6,4,2]");
### chk("M2([a,b,c,d,e])=[a,2*b,2*c,2*d,2*e]");
### # 0 indexed M2i.
### chk("M2i([5,8,6,4,2])=[5,4,3,2,1]");
### chk("M2i([a,2*b,2*c,2*d,2*e])=[a,b,c,d,e]");
### # 0 indexed NEGATE.
### chk("NEGATE([5,4,3,2,1])=[5,-4,-3,-2,-1]");
### chk("NEGATE([a,b,c,d,e])=[a,-b,-c,-d,-e]");
### # 0 indexed eultrans2. (=BINOMIALi)
### chk("eultrans2([1,1,2,3,5,8,13])=[1,0,1,-1,2,-3,5]");
### chk("eultrans2([a0,a1,a2,a3])=[a0,a1-a0,a2-2*a1+a0,a3-3*a2+3*a1-a0]");
### # 1 indexed Etrans.
### chk("Etrans([1,1,0,0,0,0,0])=[1,2,2,3,3,4,4]");
### chk("simplify(Etrans([a1,a2,a3]))=[a1,a2+a1/2+a1^2/2,a3+a1*a2+a1/3+a1^2/2+a1^3/6]");
### # no index pairtrans.
### chk("pairtrans([5,4,3,2,1])=[5,9,7,5,3]");
### chk("pairtrans([a,b,c,d,e])=[a,a+b,b+c,c+d,d+e]");
### # 1 indexed BOUS in. 0 indexed out. b0=1 always.
### chk("BOUS([5,4,3,2,1])=[1,6,15,32,83,262]");
### chk("BOUS([a1,a2,a3,a4])=[1,1+a1,1+2*a1+a2,2+3*a1+3*a2+a3,5+8*a1+6*a2+4*a3+a4]");
### # 0 indexed BOUS2.
### chk("BOUS2([5,4,3,2,1])=[5,9,16,33,84]");
### chk("BOUS2([a0,a1,a2,a3])=[a0,a1+a0,a2+2*a1+a0,a3+3*a2+3*a1+2*a0]");
### # 0 indexed BOUS2i.
### chk("BOUS2i([5,4,3,2,1])=[5,-1,0,-5,4]");
### chk("BOUS2i([a0,a1,a2,a3])=[a0,a1-a0,a2-2*a1+a0,a3-3*a2+3*a1-2*a0]");
### # 1 indexed CONTINUANT.
### chk("CONTINUANT([5,4,3,2,1])=[5,21,68,157,225]");
### chk("CONTINUANT([a,b,c,d])=[a,a*b+1,a*b*c+a+c,a*b*c*d+a*b+a*d+c*d+1]");
### # 1 indexed CONTINUANTi.
### chk("CONTINUANTi([5,21,68,157,225])=[5,4,3,2,1]");
### chk("CONTINUANTi([a,b,c,d])=[a,b/a-1/a,c/b-a/b,d/c-b/c]");
### # 1 indexed DIRICHLET.
### chk("DIRICHLET([5,4,3,2,1],[1,2,3,4,5])=[5,14,18,30,26]");
### chk("DIRICHLET([a1,a2,a3,a4],[b1,b2,b3,b4])=[a1*b1,a1*b2+a2*b1,a1*b3+a3*b1,a1*b4+a2*b2+a4*b1]");
### # 1 indexed DIRICHLETi.
### chk("DIRICHLETi([1,2,3,4,5,6])=[1,-2,-3,0,-5,6]");
### chk("DIRICHLETi([1/a,b,c,d,e])=[a,-a^2*b,-a^2*c,a^3*b^2-a^2*d,-a^2*e]");
### # no indexed lHANKEL.
### chk("lHANKEL([1,4,9,16,25,36])=[7,17,31,49]");
### chk("lHANKEL([a,b,c,d,e])=[b*b-a*c,c*c-b*d,d*d-c*e]");
### # no indexed HANKEL.
### chk("HANKEL([6,5,4,3,2,1])=[6,-1,0]");
### chk("HANKEL([a,b,c,d,e,f])=[a,a*c-b*b,a*c*e-a*d*d-b*b*e+2*b*c*d-c*c*c]");
### # no index SERIES2TOLIST.
### chk("SERIES2TOLIST(a+b*x+c*y+d*x*x+e*x*y+f*y*y,x,y,3)=[a,b,c,d,e,f]");
### # no index SERIES2TOLISTMULT.
### chk("SERIES2TOLISTMULT(a+b*x+c*y+d*x*x+e*x*y+f*y*y,x,y,3)=[a,b,c,2*d,e,2*f]");
### # no index RIGHT0.
### chk("RIGHT0([5,4,3,2,1])=[0,5,4,3,2,1]");
### chk("RIGHT0([a,b,c,d,e])=[0,a,b,c,d,e]");
### # 0 indexed SUMTABL.
### chk("SUMTABL([4,3,2,1],[1,2,3,4])=[5,4,6,3,5,7,2,4,6,8]");
### chk("SUMTABL([a0,a1,a2],[b0,b1,b2])=[a0+b0,a1+b0,a0+b1,a2+b0,a1+b1,a0+b2]");
### # 0 indexed LAH.
### chk("LAH([5,4,3,2,1])=[5,4,11,44,229]");
### chk("LAH([a0,a1,a2,a3,a4])=[a0,a1,a2+2*a1,a3+6*a2+6*a1,a4+12*a3+36*a2+24*a1]");
### # 0 indexed LAHi.
### chk("LAHi([5,4,3,2,1])=[5,4,-5,8,-11]");
### chk("LAHi([a0,a1,a2,a3,a4])=[a0,a1,a2-2*a1,a3-6*a2+6*a1,a4-12*a3+36*a2-24*a1]");
### # 1 indexed INVERSE.  numeric input only.       RETURN(b) missing bug
### chk("INVERSE([4,5,3,2,1])=[5,4,3,1,2]");
### # no index digcatL. numeric input only.
### chk("digcatL([5,43,210])=543210");
### # 0 indexed AERATE.
### chk("AERATE([5,4,3,2,1],1)=[5,0,4,0,3,0,2,0,1,0]");
### chk("AERATE([a,b,c,d],2)=[a,0,0,b,0,0,c,0,0,d,0,0]");
### 
### printf("done");
### 
### %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%