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:
