This site is supported by donations to The OEIS Foundation.
User:R. J. Mathar/transforms3
From OeisWiki
# This is a file with Maple programs. Usage from within another Maple program:
# read("transforms3") ;
# Read the b-file, named in the argument, and return the list of the
# numbers of its 2nd column. The information on the offset of the sequence
# is not transmitted: the first number in the list that is returned is just
# the topmost value in the b-file, which may be a(0), a(1) or another
# entry depending on the offset.
#
# Example to get a list of prime numbers after having downloaded b000040.txt
# from the OEIS web page:
# a := BFILETOLIST("b000040.txt") ;
#
# @param bfilename The string with the input file which is read/scanned.
# @since 2010-09-06 Delete non-ASCII characters in the input
# Richard J. Mathar, 2010-02-20
BFILETOLIST := proc(bfilename)
local a,f,bline,b,nold,n ;
# List intially empty
nold := -999999 ;
a := [] ;
f := fopen(bfilename,READ) :
bline := readline(f) :
# Loop over the entire b-file
while bline <> 0 do
# There have been files with non-ASCII characters in the past
# (BOM markers of UTF codes in particular) which are difficult to detect.
# The following is an attempt to treat this error case gently,
# deleting non-ASCII bytes completely.
bline := StringTools[Select](StringTools[IsASCII],bline) ;
bline := StringTools[TrimLeft](bline) ;
# Skip lines with comments
if StringTools[FirstFromLeft]("#",bline ) = 1 then
;
else
# interpret the line contents as two integers (index and value)
b := sscanf( bline,"%d %d") ;
# grab the second component of each line, the value; append
if nops(b) >1 then
a := [op(a),op(2,b)] ;
n := op(1,b) ;
if n <> nold+1 and nold <> -999999 then
printf("# discontinuous first column %d %d\n",nold,n) ;
end if;
nold := n ;
end if;
end if ;
bline := readline(f) ;
end do:
fclose(f) ;
RETURN(a) ;
end proc:
# Write the list "L" to the file "bfilename" with the first number
# given the offset "offs".
#
# Warning: existing files are overwritten.
#
# Example:
# L := [seq(ithprime(i),i=1..20)] ;
# LISTTOBFILE("B000040.txt",L,1) ; # note capital B in the example, not in real life.
#
# Richard J. Mathar, 2008-03-12
LISTTOBFILE := proc(bfilename,L,offs)
local n,f ;
f := fopen(bfilename,WRITE) :
# Loop over the list members
for n from 1 to nops(L) do
fprintf(f,"%d %d\n",n+offs-1, op(n,L) ) ;
end do:
fclose(f) ;
printf("%%I A%6s\n", substring(bfilename,2..7) ) ;
printf("%%H A%6s F. Lastname, <a href=\"/A%6s/%s\">Table of n, a(n) for n=%d..%d</a>\n",substring(bfilename,2..7),substring(bfilename,2..7), bfilename,offs, nops(L)+offs-1) ;
RETURN() ;
end proc:
# Read a b-file and replace tabs by spaces and <cr><lf> combinations
# by <lf>. Remove empty lines.
# The argument bfilename should be a readable b-file name b??????.txt
# (no directory component). The result will be placed in a file named
# B??????.txt (so that's just the first letter replaced by a capital-B), which is
# overwritten if existing.
# Richard J. Mathar, 2009-03-10
BFILE2BFILE := proc(bfilename)
local f,fo,bline,b,n,boutname ;
f := fopen(bfilename,READ) :
boutname := sprintf("B%6s.txt", substring(bfilename,2..7) ) ;
fo := fopen(boutname,WRITE) :
bline := readline(f) :
# Loop over the entire b-file; eliminate trailing blank lines
while bline <> 0 do
# copy lines with comments
if StringTools[FirstFromLeft]("#",bline ) <> 0 then
fprintf(fo,"%s\n",bline ) ;
else
# interpret the line contents as two integers (index and value)
b := sscanf( bline,"%d %d") ;
if nops(b) >=2 then
# grab the second component of each line, the value; append
fprintf(fo,"%d %d\n",op(1,b), op(2,b) ) ;
fi;
end if ;
bline := readline(f) ;
end do:
fclose(f) ;
fclose(fo) ;
RETURN() ;
end:
# Convert a floating point constant to a sequence of decimal digits.
# @param c the floating point constant
# @return the list with zeros retained according to OEIS guidelines.
# The length of the list depends on the value of Digits
# at the time of the call. Two digits are chopped off as we assume
# that the value of the constant would typically be rounded, not
# truncated.
#
# Examples:
# CONSTTOLIST(3.14159) ; # returns [3,1,4,1,5,9,0,0,...]
# CONSTTOLIST(59) ; # returns [5,9,0,0,0,..]
#
# Richard J. Mathar, 2008-09-19
CONSTTOLIST := proc(c)
local absc,a,dgs;
# ensure it is the positive value, if negative
absc := abs(c) ;
a := [] ;
# for numbers < 0.1, start with the zeros to maintain offset 0
while absc < 0.1 do
a := [op(a),0] ;
absc := 10*absc ;
end do:
# dgs is the number of digits of absc before the decimal point, the number that remains to be printed
dgs := 1+ilog10(absc) ;
# dgs is the number of places we want to shift absc to the right, truncating two more digits than DIGITS
dgs := Digits-2-dgs ;
absc := floor(absc*10^dgs) ;
a := [op(a),op(ListTools[Reverse](convert(absc,base,10)))] ;
RETURN(a) ;
end proc:
# Convert a floating point constant to a sequence of digits in base b.
# @param c the floating point constant
# @param b the base (2, 3,... )
# @return the list with zeros retained according to OEIS guidelines.
# The length of the list depends on the value of Digits
# at the time of the call. Two digits are chopped off as we assume
# that the value of the constant would typically be rounded, not
# truncated.
#
# Examples:
# CONSTTOLIST(evalf(Pi/3),2) ; # returns 1,0,0,0,0,1,... (See A019670)
# CONSTTOLIST(59,2) ;
#
# Richard J. Mathar, 2009-02-06
CONSTTOLISTB := proc(c,b)
local absc,a,dgs;
# ensure it is the positive value, if negative
absc := abs(c) ;
a := [] ;
# for numbers < 0.1, start with the zeros to maintain offset 0
while absc < 1/b do
a := [op(a),0] ;
absc := b*absc ;
end do:
# dgs is the number of digits of absc before the decimal point, the number that remains to be printed
dgs := 1+floor(log(absc)/log(b)) ;
# dgs is the number of places we want to shift absc to the right, truncating two more digits than DIGITS
dgs := Digits-2-dgs ;
absc := floor(absc*b^dgs) ;
a := [op(a),op(ListTools[Reverse](convert(absc,base,b)))] ;
RETURN(a) ;
end proc:
# Convert a sequence to a floating point constant.
# The arguments are the list and the (correct) offset.
# Richard J. Mathar, 2009-02-02
LISTTOCONST := proc(L,offs)
local c,n,dgs,realdigs ;
c := 0 ;
dgs := 0 ;
for n from 1 to nops(L) do
if dgs = 0 then
if op(n,L) <> 0 then
dgs := 1 ;
end if;
else
dgs := dgs+1 ;
end if;
c := 10*c+op(n,L) ;
end do:
# Set effective precision for the floating point conversion to
# what the list's length (after removing leading zeros) offers.
realdigs := Digits ;
Digits := dgs ;
c := evalf( c*10^(offs-nops(L)) ) ;
Digits := realdigs ;
RETURN(c) ;
end proc:
# Motzkin transform of the sequence provided by the list L.
# This is equivalent to replacing the independent variable x in the
# generating function for L by x*M(x), where M(x) is the generating
# function for the Motzkin numbers.
#
# Examples:
# L := [0,1,0,0,0,0,0,0,0] ; Motzkin(L) ; # returns [0,1,1,2,4,9,21,51,27]
# L := [0,1,0,1,0,1,0,1,0,1] ; Motzkin(L) ; # returns [0,1,1,3,7,19,51,141,393,1107,...]
# L := [1,0,1,0,1,0,1,0,1] ; Motzkin(L) ; # returns [1,0,1,2,6,16,45,126,357,1016,2907,...]
#
# Richard J. Mathar, 2008-11-07
MOTZKIN := proc(L)
local d, bgf, mL, maxd;
# output list intially empty
# maximum degree of output is limited by initial length of L
maxd := nops(L)-1 ;
# generate an auxiliary list of Motzkin numbers
mL := taylor( (1-x-sqrt(1-2*x-3*x^2))/2/x, x=0,maxd+2) ;
bgf := op(1,L) ;
for d from 1 to maxd do
bgf := bgf+taylor(op(d+1,L)*mL^d,x=0,maxd+1) ;
end do:
[seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:
# Inverse Motzkin transform of the sequence provided by the list L.
# The inverse of the operation given by MOTZKIN above.
#
# Richard J. Mathar, 2008-11-08
MOTZKINi := proc(L)
local d, bgf, mL, maxd;
# output list intially empty
# maximum degree of output is limited by initial length of L
maxd := nops(L)-1 ;
# generate an auxiliary list of Motzkin numbers
mL := taylor( x/(1+x+x^2), x=0,maxd+2) ;
bgf := op(1,L) ;
for d from 1 to maxd do
bgf := bgf+taylor(op(d+1,L)*mL^d,x=0,maxd+1) ;
end do:
[seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:
# Catalan transform of the sequence provided by the list L.
# This is equivalent to replacing the independent variable x in the
# generating function for L by x*C(x), where C(x) is the generating
# function for the Catalan numbers.
#
# Richard J. Mathar, 2008-12-17
CATALAN := proc(L)
local d, bgf, cL, maxd;
# output list intially empty
# maximum degree of output is limited by initial length of L
maxd := nops(L)-1 ;
# generate an auxiliary list of Catalan numbers
cL := taylor( (1-sqrt(1-4*x))/2, x=0,maxd+2) ;
bgf := op(1,L) ;
for d from 1 to maxd do
bgf := bgf+taylor(op(d+1,L)*cL^d,x=0,maxd+1) ;
end do:
[seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:
# Inverse Catalan transform of the sequence provided by the list L.
# The inverse of the operation given by CATALAN above.
# Richard J. Mathar, 2008-12-11
CATALANi := proc(L)
local d, bgf, cL, maxd;
# output list intially empty
# maximum degree of output is limited by initial length of L
maxd := nops(L)-1 ;
cL := taylor( x*(1-x), x=0,maxd+2) ;
bgf := op(1,L) ;
for d from 1 to maxd do
bgf := bgf+taylor(op(d+1,L)*cL^d,x=0,maxd+1) ;
end do:
# bug in Maple 9: the following line is necessary to get correct results
bgf := taylor(bgf,x=0,maxd+1) ;
[seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:
# Pierce sequence generated from a constant.
# The input is a number between 0 and 1, the output the list of Pierce
# expansion coefficients.
# Richard J. Mathar, 2009-01-21
PIERCE := proc(c)
local resid,L,i,an ;
resid := c ;
L := [] ;
if whattype(c) = `float` then
if c >0 and c <= 1 then
for i from 1 do
an := floor(1./resid) ;
L := [op(L),an] ;
resid := evalf(1.-an*resid) ;
if ilog10( mul(i,i=L)) > 0.7*Digits then
break ;
end if ;
end do:
end if;
end if;
L ;
end proc:
# Engel sequence generated from a constant.
# The input is a number between 0 and 1, the output the list of Engel
# expansion coefficients.
# Richard J. Mathar, 2009-07-30
ENGEL := proc(c)
local resid,L,i,an ;
resid := c ;
L := [] ;
if whattype(c) = `float` then
if c > 1 then
RETURN( [1,op(procname(c-1.))] ) ;
elif c >0 and c <= 1 then
for i from 1 do
an := ceil(1./resid) ;
L := [op(L),an] ;
resid := evalf(an*resid-1.) ;
if ilog10( mul(i,i=L)) > 0.7*Digits then
break ;
end if ;
end do:
end if;
end if;
L ;
end proc:
# Beatty sequence generated from a constant.
# Richard J. Mathar, 2009-05-20
BEATTY := proc(c)
local L,n,x ;
L := [] ;
x := evalf(c) ;
for n from 1 to 0.8*Digits do
L := [op(L),floor(n*x)] ;
end do:
L ;
end proc:
# L-Schroeder transform of the sequence provided by the list L.
# This is equivalent to replacing the independent variable x in the
# generating function for L by x*S(x), where S(x) is the generating
# function for the large Schroeder numbers.
#
# Richard J. Mathar, 2008-11-13
SCHROEDER := proc(L)
local d, bgf, mL, maxd;
# output list intially empty
# maximum degree of output is limited by initial length of L
maxd := nops(L)-1 ;
# generate an auxiliary list of Schroeder numbers
mL := taylor( (1-x-sqrt(1-6*x+x^2))/2, x=0,maxd+2) ;
bgf := op(1,L) ;
for d from 1 to maxd do
bgf := bgf+taylor(op(d+1,L)*mL^d,x=0,maxd+1) ;
end do:
[seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:
# Inverse L-Schroeder transform of the sequence provided by the list L.
# This is equivalent to replacing the independent variable x in the
# generating function for L by x*(1-x)/(1+x).
#
# Example: if L is A001850, the funtion returns a list of 1 followed
# by A001333.
# Richard J. Mathar, 2008-11-13
SCHROEDERi := proc(L)
local d, bgf, mL, maxd;
# output list intially empty
# maximum degree of output is limited by initial length of L
maxd := nops(L)-1 ;
# generate an auxiliary list of Schroeder numbers
mL := taylor( (1-x)*x/(1+x), x=0,maxd+2) ;
bgf := op(1,L) ;
for d from 1 to maxd do
bgf := bgf+taylor(op(d+1,L)*mL^d,x=0,maxd+1) ;
end do:
[seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:
# Reduce the list given by the argument L to a list in which immediately repeated
# numbers (ie, clusters of the same numbers) are reduced to a single occurrence.
#
# Example: if L contains [3,5,6,7,7,8,3,3,3,10],
# then the list returned is [3,5,6,7,8,3,10], where "7,7" has shrunk to "7"
# and "3,3,3" has shrunk to "3".
#
# Richard J. Mathar, 2008-05-03
LISTNODUPL := proc(L)
local b, a, aprev;
# output list intially empty
b := [] ;
# Note that aprev is not initialized to avoid a coincidental equality with L[1].
for a in L do
if a <> aprev then
b := [op(b),a] ;
end if ;
aprev := a ;
end do:
RETURN(b) ;
end proc:
# Reduce the list given by the argument "L" to a list which only contains clusters
# of 2 or more repeated numbers. The argument "multpl" is either a number
# larger than 0 which specifies the multiplicity (cluster length) which the
# repetitions need to have (exactly) to enter the new list, or is 0 (zero)
# to indicate that any cluster length of 2 or more will have one representative
# be put on the list returned.
#
# Example: If L contains [3,5,7,7,8,3,3,3,10,7,7,7,11]
# then the list returned by LISTDUPL(L,2) is [7], which represents the
# first of the two clusters with 7's.
# Example: If L contains [3,5,7,7,8,3,3,3,10,7,7,7,11]
# then the list returned by LISTDUPL(L,3) is [3,7], which represents the
# cluster with the 3's and the 2nd cluster with the 7's, each showing
# numbers repeated 3 times.
# Example: If L contains [3,5,7,7,8,3,3,3,10,7,7,7,11,11]
# then the list returned by LISTDUPL(L,0) is [7,3,7,11], which represents the
# cluster with the two 7's, the 2nd cluster with the 3's, the 2nd cluster
# with 7's, and the final cluster with 11's.
# Example: If L contains [3,5,7,7,8,3,3,3,10,7,7,7,11,11]
# then the list returned by LISTDUPL(L,1) is [3,5,8,10], which represents the
# isolated numbers (numbers in clusters of length 1).
#
# To obtain results which register repeated numbers without regard to
# their placement in the original list, these may be sorted before calling
# the routine: LISTDUPL(sort(L),multpl) ;
#
# Richard J. Mathar, 2008-05-03
LISTDUPL := proc(L,multpl)
local b, a, aprev, cnt;
if multpl < 0 then
ERROR("negative multiplicity ",args[2]) ;
end if ;
# output list intially empty
b := [] ;
cnt :=0 ;
aprev := op(1,L) ;
for a in L do
if a <> aprev then
if ( cnt = multpl and multpl >0 ) or (cnt > 1 and multpl = 0 ) then
b := [op(b),aprev] ;
end if ;
cnt := 1 ;
else
cnt := cnt+1 ;
end if ;
aprev := a ;
end do:
# At that time, aprev represents a number in the last cluster. Since
# the lists in the OEIS are generally incomplete at the end, the following
# check on the length of this last cluster is not executed, because
# the length of this cluster might actually be larger than the list
# representatives tell; this may lead to misinterpretations if the lists
# are "full"!
# if ( cnt = multpl and multpl >0 ) or (cnt > 1 and multpl = 0 ) then
# b := [op(b),aprev] ;
# end if ;
if cnt > 1 and multpl = 0 then
b := [op(b),aprev] ;
end if ;
RETURN(b) ;
end proc:
# Lehmer's cotangent representation of an irrational number.
# Example
# CONSTTOCOT(evalf(Pi)); # see A002667
# CONSTTOCOT(evalf(log(2))); # see A081785
# Fritz Schweiger, Metrische Ergebnisse ueber den Kotangensalgorithmus, Acta Arith. 26 (1975) 217
CONSTTOCOT := proc(x)
local L,xcut ;
xcut := x ;
L := [floor(xcut)] ;
while type(L[-1],'infinity') = false do
arccot(L[-1])-arccot(xcut) ;
xcut := cot(evalf(%) ) ;
L := [op(L),floor(xcut)] ;
end do;
return subsop(-1=NULL,L) ;
end proc:
# LOG transformation operation on o.g.f.
# The input list L represents [a0=1,a1,a2,a3,a4....] of an ordinary
# generating function, g=1+a1*x+a2*x^2+.... This is to say that the first
# element of L is tacitly assumed to be equal to 1.
#
# The output contains the coefficients of log(g) = sum_{i=1..inf} b(i)*x^i/i .
# Richard J. Mathar, 2009-03-13
OLOG := proc(L)
local x,gf,gfl,i,a;
gf := add( op(i,L)*x^(i-1),i=2..nops(L)) ;
gfl := 0 ;
for i from 1 to nops(L) do
gfl := taylor(gfl+gf^i/i,x=0,nops(L)) ;
end do:
a := [] ;
for i from 1 to nops(L)-1 do
a := [op(a), i*coeftayl(gfl,x=0,i)] ;
end do:
a;
end proc:
# Akiyama-Tanigawa algorithm, one step
# The input list L with L[0], L[1] is transformed to the output list
# A[0], A[1] as described in
# D. Merlini, R. Sprugnoli, M. C. Verri, The Akiyama-Tanigawa Transformation, Integers, 5 (1) (2005) #A05
# Richard J. Mathar 2010-12-02
AKIYATANI := proc(L)
local a,k ;
a := [] ;
for k from 1 to nops(L)-1 do
a := [op(a),k*(op(k,L)-op(k+1,L))] ;
end do;
a ;
end proc:
# Akiyama-Tanigawa algorithm, generate first column.
# The input list L is iteratively transformed to construct the first columnof the table.
# D. Merlini, R. Sprugnoli, M. C. Verri, The Akiyama-Tanigawa Transformation, Integers, 5 (1) (2005) #A05
# Richard J. Mathar 2010-12-02
AKIYAMATANIGAWA := proc(L)
local a,r,img ;
a := [op(1,L)] ;
img := L ;
for r from 1 to nops(L)-1 do
img := AKIYATANI(img) ;
a := [op(a),op(1,img)] ;
end do;
a ;
end proc:
# Inverse Akiyama-Tanigawa tabulation: generate top row given the column k=0.
# Richard J. Mathar 2010-12-02
AKIYAMATANIGAWAi := proc(L)
local a,r,img,k,row ;
img := [op(-1,L)] ;
for r from nops(L)-1 to 1 by -1 do
row := [op(r,L)] ;
for k from 2 to 1+nops(L)-r do
op(k-1,row)-op(k-1,img)/(k-1) ;
row := [op(row),%] ;
end do;
img := row ;
end do;
img ;
end proc:
# Given the left entry d(x) and the right entry h(x) of the Riordan array,
# compute the (n,k) element of the array.
# @param d a function of x
# @param h a function of x
# @param n the row index >=0
# @param k the column index >=0
# Richard J. Mathar 2011-01-09
RIORDAN := proc(d,h,n,k)
d*h^k ;
expand(%) ;
coeftayl(%,x=0,n) ;
end proc:
RIORDANEXP := proc(d,h,n,k)
RIORDAN(d,h,n,k)*n!/k! ;
end proc:
# Generate a list of multinomial partitions.
# @param n The sum of the first moment of the parts, sum_j j*p_j = n
# @param k The sum of the parts, sum_j p_j = k
# @param nparts The number of parts, sum_j 1 = nparts
# @param zIsPart If true, the constrained is 0<=p_j, other wise 1<=p_j.
# @return A list of the form [[p_1,p_2,p_3,...],[p_1,p_2,p_3,...]]
# where each of the sublists is a partition in nparts obeying the two constraints.
# There no other restrictions but p_j<=k. Some of the p_j may be equal.
# Richard J. Mathar, 2011-04-06
multinomial2q := proc(n::integer,k::integer,nparts::integer,zIsPart::boolean)
local lpar ,res, constrp;
res := [] ;
if n< 0 or nparts <= 0 then
;
elif nparts = 1 then
if n = k then
if n > 0 or zIsPart then
return [[n]] ;
end if;
end if;
else
for lpar from `if`(zIsPart,0,1) do
if lpar*nparts > n or lpar > k then
break;
end if;
# recursive call, assuming that lpar is the last part
for constrp in procname(n-nparts*lpar,k-lpar,nparts-1,zIsPart) do
if nops(constrp) > 0 then
res := [op(res),[op(constrp),lpar]] ;
end if;
end do:
end do:
end if ;
return res ;
end proc:
# Generate a list of multinomial partitions.
# @param n The sum of the first moment of the parts, sum_j j*p_j = n
# @param k The sum of the parts, sum_j p_j = k
# @param zIsPart If true, the constrained is 0<=p_j, other wise 1<=p_j.
# @return A list of the form [[p_1,p_2,p_3,...],[p_1,p_2,p_3,...]]
# where each of the sublists is a partition obeying the two constraints.
# There are no other restrictions but 1<=p_j<=k. Some of the p_j may be equal.
# The Bell polynomials of the second kind are an application of this.
# Note that ordering in the sublists matters.
# Empty partitions (in case n=k=0) are not returned as empty sublists but not returned at all.
# Example:
# multinomial2(10,5,false) returns [[1,3,1],[2,1,2]] where 1+3+1=5, 1+2*3+3*1= 10
# for the first sublist and 2+1+2=5, 2+2*1+3*2=10 for the second sublist.
# Richard J. Mathar, 2011-04-06
multinomial2 := proc(n::integer,k::integer,zIsPart::boolean)
local npart,res,constrp ;
res := [] ;
for npart from 1 do
# The smallest part is 1. So sum_j j*p_j >= npart*(npart+1)/2 =n means
# that the number of parts is restricted by n, only to a lesser degree
# by k.
if npart*(npart+1)/2 > n or npart > k then
return res;
end if;
# loop over the partitions with fixed number of parts and
# attach those solultions to the result vector
for constrp in multinomial2q(n,k,npart,zIsPart) do
if nops(constrp) > 0 then
res := [op(res),constrp] ;
end if ;
end do:
end do:
end proc:
# Bell polynomial of the second kind.
# @param n
# @param k
# @param x The vector of x1,x2,x3...
# @return sum_{j1+j2+..=k, j1+2*j2+.. =n} n!/(j1!j2!..) (x1/1!)^j1 (x2/2!)^j2 ...
# Richard J. Mathar, 2011-04-06
bell2 := proc(n::integer,k::integer,x::list)
local p,npart,b,m,j;
npart := nops(x) ;
b := 0 ;
for p in multinomial2q(n,k,npart,true) do
m := 1 ;
for j from 1 to nops(p) do
m := m*(op(j,x)/j!)^op(j,p)/op(j,p)! ;
end do:
b := b+ n!*m ;
end do:
return b;
end proc:
# Seidel tranformation with parameters alpha, beta.
# @param L Input sequence
# @param alpha First multiplier
# @param beta Second multiplier
# @return The left column T(n,0) of the matrix T(0,i)=L(i), T(n,m) = alpha*T(n-1,m)+beta*T(n-1,m+1).
# The case alpha=beta=1 is just the binomial transform.
# Richard J. Mathar, 2012-10-10
SEIDEL := proc(L,f1,f2)
local a,n ;
a := [] ;
n := nops(L) ;
for n from 0 to nops(L)-1 do
add(f1^(n-k)*f2^kbinomial(n,k)*op(k+1,L),k=0..n) ;
a := [op(a),%]
end do:
a ;
end proc:
# Inverse Seidel tranformation with parameters alpha, beta.
# @param L Input sequence
# @param alpha First multiplier
# @param beta Second multiplier
# @return The top row T(0,n) of the matrix T(i,0)=L(i), T(n,m) = alpha*T(n-1,m)+beta*T(n-1,m+1).
# Richard J. Mathar, 2012-10-13
SEIDELi := proc(L,f1,f2)
SEIDEL(L,-f1/f2,1/f2) ;
end proc:
Digits := 500 ;
# Digits := 500 ;
# First Nmax terms of the signature sequence of x.
# @param x The floating point constant to be rendered.
# @param Nmax The number of terms of the signature sequence to be generated.
# @return A list of the first Nmax terms of the signature sequence.
# Richard J. Mathar, 2013-07-30
SignSeq := proc(x,Nmax)
local Lij,i,j,k,ts,ins ;
Lij := [] ;
for i from 1 to Nmax do
for j from 1 to Nmax do
ts := i+j*x ;
ins := 0 ;
for k from 1 to nops(Lij) do
if ts > op(1,op(k,Lij))+op(2,op(k,Lij))*x then
ins := k ;
end if;
end do:
if ins > Nmax then
;
elif ins = 0 then
Lij := [[i,j],op(Lij)] ;
elif ins = nops(Lij) then
Lij := [op(Lij),[i,j]] ;
else
Lij := [op(1..ins,Lij),[i,j],op(ins+1..nops(Lij),Lij)] ;
end if;
end do:
end do:
[seq(op(1,op(k,Lij)),k=1..Nmax)] ;
end proc:
# Deleham's Delta function in A084938
Delta := proc(r::list,s::list,n,k)
option remember;
local P ,nloc,kloc,q,N,x,y;
N := min(nops(r),nops(s)) ;
P := Array(0..N,0..N) ;
for kloc from 0 to N do
P[0,kloc] := 1 ;
end do:
for nloc from 1 to min(N,n) do
for kloc from 0 to N-nloc-1 do
q := x*op(1+kloc,r)+y*op(1+kloc,s) ;
if kloc > 0 then
P[nloc,kloc] := expand(P[nloc,kloc-1]+q*P[nloc-1,kloc+1]) ;
else
P[nloc,kloc] := expand(q*P[nloc-1,kloc+1]) ;
end if;
end do:
end do:
coeff(coeff(P[n,0],x,n-k),y,k) ;
end proc:
DELTA := proc(r::list,s::list)
local P ,nloc,kloc,q,N,x,y,k;
N := min(nops(r),nops(s)) ;
P := Array(0..N,0..N) ;
for kloc from 0 to N do
P[0,kloc] := 1 ;
end do:
printf("1,\n") ;
for nloc from 1 to N do
for kloc from 0 to N-nloc-1 do
q := x*op(1+kloc,r)+y*op(1+kloc,s) ;
if kloc > 0 then
P[nloc,kloc] := expand(P[nloc,kloc-1]+q*P[nloc-1,kloc+1]) ;
else
P[nloc,kloc] := expand(q*P[nloc-1,kloc+1]) ;
end if;
if kloc = 0 then
for k from 0 to nloc do
printf("%d,",coeff(coeff(P[nloc,0],x,nloc-k),y,k) );
end do;
end if;
end do:
printf("\n") ;
end do:
return ;
end proc:
# @param r list of r
# @param s list of s
# @param delfcol If true delete first column of resulting array
# @param delfrow If true delete first row of resulting array
DELTAgf := proc(r::list,s::list,x,y,delfcol::boolean,delfrow::boolean)
local N,n,g;
N := min(nops(r),nops(s)) ;
g := 0 ;
for n from N to 1 by -1 do
(op(n,r)*x+op(n,s)*x*y)/(1-g) ;
g := factor(%) ;
end do:
g := 1/(1-g) ;
if delfcol then
g := (g-subs(y=0,g))/y;
end if;
if delfrow then
g := (g-subs(x=0,g))/x;
end if;
g := factor(g) ;
printf("\nG.f.: (%a)/(%a). - ~~~~\n",factor(numer(g)),factor(denom(g))) ;
g ;
end proc:
# Count the frequency of n in the list L
# @param L The list of elements.
# @param n The needle to search for.
# @return The number of occurences of n in L.
# @since 2016-03-06
freq := proc(L::list,n::integer)
local ct,l;
ct := 0 ;
for l in L do
if l = n then
ct := ct+1 ;
end if;
end do:
return ct;
end proc:
# The number of picking c parts of an existing set of C
# elements (optionally picking some elements more than once),
# without regard to the ordering of the elements of C.
pickP := proc(C::integer,c::integer)
binomial(C+c-1,c) ;
end proc:
# A multiset decompositon of L.
# @param L L[1] corresponds to the n=0 elements in the first column
# @param n row index >=1 into L, to element n of the root sequence.
# The elemen at index 0 of the root sequence is not in that list.
# @param f column index >=1 and <=n. If f=1, pick L[f] ;
# @example L=A000081 starting at n=1 gives A033185.
# @example L=A038055 starting at n=1 gives A271878.
# @example L=A000079 starting at n=1 gives A209406.
# @example L=A000055 starting at n=1 gives A095133.
# @example L=A001349 starting at n=1 gives A201922.
# @example L=A038059 starting at n=1 gives A271879.
# @example L=A000012 starting at n=1 gives A008284.
# @example L=A000027 starting at n=1 gives A089353.
# @see arxiv.org:1603.00077
MultiSet := proc(L::list,n::integer,f::integer)
option remember ;
local c,ct,i,cof,thisp,oldp,thisC,fr;
if n < 0 then
return 0;
elif n = 0 then
return 1;
elif f = 1 then
# refers to element n-1 (but maple index n)
return op(n,L) ;
else
ct := 0 ;
# cannot use combinat[partition](N,f) here because
# f would not mean the number of parts...
for c in combinat[partition](n) do
if nops(c) = f then
oldp := -1 ;
# print("N=",n,"f=",f,"c=",c) ;
cof := 1;
for i from 1 to nops(c) do
thisp := op(i,c) ;
thisC := procname(L,thisp,1) ;
fr := freq(c,thisp) ;
if thisp <> oldp then
cof := cof*pickP(thisC,fr) ;
# print("N=",n,"f=",f,"c=",c,"factor",pickP(thisC,fr),"from fr",fr) ;
oldp := thisp ;
end if;
end do:
ct := ct+cof ;
end if;
end do:
ct ;
end if;
end proc:
MultiSetSum := proc(L::list,n::integer)
add( MultiSet(L,n,f),f=1..n) ;
end proc:
# A labeled multiset decompositon of L.
# Luschny's Bell transform (apart from the offset)
# @param L L[1] corresponds to the n=0 elements in the first column
# @param n row index >=1 into L, to element n of the root sequence.
# The elemen at index 0 of the root sequence is not in that list.
# @param f column index >=1 and <=n. If f=1, pick L[f] ;
# Examples: A143543 generated by A001187
# Examples: A307804 generated by A123543
# Examples: A307806 generated by A003515
# Examples: A228550 generated by A033678
# Examples: A105599 generated by A000272
# @since 2019-04-29
LblSet := proc(L::list,n::integer,f::integer)
option remember ;
local c,ct,i,cof ;
if n < 0 then
return 0;
elif n = 0 then
return 1;
elif f = 1 then
# refers to element n-1 (but maple index n)
return op(n,L) ;
elif f > n then
return 0 ;
else
ct := 0 ;
# i.e. partitions of n into f positive parts
for c in combinat[composition](n,f) do
cof := combinat[multinomial](n,op(c)) ;
for i from 1 to nops(c) do
cof := cof*procname(L,op(i,c),1) ;
end do:
# print("n=",n,"f=",f,"parts=",c,"multi=",cof) ;
ct := ct+cof ;
end do:
ct/f! ;
end if;
end proc:
# @param L The list a(i), i>=1, of the number of labelled structures with 1 part and i nodes
# @param n The total number of nodes
# @return the number of not necessarily connected labelled structures with n nodes
# Examples: A000012 creates A000110
# Examples: A004109 creates A002829.
# Examples: A001710 (with A001710(1)=0) generates A001205.
# Examples: A001710 starting 1,0,1,3,12... (offset 1 here) generates A108246
# Examples: A003515 generates A003514
# Examples: A059166 generates A059167
# Examples: A272905 generates A005815
# Examples: A092430 generates A007833
# Examples: A096332 generates A066537
# Examples: A079306 generates A006024
# Examples: A053941 generates A047656
# Examples: A000272 generates A001858
# @since 2019-04-29
LblSetSum := proc(L::list,n::integer)
add( LblSet(L,n,f),f=1..n) ;
end proc:
# Inverse function of LblSetSum()
# @param L The number of lableled structures with n>=1 nodes
# @param n The index >=1 into the list of labeled structures with 1 component
# @param estim A pair of minimum andmaximum estimators for LblSetsumi()
# @return A list with a lower and upper estimate of LblSetSumi()
# @since 2019-05-05
LblSetSumiAux := proc(L::list,n::integer,estim::list)
option remember;
local a,Ltst ;
if n = 0 then
return 1 ;
elif n = 1 then
return op(n,L) ;
else
# trial and error check of the n'th component of the inverse
a := floor((op(1,estim)+op(2,estim))/2) ;
[seq(LblSetSumi(L,i),i=1..n-1),a] ;
Ltst := LblSetSum(%,n) ;
if op(1,estim) = op(2,estim) then
if Ltst = op(n,L) then
return a ;
else
error
end if;
else
if Ltst = op(n,L) then
return a ;
elif Ltst > op(n,L) then
return procname(L,n,[op(1,estim),a]) ;
else
return procname(L,n,[a,op(2,estim)]) ;
end if;
end if;
end if;
end proc:
# Inverse function of LblSetSum()
# @param L The number of lableled structures with n>=1 nodes
# @param n The index >=1 into the list of labeled structures with 1 component
# @return The number of labeled structures with 1 component and n nodes.
# @since 2019-04-29
LblSetSumi := proc(L::list,n::integer)
option remember;
local a,Ltst ;
if n = 0 then
return 1 ;
elif n = 1 then
return op(n,L) ;
else
# trial and error check of the n'th component of the inverse
# assuming that all values are non-negative in the original list
return LblSetSumiAux(L,n,[0,op(n,L)]) ;
end if;
end proc:
# Discriminator sequence derived from a strictly increasing sequence.
# @param L The input sequence
# @param checkSrtd If true, an additional check is run to ensure the input sequence
# is strictly monotonic.
# @return The initial elements of the discriminator.
# If the prameter L is not an increasing sequence, that output is the empty list.
# @since 2016-05-04
# @author R. J. Mathar
# @example with input of squares, generating A016726
# L := [seq(n^2,n=1..10)] ;
# DISCRTR(L,true) ;
#
# @example with input the cubes, generating A192429
# L := [seq(n^3,n=1..10)] ;
# DISCRTR(L,true) ;
DISCRTR := proc(L::list, checkSrtd::boolean)
local Ldisc, i ,n ,m, el1, mworks, modset;
# initial discriminator list
Ldisc := [] ;
# number of elements in input list
n := nops(L) ;
# Test that input is monotonic
if checkSrtd then
for i from 2 to n do
if op(i,L) <= op(i-1,L) then
return [] ;
end if;
end do:
end if:
# compute the i'th element of the output list by considering
# the first i elements of the input
for i from 1 to n do
# brute force examination of the moduli
for m from 1 do
# examine all ordered pairs of input sequence
[seq(modp(op(el1,L),m) ,el1=1..i)] ;
modset := convert(%,set) ;
if nops(modset) = i then
# found an m that works: append to output
Ldisc := [op(Ldisc),m] ;
break ;
end if;
end do:
end do:
return Ldisc ;
end proc:
# given a squence with g(x) find f(x) such that g(x)=( f(x^2)+f(x)^2)/2, according
# to the cycle index of the symmetric group with 2 elements .
# @since 2016-05-04
# @author R. J. Mathar
# @param L The input sequence with leading value 1,
# representings g(x) = L[1] + x*L[2] + x^2*L[3]+..
# @return the truncated generating functin f
# @examples This relates A005995 with A000217 (without leading zero), or A005993 with A000027 (offset 0)
INVRTSYM2 := proc(L)
local g,n, f,itr;
n := nops(L) ;
g := gfun[listtoseries](L,x) ;
f := 1;
for itr from 1 to n do
f := sqrt(2*g-subs(x=x^2,f)) ;
f := taylor(f,x=0,itr+2) ;
f := convert(f,polynom) ;
end do:
f ;
end proc:
# Given sequences L1 and L2, compute the exponential convolution
# @param L1 First multiplicative squence, starting at a(1)
# @param L2 Second multiplicative sequence, starting at b(1)
# @return The multiplicative sequence defined by Lelechenko in arxiv:1405.7597
CONVE := proc(L1,L2)
local L,n,a,pe,i,p,e,f,d ;
L := [1] ;
for n from 2 do
a := 1 ;
for pe in ifactors(n)[2] do
p := pe[1] ;
e := pe[2] ;
f := 0 ;
for d in numtheory[divisors](e) do
if p^d > nops(L1) or p^(e/d) > nops(L2) then
return L ;
end if ;
f := f+L1[p^d]*L2[p^(e/d)] ;
end do:
a := a*f ;
end do:
L := [op(L),a] ;
end do:
end proc:
# Given an ordered sequence L=[a1,a2,a3,..,ai] and g.f. (1+x^a1)*(1+x^a2)*(1+x^a3)..
# provide the sequence with that g.f.
# Which if a1, a2,a3 are parts in a partition, provide the number of
# ways n can be partitioned in distinct parts taken from the ai.
# @example: if L=[1,2,3,4,6..] = A029744, the result will be 1,1,1,2,2,2,... = A008620
# @example: if L=[1,2,3,4,5..] = A000027, the result will be 1,1,1,2,2,3,... = A000009
# @example: if L=[1,2,3,5,6] = A005117, the result will be 1,1,1,2,1,2,... = A087188
# @example: if L=[1,1,2,2,3,3,4,4..] = A008619, the result will be 1,2,3,6,9,.. = A022567
# @example: if L=[2,3,5,7,11,13] = A000040, the result will be 1,0,1,1,0 = A000586
# @example: if L=[2,3,5,7,11,13] = A109763, the result will be 1,0,2,2,.. = A280245
# @example: if L=[1,4,6,8,10,12] = A053012, the result will be 1,1,0,0.. = A226749
# @example: if L=[1,2,6,24,120] = A000142 without reps , the result will be 1,1,1,1,0,0 = A115944
# @example: if L=[4,6,9,10,14,..] = A001358 without reps , the result will be A112020
# @since 2023-03-01
# @author R. J. Mathar
DISPART := proc(L)
local x ;
mul( 1+x^p,p=L) ;
taylor(%,x=0,max(op(L))) ;
gfun[seriestolist](%) ;
end proc:
# Input represents list L=[1,a1,a2,a4,..] which means leading coefficient must be 1.
# Output contains a list of (possibly repeated) exponents b1, b2, b3 such
# that 1+a1*x+a2*x^2+a3*x^3+... = (1_x^b1)*(1+x^b2)*(1+x^3b)*..
# example: L= [seq(i,i=1..20)] = A000027 produces [1,1,2,2,4,4,8,8,16,16,...]
# which represents (1+x)^2*(1+x^2)^2*(1+x^4)^2*(1+x^8)^2*... =((1+x)*(1+x^2)*(1+x^3)*..)^2 and
# expresses the fact A000027 is the autoconvolution of A000012.
# @since 2023-03-01
# @author R. J. Mathar
DISPARTi := proc(L)
local Lout,g,x,hd,ld ;
Lout := [] ;
g := add(op(i,L)*x^(i-1),i=1..nops(L)) ;
while g<> 1 do
hd := degree(g) ;
ld := ldegree(g-1) ;
g := g/(1+x^ld) ;
Lout := [op(Lout),ld] ;
if ld = 0 or nops(Lout) > 200 then
break ;
end if;
g := taylor(g,x=0,hd+1) ;
g := convert(g,polynom) ;
end do:
Lout ;
end proc:
# Support function that sorts 2 lists with sort() and
# considers a list "prior" to another list if its first
# entry is numerically smaller. This is the kind of sorting values
# and i-indices in the signature sequences.
# @param l1 a list of the form [val,...] where the dots mean
# the actual number of elements is unimportant and only the first
# item is used.
# @param l2 a list of the form [val,...] where the dots mean
# the actual number of elements is unimportant and only the first
# item is used.
# @return true if the val of l1 is smaller than the val of l2.
# @since 2024-05-28
# @author R. J. Mathar
SIGSEQsort := proc(l1::list,l2::list)
if op(1,l1) < op(1,l2) then
return true ;
else
return false ;
end if
end proc:
# Construct the signature sequnce of the constant x.
# @param x The constant defining a list of i+j*x values.
# @param vmax The maximum value of i+j*x considered for
# evaluation. The length of the
# signature list that will be returned grows with vmax.
# @return the list of i-values after the values of i+j*x, i,j>=1
# i+j*,x < vmax, are sorted by the standard increasing arithmetic value.
# @example
# For the sign. seq. of log(2), A373212 run
# Digits := 100 ; SIGSEQ(log(2),50) ;
# For the sign. seq. of Pi^2, A373211, run
# Digits := 100 ; SIGSEQ(Pi^2,50) ;
# @since 2024-05-28
# @author R. J. Mathar
SIGSEQ := proc(x,vmax)
local TBsrtd,i,j ;
# TBsrtd contains a list of [[i1+j1*x,i1],[i2+j2*x,i2],...]
# i.e., pairs of i+j*x and their i-indices.
TBsrtd := [] ;
# double loop over i and j to gather the i+j*x
for i from 1 do
if i > vmax then
break ;
end if;
for j from 1 do
if evalf(i+j*x) > vmax then
break ;
end if;
TBsrtd := [op(TBsrtd),[evalf(i+j*x),i]] ;
end do:
end do:
# sort the values of i+j*x in natural increasing order
sort(TBsrtd,SIGSEQsort) ;
# extract the i-values of the sorted list
[seq(op(2,v),v=%)] ;
end proc:
digcatL2 := proc(L1,L2)
digcatL([op(L1),op(L2)]) ;
end proc:
