###################################################################### ##SVT.txt: Save this file as SVT.txt # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read `SVT.txt` # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: April 26, 2016 #This version: May 13, 2016 with(linalg): print(`Created: April 26, 2016`): print(`This version: May 13, 2016`): print(` This is SVT.txt `): print(`It accompanies the article `): print(` On the Sequences Enumerating Singular Vector Tuples of d-dimensional n by n ... by n tensors `): print(`by Shalosh B. Ekhad and Doron Zeilberger `): print(`and also available from Zeilberger's website`): print(``): print(`Please report bugs to zeilberg at math dot rutgers dot edu`): print(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg/ .`): print(`---------------------------------------`): print(`For a list of the Supporting procedures type ezra1();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`---------------------------------------`): print(`---------------------------------------`): print(`For a list of the MAIN procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`---------------------------------------`): with(combinat): Digits:=20: ezra1:=proc() if args=NULL then print(` The supporting procedures are: Asy, AsyC, A3f, B3f, A4f,B4f, CharRoots, Cnk, EmpirAsyC, findrec, Findrec, MMMT,`): print(` SVTseqSlow, SeqFromRec, Zinn `): print(``): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: Adf, Bdf, Bdr, SVT, SVTseq, SVTseq3, SVTseq4, SVTseqPC, SVTpaper `): print(` `): elif nops([args])=1 and op(1,[args])=A3f then print(`A3f(a1,a2,a3): the sum of B3(b1,b2,b3) for 0<=b1<=a1, 0<=b2<=a2, 0<=b3<=a3 `): print(`Try:`): print(` A3f(3,3,3); `): elif nops([args])=1 and op(1,[args])=A4f then print(`A4f(a1,a2,a3,a4): the sum of B4(b1,b2,b3,b4) for 0<=b1<=a1, 0<=b2<=a2, 0<=b3<=a3, 0<=b4<=a4 `): print(`Try:`): print(` A4f(3,3,3,3); `): elif nops([args])=1 and op(1,[args])=Adf then print(`Adf(a): c(a), where a is a list of non-negative integers, of length k, say, outputs the coefficient of x1^a[1]*...xd^a[d]`): print(`in the rational function 1/(1-e_2(x1, ..., xd)- 2*e_3(x1, ..., xd)- (d-1)e_d(x1, ..., xd))*(1/((1-x1)*...*(1-xd))`): print(`Try:`): print(` Adf([3,3,3]); `): elif nops([args])=1 and op(1,[args])=Asy then print(`Asy(ope,n,N,mu,M): inputs a linear recurrence operator with polynomial coefficients, `): print(`in the variable n, and where the forward shift operator is denoted by N , i.e. Nf(n):=f(n+1)`): print(`and one of its characteristic roots, mu, and a positive integer M. finds the asymptotic series`): print(`annilihating the recurrence, to order M `): print(`regardless of initial conditions`): print(`Try:`): print(`Asy(n^2*(N-1)*(N-2)*(N-3)+ n*(N^2+5*N+1)+ (N^2+7*N+1),n,N,3,4);`): elif nops([args])=1 and op(1,[args])=AsyC then print(`AsyC(ope,n,N,mu,M,INI,K): inputs a linear recurrence operator with polynomial coefficients, `): print(`in the variable n, and where the forward shift operator is denoted by N , i.e. Nf(n):=f(n+1)`): print(`and one of its characteristic roots, mu, and a positive integer M, as well as a list`): print(`of numbers, INI, whose length is the order of the recurrence (i.e., degree(ope,N)), `): print(`It also inputs a large integer K telling us how far to go in order to estimate/guess the constant in front`): print(`It outputs the asymptotic series, INCLUDING THE CONSTANT IN FRONT`): print(`annilihating the recurrence, to order M, and satisfying the given initial conditions `): print(`It returns FAIL if it is the wrong mu`). print(`Everything is rigorous except for the constant, that is estimated. `): print(`Try:`): print(`AsyC(n^2*(N-1)*(N-2)*(N-3)+ n*(N^2+5*N+1)+ (N^2+7*N+1),n,N,3,4,[1,1,1],1000);`): elif nops([args])=1 and op(1,[args])=B3f then print(`B3f(a1,a2,a3): the coefficient of x1^a1*x2^a2*x3^a3 in 1/(1-(x1*x2+x1*x3+x2*x3)-2*x1*x2*x3)`): print(`Try:`): print(` B3f(3,4,6); `): elif nops([args])=1 and op(1,[args])=B4f then print(`B4f(a1,a2,a3,a4): the coefficient of x1^a1*x2^a2*x3^a3*x4^a4 in 1/(1-e2(x1,x2,x3,x4)-2*e3(x1,x2,x3,x4)-3*e4(x1,x2,x3,x4))`): print(`Try:`): print(` B4f(3,3,3,3); `): elif nops([args])=1 and op(1,[args])=Bdf then print(`Bdf(a): f(a), where a is a list of non-negative integers, of length d, say, outputs the coefficient of x1^a[1]*...xd^a[d]`): print(`in the rational function 1/(1-e_2(x1, ..., xd)- 2*e_3(x1, ..., xd)- (d-1)e_d(x1, ..., xd))`): print(` Try:`): print(` Bdf([3,3,3]); `): elif nops([args])=1 and op(1,[args])=Bdr then print(`Bdr(a): f(a), where a is a list of non-negative integers, of length d, say, and a positive integer, r, outputs the coefficient of x1^a[1]*...xd^a[d]`): print(`in the rational function 1/(1-e_2(x1, ..., xd)- 2*e_3(x1, ..., xd)- (d-1)e_d(x1, ..., xd))^r`): print(` Try:`): print(` Bdr([3,3,3],2); `): elif nops([args])=1 and op(1,[args])=CharRoots then print(`(ope,n,N): inputs a linear recurrence operator with polynomial coefficients, ope`): print(` in the variable n, and where the forward shift operator is denoted by N , i.e. Nf(n):=f(n+1) `): print(` and outputs the ordered list`): print(`of its characteristic roots. If it has complex roots, it returns FAIL. Try:a`): print(`CharRoots((n+4)*N^2-n*N-1,n,N);`): print(``): elif nargs=1 and args[1]=Cnk then print(`Cnk(n,k): all the 0-1 vectors of length n with k 1's`): elif nargs=1 and args[1]=EmpirAsyC then print(`EmpirAsyC(L1,n,mu,t,M): inputs a list L1 of positive numbers believed to gave asympotics of the form`): print(`mu^n*n^t, finds the constant using asympotics of ofer M. Try:`): print(`EmpirAsyC([seq(5*3^i/i*(1+4/i),i=1..40)],n,3,-1,5);`): elif nargs=1 and args[1]=findrec then print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`): print(`the sequence f of degree DEGREE and order ORDER.`): print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`): elif nargs=1 and args[1]=Findrec then print(`Findrec(f,n,N): Given a list f tries to find a linear recurrence equation with`): print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `): print(`annihilating the seqeunce f. If it fails it returns FAIL. This means that either there is no recurrence`): print(`or its complexity is such that you need more terms. It does not handle solutions of recurrences with constant coefficients`): print(`e.g. try Findrec([seq(i,i=1..20)],n,N);`): elif nargs=1 and args[1]=MMMT then print(`MMMT(A,x): the generating function given by MacMahon's master theorem in the list of variables x.`): print(`Try:`): print(`MMMT([[0,1,1],[1,0,1],[1,1,0]],[x1,x2,x3]); `): elif nargs=1 and args[1]=SeqFromRec then print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`): print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`): print(`satisfied by the recurrence ope(n,N)f(n)=0`): print(`extends it to the first K values`): print(`For example, try:`): print(`SeqFromRec(N-n-1,n,N,[1],10);`): elif nops([args])=1 and op(1,[args])=SVT then print(`SVT(N):Inputs a list of non-neg. integers N=[n1,..., nd], and outputs the number of singular vector tuples `): print(`for a general n1 x n2 x ... x nd tensor.`): print(`it uses the constant-term formula given in the article`): print(`Shmuel Friedland and Giorgio Ottaviani, The number of singular vector tuples and uniqueness of best rank-one approximation of tensors,`): print(` Found. Comput. Math. 14 (2014), no. 6, 1209-1242.`): print(`Try: SVT([3,3,3]); `): elif nops([args])=1 and op(1,[args])=SVTpaper then print(`SVTpaper(d,n,Kama,M,K): inputs an integer d, a symbol n, a positive integer Kama, a positive integer M, and a large positive integer K`): print(`outputs an article, that tells you first the first Kama terms of the sequence enumerating`): print(`the number of singular vector tuples of an n by n by ... by n (d times) tensor, using the constant-term formula in the article:`): print(`Shmuel Friedland and Giorgio Ottaviani, The number of singular vector tuples and uniqueness of best rank-one approximation of tensors`): print(` Found. Comput. Math. 14 (2014), no. 6, 1209-1242.`): print(` Then it tries to find a homogeneous linear recurrence equation satisfies by that sequence using that many terms. If it fails `): print(`the paper is finished there. If it succeeds, it states it, and tries to use it to find an asympotic formula.`): print(`if it succeeds it tells you so. Try:`): print(`SVTpaper(3,n,70,5,2000):`): elif nops([args])=1 and op(1,[args])=SVTseq then print(`SVTseq(d,N): inputs a positive integer d, and a positive integer N, outputs the first N terms of SVT([n$d]) from n=1 to n=N.`): print(`Warning, very slow, done directly from the definition. If d=3 do SVTseq3(N), If d=4, do SVTseq4(N),`): print(`Try: `): print(`SVTseq(3,30); `): elif nops([args])=1 and op(1,[args])=SVTseqSlow then print(`SVTseqSlow(d,N): inputs a positive integer d, and a positive integer N, outputs the first N terms of SVT([n$d]) from n=1 to n=N.`): print(`Warning, very slow, done directly from the definition. If d=3 do SVTseq3(N), If d=4, do SVTseq4(N),`): print(`Try: `): print(`SVTseqSlow(3,30); `): elif nops([args])=1 and op(1,[args])=SVTseq3 then print(`SVTseq3(N): inputs a positive integer N, outputs the first N terms of SVT([n,n,n]) from n=1 to n=N.`): print(`Try: `): print(`SVTseq3(30); `): elif nops([args])=1 and op(1,[args])=SVTseq4 then print(`SVTseq4(N): inputs a positive integer N, outputs the first N terms of SVT([n,n,n,n]) from n=1 to n=N.`): print(`Try: `): print(`SVTseq4(6); `): elif nops([args])=1 and op(1,[args])=SVTseqPC then print(`SVTseqPC(d,N): like SVTseq(d,N), but for d=3 and N<=70 pre-computed to save time.`): print(`if not d=3 and N<=70 it uses SVTseq(d,N) `): print(`Try: `): print(`SVTseqPC(3,30); `): elif nops([args])=1 and op(1,[args])=Zinn then print(` Zinn(L): inputs a list of integers, uses Zinn-Justin's method to stimate`): print(` the mu, and theta such that`): print(`resh[i] is appx. Const*mu^i*i^theta`): print(`For example, try:`): print(`Zinn([seq(5*i*2^i,i=1..30)]);`): else print(`There is no ezra for`,args): fi: end: ######start from FindRec sn:=proc(resh,n1): -1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2): end: #Zinn(resh): Zinn-Justin's method to estimate #the C,mu, and theta such that #resh[i] is appx. Const*mu^i*i^theta #For example, try: #Zinn([seq(5*i*2^i,i=1..30)]); Zinn:=proc(resh) local s1,s2,theta,mu,n1,i: if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then RETURN(FAIL): fi: n1:=nops(resh)-1: s1:=sn(resh,n1): s2:=sn(resh,n1-1): theta:=evalf(2*(s1+s2)/(s1-s2)^2): mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))): [theta,mu]: end: #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrecVerbose:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then print(`There is some slack, there are `, nops(ope)): print(ope): RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrec:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: option remember: if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then RETURN(FAIL): fi: if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then RETURN(FAIL): fi: if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then RETURN(FAIL): fi: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: #FindrecOld(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with #poly coffs. #of maximum DEGREE+ORDER<=MaxC #e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2); FindrecOld:=proc(f,n,N,MaxC) local DEGREE, ORDER,ope,L: for L from 0 to MaxC do for ORDER from 0 to L do DEGREE:=L-ORDER: if (2+DEGREE)*(1+ORDER)+4>=nops(f) then print(`Insufficient data for degree`, DEGREE, `and order `,ORDER): RETURN(FAIL): fi: ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: #Findrec(f,n,N): Given a list f tries to find a linear recurrence equation with #poly coeffs. Findrec:=proc(f,n,N) local DEGREE, ORDER,ope: option remember: for ORDER from 1 while 2*(1+ORDER)+5+ORDER<=nops(f) do for DEGREE from 1 while (1+DEGREE)*(1+ORDER)+5+ORDER<=nops(f) do ope:=findrec(f,DEGREE,ORDER,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: #SeqFromRec(ope,n,N,Ini,K): Given the first L-1 #terms of the sequence Ini=[f(1), ..., f(L-1)] #satisfied by the recurrence ope(n,N)f(n)=0 #extends it to the first K values SeqFromRec:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: #end Findrec with(linalg): #findrecEx(f,DEGREE,ORDER,m1): Explores whether thre #is a good chance that there is a recurrence of degree DEGREE #and order ORDER, using the prime m1 #For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003); findrecEx:=proc(f,DEGREE,ORDER,m1) local ope,var,eq,i,j,n0,eq1,a,A1, D1,E1,Eq,Var,f1,n,N: option remember: f1:=f mod m1: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1: od: eq:= eq union {eq1}: od: Eq:= convert(eq,list): Var:= convert(var,list): D1:=nops(Var): E1:=nops(Eq): if E10 then RETURN(false): fi: if E1-D1>=1 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+1],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: if E1-D1>=2 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+2],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: true: end: ###end from FindRec #CharRoots(ope,n,N): inputs a linear recurrence operator with polynomial coefficients, #in the variable n, and where the forward shift operator is denoted by N , i.e. Nf(n):=f(n+1) #and outputs the ordered list #of its characteristic roots. Try: #CharRoots((n+4)*N^2-n*N-1,n,N); CharRoots:=proc(ope,n,N) local ope1,lu: ope1:=lcoeff(expand(ope),n): if degree(ope1,N)=0 then RETURN(FAIL): fi: lu:=[solve(ope1,N)]: lu:=sort(lu): end: #Asy(ope,n,N,mu,M): inputs a linear recurrence operator with polynomial coefficients, #in the variable n, and where the forward shift operator is denoted by N , i.e. Nf(n):=f(n+1) #and one of its characteristic roots, mu, finds the #an asymptotic series whose leading term is mu^n that satisfies the recurrence, regardless of initial conditions. Try: #Asy((n+5)*N^2-(3*n+4)*N+2*n+2,n,N,2); Asy:=proc(ope,n,N,mu,M) local Ope,ope1,gu,t, i,x,d,eq,var,c,f: Ope:=expand(numer(ope)): d:=degree(Ope,n): ope1:=coeff(expand(Ope),n,d): if simplify(subs(N=mu,ope1))<>0 then RETURN(FAIL): fi: gu:=add( x^d*subs(n=1/x,coeff(Ope,N,i)) * mu^i*(1+i*x)^t , i=ldegree(Ope,N)..degree(Ope,N)): gu:=expand(gu): gu:=taylor(gu,x=0,2): if expand(coeff(gu,x,1))=0 then RETURN(FAIL): fi: t:=solve(coeff(gu,x,1),t): f:=1+add(c[i]/n^i,i=1..M): var:={seq(c[i],i=1..M)}: gu:=add( x^d*subs(n=1/x,coeff(Ope,N,i)) * mu^i*(1+i*x)^t*subs(n=(1+i*x)/x,f), i=ldegree(Ope,N)..degree(Ope,N)): gu:=expand(gu): gu:=taylor(gu,x=0,M+3): eq:={seq(expand(coeff(gu,x,i)),i=2..M+1)}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: mu^n*n^t*subs(var,f): end: #AsyC(ope,n,N,mu,M,INI,K): See ezra AsyC:=proc(ope,n,N,mu,M,INI,K) local Ope,gu,L1,C: Ope:=expand(numer(ope)): gu:=Asy(ope,n,N,mu,M): if gu=FAIL then RETURN(FAIL): fi: if nops(INI)<>degree(Ope,N) then print(INI, `should be of length`, degree(Ope,N)): RETURN(FAIL): fi: L1:=SeqFromRec(Ope,n,N,INI,K): if abs(evalf(L1[K]/subs(n=K,gu)))<10^(-10) then RETURN(FAIL): fi: C:=evalf(L1[K]/subs(n=K,gu)): C:=identify(evalf(C,10)): C*gu: end: #SVT(N):Inputs a list of non-neg. integers N=[n1,..., nd], and outputs the number of singular vector tuples #for a general n1 x n2 x ... x nd tensor. #Warning, very slow, done directly from the definition. If N is 3-dim do SVT3 SVT:=proc(N) local z,d,i,gu,lu: option remember: d:=nops(N): lu:=add(z[i],i=1..d): gu:=1: for i from 1 to d do gu:=gu*normal( ((lu-z[i])^N[i]-z[i]^N[i])/(lu-2*z[i])): od: gu:=expand(gu): for i from 1 to d do gu:=coeff(gu,z[i],N[i]-1): od: gu: end: #SVTseqSlow(d,N): inputs a positive integer d, and a positive integer N, outputs the first N terms of B([n$d]) from n=1 to n=N. #Warning, very slow, done directly from the definition. If N is 3-dim do SVTseq3(N),If N is 4-dim do SVTseq4(N), #Try: #SVTseqSlow(3,30); SVTseqSlow:=proc(d,N) local n1: [seq(SVT([n1$d]),n1=1..N)]: end: #SVTseq(d,N): inputs a positive integer d, and a positive integer N, outputs the first N terms of B([n$d]) from n=1 to n=N. ##Try: #SVTseq(3,30); SVTseq:=proc(d,N) local n1: [seq(Adf([n1$d]),n1=0..N-1)]: end: #SVTseqPC(d,N): Like SVTseq(d,N) but if d=3 and N<=70 then it gives pre-computed values #Try: #SVTseqPC(3,30); SVTseqPC:=proc(d,N) local n1,L1,L2: L1:= [1, 6, 37, 240, 1621, 11256, 79717, 572928, 4164841, 30553116, 225817021, 1679454816, 12556853401, 94313192616, 711189994357, 5381592930816, 40848410792017, 310909645663332, 2372280474687277, 18141232682656320, 139010366280363601, 1067160872528170536, 8206301850166625797, 63203453697218605440, 487480961825645988721, 3764885111909676266856, 29112650384540312124397, 225377607037294291831008, 1746648730442554074446041, 13549903631370307542264936, 105214484529950220234997141, 817709228332592281063815168, 6360391464220968384566019745, 49511928791684034124645730676, 385707501237368849206626328717, 3006846684379174364465426255616, 23455977069620091841962275705089, 183092319857177553953377636176456, 1430037723632571918813195121704997, 11175656244135989347042802862284928, 87384630486487451749353666860367841, 683632453962320044584305182576200288, 5350886598153793514481537897337176493, 41902005901674406571033924405839308768, 328277426887267454292261947436072998521 , 2572978040988090019983161573900671390056, 20174938455027696505833506309037744673717, 158256711967850325867904932175164613805568, 1241878121565109721305598151979300621187041, 9748914697534511954263668793857706936885416, 76557533012847845401553482298831558586858781, 601406624861928916400901634219448551509045216, 4725977743774598911951207466985442875802004041, 37149518272412304070221079042862933387873130216, 292110984170941385701186299397195372784867345221, 2297582777856431871486896289498367921392694659456, 18076679242941480643112276737449814240450399614961, 142261109059091890572028251448021584831777493283976, 1119875738825722199969846898012482760095108120469421, 8817913547145642408782946965788762982054339519598816, 69449739062865201034991808469623104398486763332361881, 547117715028973597496117057639950244728453193275378536, 4311149329690424172211655089600627818487132234309850581, 33978514528762393389117490950553726381113420890825711616, 267862418462871260646089896898133632801592296425890350401, 2112090051496950649407979640746003519860385829377421719956, 16657262165394749047223410973405780093084150560578941606221, 131396189046352933535318043187955359702871390804208939532416, 1036686394280916864442473886582846818530956264890966401359201, 8180787959749319986573210152825554401372254551057541041233416]: L2:= [1, 24, 997, 51264, 2940841, 180296088, 11559133741, 765337680384, 51921457661905, 3590122671128664, 252070718210663749, 17922684123178825536, 1287832671004683373753, 93368940577497932331288, 6821632357294515590873917, 501741975445243527381995520, 37121266623211130111114816929, 2760712710223967190110979892824, 206267049696409355312012281872181, 15475371710788303634337687336412224, 1165407938272751946314416465834086889, 88061985234092428001001009986265842520, 6674841972181830152646239523454436828797, 507365393175531718935697292222219542209024, 38665843085507634007413709865882990263469041, 2953740030556494206297508806673140156327404248, 226139608318389776420643582181314917081622957061, 17348847768732971801911980996414434839471856232768, 1333495999756767180730069621847180744764824179044825, 102679492043788645453302773152678289049925015355480088, 7919467782979625176350608079488125284094134882129213277, 611762406707957516785060447056716898028654120149205549056, 47326297794171750329395108432285966368588072502073467990849, 3666210684055121798038265960682272734448162887732387739437656, 284375989192008571736288450556047894658675111432554006369198101, 22085002503311452557061359446412536526728124547277566177787892288, 1717125381944152148262069150947173852066288121156424979958627271625, 133653371349899710142461127502334621416505147105558621766272672907736, 10413736715423917925250315399202563233072956051319130146752533719981053, 812193134071724372030028233213580873621215660428760806777376397599783424, 63404056487163136402118230809098571181603513366964484046606874882278962769, 4954048702957761157113524005575838540485231482644154318297501579169017363736, 387410031814835435731535807396509767481143902823175117339431737579621537848213, 3032019104953631935927245935785181752925383265226972671661881411766405320238624\ 0, 2374803834313774545027493708891468517565412197712400933613898364537188007266\ 415065, 18614165945764264383509075035019318709521380225496240889482096963404716\ 0142280221656, 1460043170995964708856156528528938389222377695391471169729078084\ 4366152779284725947757, 1145990988900840641498974352772435574812691486580998832\ 888840442789761192955239298207744, 90007356802866260340681675743343295423352567\ 419796418449737235391358299924670403492251105, 70736693805083529397480221136055\ 92581392110691450778592637430119392203264791502783256731864, 556250953151679190\ 326199484492492324591698905227339332845928022594121900029904848181878033749, 43\ 7669196230078776788737597365736295916093100444453781362348341558871850881876505\ 01694445489216, 344556869088661514142383463598133274533590051098179797634568099\ 3602708821591585237841811565673065, 2713979783990190110602974150069380468418799\ 26951597708635047601874185952373678877541675978258404568, 213882183493150066670\ 6483254859097343585216990790743086036561160356374638779426179264517198748656322\ 9, 1686385703886231464111507816879306616123474300541011700087381363447966207871\ 114366608478451366840144384, 13302888122049843048583783811778374771716769379195\ 1549233571332531946972464580592254702500946111712214577, 1049867536127299365636\ 5497227516437477855862285648134225237865420533609968529304268060706932598043796\ 535256, 82892607577994651456907533699515307250159154916797481720621027970596587\ 8763162178927441457954079168327435141, 6547618275154417236807627493509617188339\ 0561330168239058377228182666146594166925063393822640826631362216652096, 5174054\ 6109669857407254864558010313843670271715836656387781610740275475557105777415128\ 17639228972097875063728153, 409027887749565473877643124644902175375544568079990\ 109075128262610563868460599234291869165223723204244596556580376, 32347708141131\ 6703815112633244396816051373975539285755764374214440614410184428935385384474359\ 17662519366143191410589, 255916030760866495394527698951579075886076340724451964\ 0658122483760011731366878142308998366812224871870032223390662656, 2025395240895\ 1603246698248194238532846559543259009791531066210961323946905292332059644843530\ 9068417769429529956502685313, 1603524387663098386262074342517422638718342567678\ 2167996633336356019580148925904327992960651893088013216529811030500811352, 1269\ 9605587989209055001298402406398177669502575479970411995455982570766601035104530\ 61245176163841999507880978168128227321301, 100611904723621504989197065102687600\ 5653465553050036319143930926830585636462050949830520738795678147410970775185526\ 58278518336, 797349578873080205697751714231745229204003541172569345165375061133\ 1139698481136957917921419063990733623143125922769418501460361, 6320979719089295\ 5360382815502363297770861920831751714025080216909873677318083716624321039453174\ 6401261452355825000537762507812568, 5012477014633297839861201240291319860394600\ 7439132120509864627446835189634880231488902214692421032853360665076770619392321\ 171964157, 39760246007434002263122771664553011899184436308132273433144559413030\ 99745331143603268930380461032113850457743160958373952262009784832, 315479292291\ 3956577203778130688262135870108448404747784377827959627670582691460350012918017\ 35333187359438736270047927662791312330794513, 250388500221019531312713421084697\ 4108947429561673585282648831960094871813676160195925602880976709018133779466713\ 5887587705734069647224728, 1987816823728556596123625241672421825758958175776974\ 8598555938762939620589203386825268137188450332215916917816807530042965920535983\ 12894517, 157853302741847454273439691565911926574130298546204629589813585802776\ 508147757855151240838895127044522108017538593152066002261339870960194880, 12538\ 4340752346336853109321776495015359570646919793393138881734833630652099581700596\ 32146740282752565168203907489021983261477646662277381157849, 996190396276612251\ 1705647985861203267338110511349611130284466638206704819483647150478457066892315\ 24090521025556311814462365392336365072956287320, 791677018490697051619705534671\ 0226992773793074176044840611548011015860153453335364163605644990813188457321020\ 0459832117005673118892374654134729741, 6292999305076221635289521103045702396556\ 7376258264950085614014090402431439242581649692447928776432672043049984187261251\ 24412095339186826385979822080, 500344014007587203400001211826809049233422877132\ 7641155252146099007741056540797243004321781335218618839157977192712385019667895\ 41245938813176876065441, 397904295110177421949295675871828872974954018166295543\ 9186408652489360565644285599782161186946670801542126458929974693277169748109675\ 4808754698102758680, 3165082579580981019251700832149775222812744145017531963735\ 1662364725334284986046134953645539434430463186131518819367838854558031474818175\ 64155533105479557, 251817353557598225180603802212268857541397484886337857197230\ 2279551986452480708827688768932916629056759166300223219344311793862722539742397\ 66085268872234560, 200391008125539684200079917324872930227709966655258060255824\ 6613861124071281344613410787763539498777520450453286440611775838083231201347796\ 8123653989200572169, 1594999756584197937091376922485656814469152689497999394432\ 8298668163227068670231581456699380322714150785801699035367785590624403429191411\ 25154923097178694101784, 126978664303800230048248365830511062083289620698244402\ 6029228239668793253643263609196501277414985310316356995046620271553280116065881\ 03981425934300726451012014445, 101108257767235959829205622808187626030087725574\ 0083953925780215620136332179178655408085422816846337203330708366098913467279944\ 8372382157581384048668873433076978176, 8052417864994305539040572420568131583499\ 9730642478434734308748734687047297220429693220672311182622622531466083639080822\ 1475211272627518645358168693593115314639650865, 6414280001915575133736154854176\ 2126269147211399946165883923705134557821649378820742943693615808568491307650353\ 408642911901481644815713508372697983466547381624580088984, 51103382899294247295\ 3254172102565644737684958492382612372875509431010005371969091516202908232303526\ 7467941530833349312464046507487025779488297520521107193830389633201749, 4072206\ 2143744200374032801101071020066037586589924393678952402729912329661888558484672\ 2264183193238403129380380290080150371748126599234960329602439950125606394240547\ 435328, 32455369529124983521698348415787175273222918527213635107991789844112393\ 2936474026634592819859895767430163620142373287545847142914231553647287608110560\ 38156251482959017915577, 258713088005783581345986148647833214410935763755425067\ 8489673329018435731802025930947928694299728822248249096458582782759970039327392\ 397322391481123153017756601550869712717016, 20626412580361576753969781335412891\ 1746799490826642473542333196922308931176801844910842778549891743970867718151171\ 563738206902856248717386966184856153669522720001199393298593101, 16447538968034\ 7639937459460201287562707660861565560865034009374275547406662591987035247125543\ 1890673538137993003726605995976897147118933198139273849479582463692850041026480\ 0493568, 1311742475245668886089440516090460879217389738410928738340599513010903\ 3702245389167597855634031489176850524751528156463035826886634819842524312077986\ 12394112462467087561914459972545, 104632169460506745543037354025742426410639372\ 5611297068124324998083522514289486856007612900419270319119779715818648862820481\ 94024504765585587462025050901069089290189485679471464178904, 834736628061662143\ 0072352616934982796497241260060485638513673536445597251986895370586080624266250\ 2316702764017941400119600532512199664174644505379064037808273162077272192460332\ 39561461, 666039438770326033747497326351168176537246337683607561298874676981936\ 3305591339385803585921408064482807876740189356283370522312105693971947064238639\ 81366352404470368586638132792642025024, 531514773440676360168130393354967776826\ 7528050086635907987786001741901310843264343107577234932454090569312661525257204\ 3833627467420865066902434459184172282414583999076165755373790928695049, 4242231\ 6052795413113172460739819368710758076902717876886789795582992846251243404121309\ 8386912180274694681833123444117101483385099407291468896319209536042224479872146\ 5346137039825333940805720, 3386380938240000283548844348656961586868383333620671\ 8053814194190770673931801651648103684148227782716806488695026791990480815721112\ 6302747981626037303254014002121929535185789724389130405980445, 2703575130703808\ 7514099317793796184145755757661734868175305795591792813036749399071326545031768\ 1600515636653074082234098163769402728927754852637057359871864909102466484033943\ 26786860528631069184, 215874379089302750363434099339956856165439763543601478254\ 8254081295799471159730098573193501379982496670402592422254595871368004176624119\ 623107904933101646550373409097151531725469605047871702334289, 17239419984865383\ 2578851474687227714626127688851198213618562023982474832941210578405158001998691\ 5507889321175240855981695826859982216933686899920181463978651507754984377988758\ 57094539817536031491224, 137689868153063928099873333533128517858367534535863772\ 7292604837460531139474545320078959898100562293194691352804447704120110227180638\ 5791978778226851094950539722271728058633084683603811083199796810933, 1099861575\ 2100832605395344771452142626486926763354545797087130977669518521307026996298176\ 5578888814716562889155614549732572477892130623640328748890689675678139417224122\ 8077994858349916559930497196726592, 8786779802522957580011353673943727787972540\ 6144613920377738046535236628937979712714677511935324724006177134550964707202769\ 0155528580474501859184149621594961451829759048495015594607481368258200517339691\ 77, 702063036987849353607209476891253472021486502536546242088832391600897387166\ 5558806008653685082823757772734359939779427221609358458166972237222212991391831\ 674267898210779250381540724993215430710809091637080, 56101715022433697024041069\ 7517664062038838213270295968405053561236911744923709706505124711432911080090546\ 8860981638271595028171871478531471802047757125291790458131645292406682330370277\ 99221994448419271992141, 448362074377990766195021286896356851920786439358786160\ 1237509388045184203863238979103610871462453373402268947216757552120565243118166\ 3516144458210435782115810038749325395856156029804972765105063247818630696960, 3583714208831485840578896294599155651628619947433246571053255848726194672796341\ 3797051012092580679764655712393855284036697524319773457805027103324323042817776\ 31855744946231369871312625188116653155965568635636321, 286476315985921240993141\ 4870223879677228862456446998356026407059298275383926229624749994926096428088666\ 7062656625857394777545994850016634523095135187169761380329534674073331629828580\ 8885978063233356671176929732056, 2290308822790051854754677990629701482655958335\ 2616908139295263144823788880951028738965830482073085086439042037053313923924317\ 2789644978769687267557733981383524338071233603547722299112132835801809387507985\ 80068155797, 183125352334330098712122561662980426921208778350383462145410463826\ 4680179881143231661424395130380893017405128221705229819522338336181968817515748\ 090199076848180385444685918581083991272970802462712881475178279087266880, 14643\ 7114518221480891687503635204462755058110770323827130354247029644748925830135536\ 3856305109230321708681548409807336076079089089009058628912165372672757248446099\ 39587950862100666386129499533089258783162245410036572009, 117111971969687057427\ 0210254881556032962353495981146853207702229157516642153370685912737642204566606\ 2622071936538669116702420910986977180382935134321807589362243345077789896456362\ 637854345505219281498852113714539128858584, 93669470859584928227586717292422075\ 4015683270581229363183561062069985279944377170279872066780284103646223763907752\ 9534538659774615653224270142107012981683869824377149984515004313945260364804197\ 62859144179496033914773546461, 749274068105905168876166859985374134357379312133\ 1831205101121921103303505234356840444257709229051392165161228233718375196837262\ 6331008832925675985575471372334720937227594844354248627281156257248212956975928\ 225249370721513984, 59941614933679463626504967326006653987145863699362348228589\ 6743612342900420147848389655133228300405590320581518634329459320050047004324130\ 1102778842845465216882867492527080317380210940911812087410683452703985854156802\ 265119409, 47957941878972464875437332732908289794308358487242825334257123056236\ 4961222307155241434265194530226984658295860898472846937787947315279893702655155\ 6400918466165360488699334143543149321203782032890089033510719518011351811193983\ 60, 383739293435178647672194769351789761649874413148601781731438951369594823344\ 8372925665904260246149082142480192679543186361849077910944584442810671478465015\ 4930274080324523686547229226676317207514738179545473053440563981693915851397, 3070824078000560388022474386115571855815928223911271993834520574732080220966078\ 7046194765660857372819020127474436174740212536376118680458802696278073129149363\ 29550638122220941587080990682733076498066391855408106775628122493604885824, 245\ 7626295136293652765756383256173311519115250867946831355236480001623985468918024\ 7541761193098509227952733181243998219759311118764070837071077317520055598504976\ 0363175987350680787500961016980930193166001019553647945143818583004710041, 1967\ 0632412423664537517753243083900445150435222365273043230045469432966150066154363\ 7065432946903757213919813514866071320953991101387069752544874939092570501758286\ 11133644152352531051232026965651420729669238667800882644205567654807233048, 157\ 4569061362861436030957981531576384763064351327180209417047625839369250751169633\ 9781484516929509891676918490128446446910000111440785998433443895571639861378908\ 13011820117846846999782962934533532469073189759728293070979551374424127301661, 1260507334516289122259511103636584953267522595758928315947565614780956578401661\ 6505972324038623764096423084659721998549557305618833211923474534729772635090193\ 8913017468751723915570829631330348441794960905063977418799793973100218202380566\ 528, 10091801196146249416397185495065889785550803661735153048983813791346720333\ 8380607669465873482085458855696141620594927297202455419978448727416955356621451\ 4341931018279416583887832077261173021621426832566164280973840827279782709915881\ 9768634625, 8080365962873125590289273560015133494621139501528831264195251674800\ 6888030607018726196294803444801424711814516634898530441502287861589234133381499\ 6730407840168867516239914473544772709590295699252067342620591514214888711947056\ 411007154106831448, 64704101621165712779466560043844027692907736782934228665968\ 7998933635985154050447370955969795533625097164457548305899477153353204089003209\ 3412798336632007585629459073800745306112029288863444260384461706308392819496918\ 1989962985248175802394362197, 5181678289483604607342187357659336100505698542310\ 9800169385855952079597119494109550848237661697125071413229910351846096470159945\ 7963785291252644667698397591227305293023591760499041844767822044043674634784240\ 8235077273434718194193183444541313910336, 4149983426538839925133935842999297763\ 8321569805787844674008655401509410967557984593713085153516812441295507274678731\ 3211194750275985284178215063522807185289251960128205437333660064686803923170520\ 765926183381303059270496138382638986953643802160067849, 33239846015519822683306\ 8874882951732051053161566176721634834861840558723304827047817831190746729377835\ 6573697283804832766085983266341857238633994660389821971658637682706085613816956\ 7901709342224678679131697482990930124402494491175659764349171116001496, 2662611\ 6206961049933674054716288921263260412313763748185026192913979934486383794171939\ 6907653397366631306602753694635742255777641921451745626927350232025498947145950\ 7349731896247377018272225061932350358881326025049628055049294153620360071973035\ 452359421, 21330070383484011501444379604071250784758660040858944722511781659462\ 0578628346743411220016706711041207616564731106269275314118505045096308106516098\ 0601768807052950841353100452877348965666483894408952911542197950547516070896596\ 31710194229211637340011145728, 170888122260835095803664806950007981332418643767\ 0130899871808284558483456341155972790761191787313800753165517615682325305499655\ 8560488602713501945634537752655290914067406207106327571417424687745498036674456\ 436026110175479092971218817558151261014852444227985, 13691974461984292604199859\ 2500440237429016722969734674682495200305948512125237053588655918399266558912759\ 8336003885301047641364280502532848875893985387816005904252155780050180296708526\ 016302056059801200934369849295608699506394112867790002072071901145603212440, 10\ 9712066997091643087962492545894433468584647525595467147826243328462876878091882\ 5409126630769399879068702181681647729352218057874287486921539330989309124073742\ 7443109107098338767259136560209011851133702916365999653851610392257962829105947\ 6501012566648113551733, 8791770441383329722422199075589339407048153964657858407\ 3775793131444315340887615388658638101539598535665876552050544437365623809450881\ 7481046492625387329122489476265958925319991125055825626022122266552693941575094\ 2211546581403277872249755763269032479650216035648, 7045818148128715655179262090\ 1338146532299012569656685251683136139034245273043730899393841325305302890537626\ 0476364594782294008116351395074540185565936439590960504061975751806020623499179\ 362072794647713891290957677976864379507873513872925671761608835968975067155417, 5647018355506196211433150770749725858723988418945725075385722208891431825215048\ 9752855200010786526875595270408856086659398979705714350497798116225995647807008\ 0838986531134345875770558448108272203316727690403693767351327613317059299970544\ 85321548080578040126389617752, 452625679941866052315328612569957965184259976613\ 1877776905462517883943306453127866097289465750974810223077276845361631712830147\ 8681838860550309485737742949554244089961965853414689344606781346106403917691208\ 47668998927540051650523439446495754992793137954761845092737357, 362819778431962\ 0877695455275676808393485252814576790645350032975740618415506466947640710082246\ 9432673200842956422108576525938915517082032333501784600396942523683415565652038\ 8635790560831843097652659750221236352067763772160994002321949668588618666075833\ 200362408246759424, 29085333690455109056766890138687902816039160548377242020743\ 4388108747579284908367718561383750634783163628965271409549448420642955971093270\ 9552682935487208040464690460915410135527014321908620429267355555141907299749320\ 6864596474705930825520367253889353230242128077516636705, 2331782565963045777760\ 3800369868519761525181093510663144577018641685951004226867998683776135919890032\ 9773543582121822432143662209424211630823604475459521461424760096282119831077025\ 5917307559573105175340155404403425256546729207824197918042671122884712090440825\ 177006091934616, 18695304054438320801911884354338651452197392631824802633379680\ 5417290118378543551933439431406495648859209795166356250843178090409614705590611\ 2564858066201480115169784820626494457415687659766136161570406185020689717821199\ 13890382859134348205509249965939756878112474861232710757, 149901881352968983939\ 5522293397990426683721885460698662768129402438301423576720565319872014536955837\ 3066659006489047588654498373565283644775207970467963168570785059349283102910438\ 2762707794775849496738009350800508088399259623672920147234099489517128259954712\ 60838186647916890688, 120201890614396043292152088357771575561484315957331258199\ 8232001840450729584949447275001419398692935370254667854056646885050369077616181\ 1985913731122408814782312585636352682983951465375670988396362007945190747369011\ 31999482651052487460650274053445226538455297816341370792168019785, 963928412350\ 6855464449809920875436150950158978724735593078362062580367081608324773260416447\ 5537240915090867972917971623929500252975810303606489303752372534473043627885386\ 4436039209313509777765737707229688131987871984640935269854081481874480275117805\ 84740034918665999836819456968088, 773049215619734711737477996636737093551836158\ 4887957327539893447470040317300038938247284501679175881750994473384343419957541\ 0241424107825147360729741931073643856081804626725297103435345215948991331623351\ 9359081081520431387822205346992144971806148707416350933510784091736410452709031\ 7, 6200090155802414574963002186363732062698143727106792288578063930961372304129\ 2210339914169979839281968071600352254683404531216408454522578254397460637422561\ 6103916174017158718088276984133755409680874526577812779758964226036196777184045\ 558546786851600776591224549613246799646881870762496, 49729831827414204843565862\ 0056484566753505119393736165263259760101209882504232142360486868719906919764798\ 2901463034901710678025294227349396119423337305744185047609606360154982540346481\ 4997710236263669257861037135299341744902323155194040188376598749728900583148262\ 455830225353552724863025, 39889972375891613617936323283599955783877853073347496\ 3316202836114740806316737052061221172237832297217362159812075964566983373912705\ 9734323381673162519969545306726055994300165533323128690406330285832859378045967\ 203810200116252772157122341467047827164471568791050174169914928119764117138456, 3199910862166606912277039254677820822698798048582346913021255972216009631801265\ 7076644870788999105453345753166533590120923924011680346876235422526513632842936\ 3680518530708988284274394878347930726829330777235558713466009270450380601887413\ 641581818448251982158977339359134500219882489714959477, 25670780528212513771497\ 5786239322500941241903338947911715035803494535811745645970033249431923813601477\ 9397804399480031510560882509634808745840071780396315018739988775523286994877341\ 0144298625174323900333314829399166065848597608588698203562236466485304667688100\ 948299959508762765145378143752000, 20595245326726584758610229212329448286354824\ 5862087620972762811149921147073211253255414204396928994324544086125148245694571\ 4967366135155771605346006752890211231889808217867801046918983421527497699318215\ 6639628472974269408006773362412833486225942350295065998951745853259364139791080\ 47329992029945, 165242299274620869147709542011007038985041685512469910763517631\ 2666259162833465228639426153414841113494546131429806952390996776391367704742357\ 9959493293373707538351484171796994680284131697107609412683574659372894239132405\ 6340683245461602523496132126071400038458249536013876688274467954743726131288, 1325871760716780556386436615872058537411437268448479124741252070527340936586750\ 6067893265399350698138529271459899897810509808714415544124534247214517903067705\ 3919990127850311769321814382469491028176443145312166226059110810057814106466896\ 44956319603742937831916338068889473123005840642313569999898157, 106391643517342\ 7606703119661956209811986342604111601457596843271127467762084662827790711336845\ 0033727607117475889803378026638962040474731621640569943321407714941476540708075\ 6734955812396284518795976259267514901874237355504728421602907948342743988391929\ 0358715337467434742316043994769648852923099283456]: if d=3 and N<=nops(L1) then RETURN([op(1..N,L1)]): elif d=4 and N<=nops(L2) then RETURN([op(1..N,L2)]): else SVTseq(d,N): fi: end: #EmpirAsyC(L1,n,mu,t,M): inputs a list L1 of positive numbers believed to gave asympotics of the form #mu^n*n^t, finds the constant using asympotics of ofer M. Try: #EmpirAsyC([seq(5*3^i/i*(1+4/i),i=1..40)],3,-1,5); EmpirAsyC:=proc(L1,n,mu,t,M) local eq,var,c,lu,i: if nops(L1)<20 or M>nops(L1)/2 then print(`List too short`): RETURN(FAIL): fi: lu:=mu^n*n^t*add(c[i]/n^i,i=0..M): var:={seq(c[i],i=0..M)}: eq:={seq(L1[i]-subs(n=i,lu),i=nops(L1)-M..nops(L1))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: evalf(subs(var,c[0]),10): end: #SVTpaper(d,n,Kama,M,K): inputs an integer d, a symbol n, a positive integer Kama, a positive integer M, and a large positive integer K`): #print(`outputs an article, that tells you first the first Kama terms of the sequence enumerating #the number of singular vector tuples of an n by n by ... by n (d times) tensor, using the constant-term formula in the article: #Shmuel Friedland and Giorgio Ottaviani, The number of singular vector tuples and uniqueness of best rank-one approximation of tensors # Found. Comput. Math. 14 (2014), no. 6, 1209-1242. #Then it tries to find a homogeneous linear recurrence equation satisfies by that sequence using that many terms. If it fails #the paper is finished there. If it succeeds, it states it, and tries to use it to find an asympotic formula. #if it succeeds it tells you so. SVTpaper:=proc(d,n,Kama,M,K) local gu,a,ope,t0,Ope,ru,mu,lu,N,i,ku,C: t0:=time(): print(``): gu:=SVTseqPC(d,Kama): print(`The Number of Singular Vector Tuples of a general`, [n$d], `tensor . `): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`Theorem: Let a(n) be the number of singular vector tuples of a`, d , `-dimensional tensor whose dimensions are all`, n): print(``): print(`The first`, Kama, `terms of that sequence are`): print(``): print(gu): ope:=Findrec(gu,n,N): if ope=FAIL then print(`While we do know that the sequence is P-recursive, i.e. satisfies SOME homogeneous linear recurrence equation with constant coefficients`): print(`the number of terms in the input`, Kama, `does not suffice. If you are really interested, try to make it larger. `): print(``): print(`nevertheless let's do some non-rigorous estimates of the asymptotics.`): ku:=Zinn(gu): print(`It seems that the sequence behaves like`, C*ku[2]^n*n^ku[1], `for some constant C`): if evalf(ku[2]/(d-1)^d-1)<0.1 and evalf(ku[1]+(d-1)/2)<0.2 then print(`These are close to our conjectured asymptotics`): print(C*((d-1)^d)^n*n^(-(d-1)/2) ): print(``): print(`for some constant C. `): # print(`estimating the constant we get that it is probably approximately`): # print(``): # C0:=EmpirAsyC(gu,n,(d-1)^d,-(d-1)/2,M): # print(C0): # print(``): # print(`Hence a non-rigorous estimate of the asymptotics is`): # print(``): # print(C0*((d-1)^d)^n*n^(-(d-1)/2) ): #print(``): fi: print(`This ends this paper that took`, time()-t0, `seconds to generate. `): print(`Have a good day. `): RETURN(): fi: print(`Theorem 1: The sequence a(n) satisfies the following linear recurrence equation with polynomial coefficients`): print(``): print(add(coeff(ope,N,i)*a(n+i),i=0..degree(ope,N))=0): print(``): print(`subject to the initial conditions`): print(``): lprint(seq(a(i)=gu[i],i=1..degree(ope,N))): print(``): print(`and in Maple notation`): print(``): lprint(add(coeff(ope,N,i)*a(n+i),i=0..degree(ope,N))=0): print(``): print(`Just for fun, the`, K, `-th term of this sequence, i.e. the number of singular vector tuples of`): print(`a `, d, `-dimensional tensor, all whose dimensions are equal to`, K, `is equal to `): print(``): lprint(SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],K)[K]): print(``): Ope:=expand(numer(ope)): ru:=sort([solve(lcoeff(Ope,n),N)]): mu:=ru[nops(ru)]: lu:=AsyC(Ope,n,N,mu,M,[op(1..degree(Ope,N),gu)],K): if lu=FAIL and nops(ru)>1 then mu:=ru[nops(ru)-1]: lu:=AsyC(Ope,n,N,mu,M,[op(1..degree(Ope,N),gu)],K): fi: if lu=FAIL then print(`We were unable to find the asymptotic expansion, but we know it exists. You are welcome to try it yourself.`): print(``): print(`This ends this paper that took`, time()-t0, `seconds to generate. `): print(`Have a good day. `): RETURN(): fi: print(`Theorem 2: The asymptotic formula for a(n), to order`, M, ` is given by `): print(lu): print(``): print(`and in Maple notation`): print(``): lprint(lu): print(``): print(``): print(`This ends this paper that took`, time()-t0, `seconds to generate. `): print(`Have a good day. `): end: #MMMT(A,x): the generating function given by MacMahon's master theorem in the list of variables x. #Try: #MMMT([[0,1,1],[1,0,1],[1,1,0]],[x1,x2,x3]): MMMT:=proc(A,x) local M,M1,i,j: M:=[]: for i from 1 to nops(A) do M1:=[]: for j from 1 to nops(A[i]) do if i=j then M1:=[op(M1),1-A[i][j]*x[i]]: else M1:=[op(M1),-A[i][j]*x[i]]: fi: od: M:=[op(M),M1]: od: 1/det(M): end: SVTseq3:=proc(N) local i: [seq(A3f(i,i,i),i=0..N-1)]: end: SVTseq4:=proc(N) local i: [seq(A4f(i,i,i,i),i=0..N-1)]: end: S1:=proc(): op(sort([args])):end: #B3f(a1,a2,a3): the coefficient of x1^a1*x2^a2*x3^a3 in 1/(1-2*(x1*x2+x1*x3+x2*x3)-2*x1*x2*x3) B3f:=proc(a1,a2,a3) option remember: if a1<0 or a2<0 or a3<0 then 0: elif a1=0 and a2=0 and a3=0 then 1: else B3f(S1(a1-1,a2-1,a3) )+B3f(S1(a1-1,a2,a3-1))+ B3f(S1(a1,a2-1,a3-1))+2*B3f(S1(a1-1,a2-1,a3-1)): fi: end: #A3f(a1,a2,a3): the Sum of B3(b1,b2,b3) for 0<=b1<=a1,0<=b2<=a2,0<=b3<=a3, A3f:=proc(a1,a2,a3) option remember: if a1<0 or a2<0 or a3<0 then 0: elif a1=0 and a2=0 and a3=0 then 1: else A3f(S1(a1-1,a2,a3))+A3f(S1(a1,a2-1,a3))+A3f(S1(a1,a2,a3-1)) -A3f(S1(a1-1,a2-1,a3))-A3f(S1(a1-1,a2,a3-1))- A3f(S1(a1,a2-1,a3-1))+ A3f(S1(a1-1,a2-1,a3-1))+B3f(S1(a1,a2,a3)): fi: end: #B4f(a1,a2,a3,a4): the coefficient of x1^a1*x2^a2*x3^a3*x4^a4 in 1/(1-2*e_2(x1,x2,x3,x4)-3*e3(x1,x2,x3,x4)-4e4(x1,x2,x3,x4)) B4f:=proc(a1,a2,a3,a4) option remember: if a1<0 or a2<0 or a3<0 or a4<0 then 0: elif a1=0 and a2=0 and a3=0 and a4=0 then 1: else B4f(S1(a1-1,a2-1,a3,a4))+B4f(S1(a1-1,a2,a3-1,a4))+ B4f(S1(a1,a2-1,a3-1,a4)) +B4f(S1(a1-1,a2,a3,a4-1))+B4f(S1(a1,a2-1,a3,a4-1))+B4f(S1(a1,a2,a3-1,a4-1)) +2*(B4f(S1(a1-1,a2-1,a3-1,a4))+B4f(S1(a1-1,a2-1,a3,a4-1))+B4f(S1(a1-1,a2,a3-1,a4-1))+B4f(S1(a1,a2-1,a3-1,a4-1))) +3*B4f(S1(a1-1,a2-1,a3-1,a4-1)): fi: end: #A4f(a1,a2,a3,a4): the Sum of B4(b1,b2,b3,b4) for 0<=b1<=a1,0<=b2<=a2,0<=b3<=a3, 0<=b4<=a4 A4f:=proc(a1,a2,a3,a4) option remember: if a1<0 or a2<0 or a3<0 or a4<0 then 0: elif a1=0 and a2=0 and a3=0 and a4=0 then 1: else A4f(S1(a1-1,a2,a3,a4))+A4f(S1(a1,a2-1,a3,a4))+A4f(S1(a1,a2,a3-1,a4))+ A4f(S1(a1,a2,a3,a4-1)) -A4f(S1(a1-1,a2-1,a3,a4))-A4f(S1(a1-1,a2,a3-1,a4))- A4f(S1(a1,a2-1,a3-1,a4))- A4f(S1(a1-1,a2,a3,a4-1))-A4f(S1(a1,a2-1,a3,a4-1))-A4f(S1(a1,a2,a3-1,a4-1)) +A4f(S1(a1-1,a2-1,a3-1,a4))+A4f(S1(a1-1,a2-1,a3,a4-1))+A4f(S1(a1-1,a2,a3-1,a4-1))+A4f(S1(a1,a2-1,a3-1,a4-1)) -A4f(S1(a1-1,a2-1,a3-1,a4-1))+B4f(S1(a1,a2,a3,a4)): fi: end: #Cnk(n,k): all the 0-1 vectors of length n with k 1's Cnk:=proc(n,k) local gu,gu1,gu11: option remember: if n=0 then if k=0 then RETURN({[]}): else RETURN({}): fi: fi: if k>n or k<0 then RETURN({}): fi: gu1:=Cnk(n-1,k): gu:={seq([op(gu11),0],gu11 in gu1)}: gu1:=Cnk(n-1,k-1): gu:=gu union {seq([op(gu11),1],gu11 in gu1)}: end: #Bdf(a): f(a), where a is a list of non-negative integers, the coefficient of x1^k1*...xd^kd #in the rational function 1/(1-e_2(x1, ..., xd)- 2*e_3(x1, ..., xd)- (d-1)e_d(x1, ..., xd)) #Try: #Bdf([3,3,3]); Bdf:=proc(a) local d,gu,gu1,k,su: option remember: d:=nops(a): if min(op(a))<0 then RETURN(0): fi: if a=[0$d] then RETURN(1): fi: su:=0: for k from 2 to d do gu:=Cnk(d,k): su:=su+(k-1)*add(Bdf(sort(a-gu1)),gu1 in gu): od: su: end: #Adf(a): c(a), where a is a list of non-negative integers, the coefficient of x1^k1*...xd^kd #in the rational function 1/(1-e_2(x1, ..., xd)- 2*e_3(x1, ..., xd)- (d-1)e_d(x1, ..., xd))*(1/((1-x1)*...*(1-xd)) #Try: #Adf([3,3,3]); Adf:=proc(a) local d,gu,gu1,k,su: option remember: d:=nops(a): if min(op(a))<0 then RETURN(0): fi: if a=[0$d] then RETURN(1): fi: su:=Bdf(a): for k from 1 to d do gu:=Cnk(d,k): su:=su+(-1)^(k-1)*add(Adf(sort(a-gu1)),gu1 in gu): od: su: end: #Bdr(a,r): f(a), where a is a list of non-negative integers of length d, say, outputs the the coefficient of x1^a[1]*...xd^a[d] #in the rational function 1/(1-e_2(x1, ..., xd)- 2*e_3(x1, ..., xd)- (d-1)e_d(x1, ..., xd))^r #Try: #Bdr([3,3,3]); Bdr:=proc(a,r) local d,gu,gu1,k,su: option remember: if r=1 then RETURN(Bdf(a)): fi: d:=nops(a): if min(op(a))<0 then RETURN(0): fi: if a=[0$d] then RETURN(1): fi: su:=Bdr(a,r-1): for k from 2 to d do gu:=Cnk(d,k): su:=su+(k-1)*add(Bdr(sort(a-gu1),r),gu1 in gu): od: su: end: