/* Copyright (C) 2007 Christian G. Bower transforms_pari.txt Version 0.9 This file is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY WHATSOEVER. Check the License for details. You can find it at: http://www.gnu.org/licenses/gpl-3.0.txt */ /* This is a file of sequence transforms (or functions). Sequences are encoded as vectors. If vector A represents sequence a, then A[1] = a_0, A[n+1] = a_n Some transforms operate on matrices. */ /* Core transforms and operations * First are the most fundamental and most used functions */ /* binomial transform */ { trv_bin(A) = forstep(i=#A-1,0,-1,A[i+1]=sum(j=0,i,binomial(i,j)*A[j+1])); A } /* inverse binomial transform */ { trv_i_bin(A) = for(i=1,#A-1,for(j=0,i-1,A[i+1]-=binomial(i,j)*A[j+1])); A } /* characteristic function of A * output sequence b_n = 1 if n if exist k : a_k == n * otherwise b_n = 0 * 2nd parameter is length of output sequence */ { trv_char(A, n) = local(B=vector(n),v); for(i=1, #A, v=A[i]; if(type(v)=="t_INT" && v>=0 && v1, A[2]=A[2]*A[1]+1); for(i=3,#A, A[i]=A[i]*A[i-1]+A[i-2]); A } /* continued fraction denominator transform * output sequence of denominators of fraction estimates for continued * fraction coefficients in input sequence */ { trv_cofr_den(A) = A[1]=1; for(i=3,#A, A[i]=A[i]*A[i-1]+A[i-2]); A } /* cycle transform * apply cycle structure to sequence of unlabeled objects * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 * function uses trv_i_moebius and trv_lyndon * function indirectly uses trv_i_euler */ { trv_cycle(A) = A = trv_i_moebius(trv_lyndon(A)); A[1] = 1; A } /* inverse cycle transform * this function ignores a_0 (A[1]) and assumes it to be 1 * the output sequence always has b_0=0 * function uses trv_moebius and trv_i_lyndon * function indirectly uses trv_euler */ { trv_i_cycle(A) = trv_i_lyndon(trv_moebius(A)); } /* difference transform * b_n = a_n - a_{n-1} * inverse of trv_psum */ { trv_diff(A) = vector(#A,i,if(i==1,A[1],A[i]-A[i-1])) } /* euler transform * apply set structure to sequence of unlabeled objects * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 */ { trv_euler(A) = local(B=vector(#A-1,n,1/n),C); B = dirmul(vecextract(A,"2.."),B); C = exp(Ser(concat(0,B))); for(i=1, #A, A[i] = polcoeff(C,i-1)); A } /* inverse euler transform * this function ignores a_0 (A[1]) and assumes it to be 1 * the output sequence always has b_0=0 */ { trv_i_euler(A) = local(B=vector(#A-1,n,1/n),C); A[1] = 1; C = log(Ser(A)); A=vecextract(A,"2.."); for(i=1, #A, A[i] = polcoeff(C,i)); A = dirdiv(A,B); concat(0,A) } /* euler_2 transform * apply (set of 2 elements) structure to sequence of unlabeled objects * aka symmetric square * this function is not invertible */ { trv_euler_2(A) = local(x='x, B=Ser(A,x)); B = (B^2+subst(B,x,x^2))/2; for(i=1, #A, A[i] = polcoeff(B,i-1)); A } /* euler_2 transform for matrices */ { trm_euler_2(A) = local(x='x,y='y,m=matsize(A)[1],n=matsize(A)[2],B=O(x^m)+O(y^n), B2=O(x^m)+O(y^n)); for(i=0, m-1, for(j=0, n-1, B+=A[i+1,j+1]*x^i*y^j); ); for(i=0, (m-1)\2, for(j=0, (n-1)\2, B2+=A[i+1,j+1]*x^(2*i)*y^(2*j)); ); B = (B^2+B2)/2; for(i=0, m-1, for(j=0, n-1, A[i+1,j+1] = polcoeff(polcoeff(B,j,y),i,x)); ); A } /* euler product transform * apply set structure to sequence of unlabeled objects * to assemble them by the product of their sizes * this function ignores a_0, a_1 (A[0], A[1]) and assumes them to be 0 * the output sequence always has b_0=0, b_1=1 */ { trv_euler_p(A) = local(B=vector(#A-1),j,s); A=vecextract(A,"2.."); for(i=2, #A, s = bigomega(i)*A[i]; j=i; while(j<=#A, B[j] += s; j*=i; ); ); A[1]=1; for(i=2, #A, A[i] = sumdiv(i, j, A[j]*B[i/j])/bigomega(i); ); concat(0,A) } /* inverse euler product transform * this function ignores a_0, a_1 (A[0], A[1]) and assumes them to be 0, 1 * the output sequence always has b_0=b_1=0 */ { trv_i_euler_p(A) = local(B=vector(#A-1),j); A=vecextract(A,"2.."); for(i=2, #A, B[i] = bigomega(i)*A[i]-sumdiv(i, j, A[j]*B[i/j]); ); for(i=2, sqrtint(#A), j=i^2; while(j<=#A, B[j] -= B[i]; j*=i; ); ); for(i=2, #A, A[i]=B[i]/bigomega(i)); A[1]=0; concat(0,A) } /* exponential transform * apply set structure to sequence of labeled objects * inverse is trv_log (logarithmic transform) * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 */ { trv_exp(A) = local(x='x,B); B = sum(i=1, #A-1, A[i+1]*x^i/i!); B=exp(B+O(x^#A)); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* invert transform * apply list structure to sequence of unlabeled objects * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 * trv_invert_x is a version not making that assumption */ { trv_invert(A) = local(x='x, B=Ser(A,x)-A[1]); B = 1/(1-B); for(i=1, #A, A[i] = polcoeff(B,i-1)); A } /* inverse invert transform * this function ignores a_0 (A[1]) and assumes it to be 1 * the output sequence always has b_0=0 */ { trv_i_invert(A) = local(x='x, B=Ser(A,x)-A[1]+1); B = 1-1/B; for(i=1, #A, A[i] = polcoeff(B,i-1)); A } /* invert product transform * apply list structure to sequence of unlabeled objects * to assemble them by the product of their sizes * this function ignores a_0, a_1 (A[0], A[1]) and assumes them to be 0 * the output sequence always has b_0=0, b_1=1 */ { trv_invert_p(A) = A=-vecextract(A,"2.."); A[1]=1; A=dirdiv(vector(#A,n,n==1),A); concat(0,A) } /* inverse invert product transform * this function ignores a_0, a_1 (A[0], A[1]) and assumes them to be 0, 1 * the output sequence always has b_0=b_1=0 */ { trv_i_invert_p(A) = A=vecextract(A,"2.."); A[1]=1; A=-dirdiv(vector(#A,n,n==1),A); A[1]=0; concat(0,A) } /* invert labeled transform * apply list structure to sequence of labeled objects * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 */ { trv_invert_l(A) = local(x='x, B); B = sum(i=1, #A-1, A[i+1]*x^i/i!); B = 1/(1-B+O(x^#A)); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* inverse invert labeled transform * this function ignores a_0 (A[1]) and assumes it to be 1 * the output sequence always has b_0=0 */ { trv_i_invert_l(A) = local(x='x, B); B = sum(i=1, #A-1, A[i+1]*x^i/i!); B = 1-1/(1+B+O(x^#A)); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* least inverse transform * b_n is the least k : a_k == n or fill value if no such k exists * 2nd parameter is length of output sequence * 3rd parameter is fill value */ { trv_leastinv(A, n, f=-1) = local(B=vector(n,x,-1),v); for(i=1, #A, v=A[i]; if(type(v)=="t_INT" && v>=0 && vb_{k-1} * this function is not invertible * output length is determined by the function */ { trv_record(A) = local(B=[A[1]]); for(i=2, #A, if(A[i]>B[#B], B=concat(B,A[i])); ); B } /* indices of record values of A * output sequence b_n is indices of a where the record values * of trv_record occurred * this function is not invertible * output length is determined by the function */ { trv_record_ix(A) = local(B=[0],j=1,k=A[1]); for(i=2, #A, if(A[i]>k, B=concat(B,i-1); j=i; k=A[i]; ); ); B } /* revert transform * apply reversion of series as ogf * this function is its own inverse * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=0 */ { trv_revert(A) = local(x='x, B=Ser(A,x)-A[1]); B = serreverse(B); for(i=1, #A, A[i] = polcoeff(B,i-1)); A } /* revert_egf transform * apply reversion of series as egf * this function is its own inverse * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=0 */ { trv_revert_egf(A) = local(x='x,B); B = sum(i=1, #A-1, A[i+1]*x^i/i!); B = serreverse(B+O(x^#A)); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* stirling transform * apply given labeled structure to non-empty sets */ { trv_stirling(A) = local(x='x,B); B = sum(i=0, #A-1, A[i+1]*x^i/i!); B = subst(B+O(x^#A),x,exp(x+O(x^#A))-1); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* inverse stirling transform */ { trv_i_stirling(A) = local(x='x,B); B = sum(i=0, #A-1, A[i+1]*x^i/i!); B = subst(B+O(x^#A),x,log(1+x+O(x^#A))); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* weigh transform * apply set structure to sequence of unlabeled objects * injectively (or asymmetrically) * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 * trv_weigh_x is a version not making that assumption */ { trv_weigh(A) = local(i=-1,B=vector(#A-1,n,i=-i;i/n),C); B = dirmul(vecextract(A,"2.."),B); C = exp(Ser(concat(0,B))); for(i=1, #A, A[i] = polcoeff(C,i-1)); A } /* inverse weigh transform * this function ignores a_0 (A[1]) and assumes it to be 1 * the output sequence always has b_0=0 * trv_i_weigh_x is a version not making that assumption */ { trv_i_weigh(A) = local(i=-1,B=vector(#A-1,n,i=-i;i*1/n),C); A[1] = 1; C = log(Ser(A)); A=vecextract(A,"2.."); for(i=1, #A, A[i] = polcoeff(C,i)); A = dirdiv(A,B); concat(0,A) } /* weigh labeled transform * apply set structure to sequence of labeled objects * injectively (or asymmetrically for L-species) * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 * function uses trv_exp */ { trv_weigh_l(A) = local(B=vector(#A-1),i,s); A=vecextract(A,"2.."); for(k=1, #A, s=1; i=1; forstep(j=k,#A,k, if(i>1, s*=(1-i)*binomial(j-1,j-k)); B[j] += s*A[k]; i++; ); ); trv_exp(concat(0,B)) } /* inverse weigh labeled transform * this function ignores a_0 (A[1]) and assumes it to be 1 * the output sequence always has b_0=0 * function uses trv_log */ { trv_i_weigh_l(A) = local(i,s); A=trv_log(A); A=vecextract(A,"2.."); for(k=1, #A\2, s=1; i=2; forstep(j=k+k,#A,k, s*=(1-i)*binomial(j-1,j-k); A[j] -= s*A[k]; i++; ); ); concat(0,A) } /* egf to sequence * output vector of sequence whose exponential generating function * is given by the input series ("t_SER") */ { trsv_egf(A) = vector(#A, i, polcoeff(A,i-1)*(i-1)!) } /* ogf multiplication operation * multiplies ogf of sequences * sequence length is the shorter of the 2 */ { opv_mul_ogf(A, B) = local(C=Ser(A,x)*Ser(B,x)); A=vector(min(#A,#B)); for(i=1, #A, A[i] = polcoeff(C,i-1)); A } /* ogf division operation * divides ogf of sequences * sequence length is the shorter of the 2 */ { opv_div_ogf(A, B) = local(C=Ser(A,x)/Ser(B,x)); A=vector(min(#A,#B)); for(i=1, #A, A[i] = polcoeff(C,i-1)); A } /* egf multiplication operation * multiplies egf of sequences * sequence length is the shorter of the 2 */ { opv_mul_egf(A, B) = local(C=vector(min(#A,#B))); forstep(i=#C-1,0,-1,C[i+1]=sum(j=0,i,binomial(i,j)*A[j+1])*B[i-j+1]); C } /* egf:dgf multiplication operation * returns sequence whose egf is b_1*A(x) + b_2*A(x^2/2!) + ... * sequence length is the shorter of the 2 * note that this operation is not commutative */ { opv_mul_ed(A, B) = local(C,i,s); C=vector(min(#A,#B)); for(k=1, #C-1, s=1; i=1; forstep(j=k,#C-1,k, s*=binomial(j-1,j-k); C[j+1] += s*A[i+1]*B[k+1]; i++; ); ); C } /* Other transforms * Following are additional more specialized functions */ /* boustrophedon transform */ { trv_bous(A) = local(x='x,B); B = sum(i=0, #A-1, A[i+1]*x^i/i!); B=(tan(x+O(x^#A))+1/cos(x+O(x^#A)))*(B+O(x^#A)); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* inverse boustrophedon transform */ { trv_i_bous(A) = local(x='x,B); B = sum(i=0, #A-1, A[i+1]*x^i/i!); B=(B+O(x^#A))/(tan(x+O(x^#A))+1/cos(x+O(x^#A))); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* carlitz transform * apply list structure to sequence of unlabeled objects * insisting that neighbors be different object types * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 * trv_carlitz_x is a version not making that assumption * function uses trv_invert */ { trv_carlitz(A) = local(i=-1,B=vector(#A-1,n,i=-i;i)); A = concat(0,dirmul(vecextract(A,"2.."),B)); trv_invert(A) } /* inverse carlitz transform * this function ignores a_0 (A[1]) and assumes it to be 1 * the output sequence always has b_0=0 * function uses trv_i_invert */ { trv_i_carlitz(A) = local(i=-1,B=vector(#A-1,n,i=-i;i),C); A = trv_i_invert(A); concat(0,dirdiv(vecextract(A,"2.."),B)); } /* extended version of carlitz transform * it does not assume a_0 = A[1] == 0 * function uses trv_invert_x */ { trv_carlitz_x(A) = local(i=-1,B=vector(#A-1,n,i=-i;i)); A = concat(A[1]/2,dirmul(vecextract(A,"2.."),B)); trv_invert_x(A) } /* chain transform * apply chain (undirected list) structure to sequence of * unlabeled objects * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 */ { trv_chain(A) = local(x='x, B=Ser(A,x)-A[1], B2 = subst(B,x,x^2)); B = (2-B^2)/(2*(1-B)) + B2*(1+B)/(2*(1-B2)); for(i=1, #A, A[i] = polcoeff(B,i-1)); A } /* chain transform for matrices * this function ignores a_0,0 (A[1,1]) and assumes it to be 0 * the output sequence always has b_0,0=1 */ { trm_chain(A) = local(x='x,y='y,m=matsize(A)[1],n=matsize(A)[2],B=O(x^m)+O(y^n), B2=O(x^m)+O(y^n),j=1); for(i=0, m-1, for(j=j, n-1, B+=A[i+1,j+1]*x^i*y^j); j=0; ); for(i=0, (m-1)\2, for(j=j, (n-1)\2, B2+=A[i+1,j+1]*x^(2*i)*y^(2*j)); j=0; ); B = (2-B^2)/(2*(1-B)) + B2*(1+B)/(2*(1-B2)); for(i=0, m-1, for(j=0, n-1, A[i+1,j+1] = polcoeff(polcoeff(B,j,y),i,x)); ); A } /* chain injective transform * apply chain (undirected list) structure to sequence of * unlabeled objects injectively (or asymmetrically) * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 */ { trv_chain_j(A) = local(x='x, B=Ser(A,x)-A[1]); B2 = subst(B,x,x^2); B = (2-B^2)/(2*(1-B)) - B2*(1+B)/(2*(1-B2)); for(i=1, #A, A[i] = polcoeff(B,i-1)); A } /* chain injective transform for matrices * this function ignores a_0,0 (A[1,1]) and assumes it to be 0 * the output sequence always has b_0,0=1 */ { trm_chain_j(A) = local(x='x,y='y,m=matsize(A)[1],n=matsize(A)[2],B=O(x^m)+O(y^n), B2=O(x^m)+O(y^n),j=1); for(i=0, m-1, for(j=j, n-1, B+=A[i+1,j+1]*x^i*y^j); j=0; ); for(i=0, (m-1)\2, for(j=j, (n-1)\2, B2+=A[i+1,j+1]*x^(2*i)*y^(2*j)); j=0; ); B = (2-B^2)/(2*(1-B)) - B2*(1+B)/(2*(1-B2)); for(i=0, m-1, for(j=0, n-1, A[i+1,j+1] = polcoeff(polcoeff(B,j,y),i,x)); ); A } /* chain labeled transform * apply chain (undirected list) structure to sequence of labeled objects * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 */ { trv_chain_l(A) = local(x='x, B); B = sum(i=1, #A-1, A[i+1]*x^i/i!); B = (1/(1-B+O(x^#A))+1+B)/2; for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* inverse chain labeled transform * this function ignores a_0 (A[1]) and assumes it to be 1 * the output sequence always has b_0=0 */ { trv_i_chain_l(A) = local(x='x, B); B = sum(i=1, #A-1, A[i+1]*x^i/i!); B = 1+B-sqrt(1+B^2+O(x^#A)); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* chain "boy-girl" transform * apply chain (undirected list) structure to sequence of * unlabeled objects * insisting that neighbors be different object types * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 * function uses trv_carlitz, opv_mul_ogf */ { trv_chain_bg(A) = local(B,C,D=[0,1,0,-1]); B = trv_carlitz(A); C = dirmul(vecextract(B,"2.."), vector(#B-1,n,n==2)); A = vecextract(A,"2.."); A = dirmul(A, vector(#A,n,D[1+bitand(n,3)])); A = opv_mul_ogf(concat(0,A), concat(1,C)); A[1] = 1; (A+B)/2 } /* chebyshev transform * G.f. B(x) = (1-x^2)/(1+x^2)*A(x/(1+x^2)) * Used by Paul Barry for A100047 */ { trv_chebyshev(A) = local(x='x, B=Ser(A,x)); B = subst(B,x,x/(1+x^2+O(x^#A)))*(1-x^2+O(x^#A))/(1+x^2+O(x^#A)); for(i=1, #A, A[i] = polcoeff(B,i-1)); A } /* inverse chebyshev transform */ { trv_i_chebyshev(A) = local(x='x, B=Ser(A,x)); B = subst(B*(1+x^2+O(x^#A))/(1-x^2+O(x^#A)),x,serreverse(x/(1+x^2+O(x^#A)))); for(i=1, #A, A[i] = polcoeff(B,i-1)); A } /* continuant transform */ { trv_continuant(A) = local(B=vector(#A+1)); B[1]=1; B[2]=A[1]; for(i=2,#A, B[i+1]=A[i]*B[i]+A[i-1]); B } /* inverse continuant transform */ { trv_i_continuant(A) = local(B=vector(#A-1)); A[1]=1; B[1]=A[2]; for(i=2,#A-1, B[i]=(A[i+1]-B[i-1])/A[i]); B } /* cycle transform for matrices * this function ignores a_0,0 (A[1,1]) and assumes it to be 0 * the output sequence always has b_0,0=1 * function uses trm_i_moebius and trm_lyndon * function indirectly uses trm_i_euler */ { trm_cycle(A) = A = trm_i_moebius(trm_lyndon(A)); A[1,1] = 1; A } /* inverse cycle transform for matrices * this function ignores a_0,0 (A[1,1]) and assumes it to be 1 * the output sequence always has b_0,0=0 * function uses trm_moebius and trm_i_lyndon * function indirectly uses trm_euler */ { trm_i_cycle(A) = trm_i_lyndon(trm_moebius(A)); } /* cycle product transform * apply cycle structure to sequence of unlabeled objects * to assemble them by the product of their sizes * this function ignores a_0, a_1 (A[0], A[1]) and assumes them to be 0 * the output sequence always has b_0=0, b_1=1 * function uses trv_i_moebius_p and trv_lyndon_p * function indirectly uses trv_i_euler_p */ { trv_cycle_p(A) = A = trv_i_moebius_p(trv_lyndon_p(A)); A[2] = 1; A } /* inverse cycle product transform * this function ignores a_0, a_1 (A[0], A[1]) and assumes them to be 0, 1 * the output sequence always has b_0=b_1=0 * function uses trv_moebius_p and trv_i_lyndon_p * function indirectly uses trv_euler_p */ { trv_i_cycle_p(A) = trv_i_lyndon_p(trv_moebius_p(A)); } /* cycle labeled transform * apply cycle structure to sequence of labeled objects * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 * function uses trv_log */ { trv_cycle_l(A) = A=-trv_log(-A); A[1]=1; A } /* inverse cycle labeled transform * this function ignores a_0 (A[1]) and assumes it to be 1 * the output sequence always has b_0=0 * function uses trv_exp */ { trv_i_cycle_l(A) = A=-trv_exp(-A); A[1]=0; A } /* cycle "boy-girl" transform * apply cycle structure to sequence of unlabeled objects * insisting that neighbors be different object types * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=1 * function uses trv_carlitz, trv_i_euler and trv_i_moebius */ { trv_cycle_bg(A) = local(B = dirmul(vecextract(A,"2.."),vector(#A-1,n,(n==1)-(n==2)))); A + trv_i_moebius(trv_i_euler(trv_carlitz(A)) - concat(0,B)) } /* euler transform for matrices * this function ignores a_0,0 (A[1,1]) and assumes it to be 0 * the output sequence always has b_0,0=1 */ { trm_euler(A) = local(x='x,y='y,m=matsize(A)[1],n=matsize(A)[2],B=O(x^m)+O(y^n), j=1,k,ki,kj); for(i=0, m-1, for(j=j, n-1, k=1; ki=i; kj=j; while(ki1, M[j,i]+=M[j-1,i-1]); ); ); M } /* gauss (sequence) transform * sequence of row sums of previous transform * outputs a sequence 1 longer than input */ { trv_gauss(A) = local(M=trv_gauss_m(A)); vector(#A+1, i, sum(j=1, i, M[i,j])) } /* inverse gauss (sequence) transform * It treats expects the first value of the input vector * to be 1. (It treats it that way and ignores the value) * outputs a sequence 1 shorter than input */ { trv_i_gauss(A) = local(M=matid(#A),B=vector(#A-1)); for(i=1, #A-1, B[i] = M[i+1,i] = A[i+1] - sum(j=1, i-1, M[i+1,j]) - 1; if(i>1, B[i]-=M[i,i-1]); for(j=i+2, #A, M[j,i] = B[i]*M[j-1,i]; if(i>1, M[j,i]+=M[j-1,i-1]); ); ); B } /* hankel transform */ { trv_hankel(A) = local(B=vector(1+#A\2)); for(i=1, #B-1, B[i+1]=matdet(matrix(i,i,x,y,A[x+y-1])); ); B } /* little hankel transform */ { trv_hankel_little(A) = local(B=vector(#A-1)); B[1]=0; for(i=2, #A-1, B[i]=A[i]^2-A[i+1]*A[i-1]; ); B } /* invert transform for matrices * this function ignores a_0,0 (A[1,1]) and assumes it to be 0 * the output sequence always has b_0,0=1 */ { trm_invert(A) = local(x='x,y='y,m=matsize(A)[1],n=matsize(A)[2],B=O(x^m)+O(y^n),j=1); for(i=0, m-1, for(j=j, n-1, B+=A[i+1,j+1]*x^i*y^j); j=0; ); B = 1/(1-B); for(i=0, m-1, for(j=0, n-1, A[i+1,j+1] = polcoeff(polcoeff(B,j,y),i,x)); ); A } /* inverse invert transform for matrices * this function ignores a_0,0 (A[1,1]) and assumes it to be 1 * the output sequence always has b_0,0=0 */ { trm_i_invert(A) = local(x='x,y='y,m=matsize(A)[1],n=matsize(A)[2],B=1+O(x^m)+O(y^n),j=1); for(i=0, m-1, for(j=j, n-1, B+=A[i+1,j+1]*x^i*y^j); j=0; ); B = 1-1/B; for(i=0, m-1, for(j=0, n-1, A[i+1,j+1] = polcoeff(polcoeff(B,j,y),i,x)); ); A } /* extended version of invert transform * it does not assume a_0 = A[1] == 0 */ { trv_invert_x(A) = local(x='x, B=Ser(A,x)); B = 1/(1-B); for(i=1, #A, A[i] = polcoeff(B,i-1)); A } /* lah transform * apply given labeled structure to non-empty sets */ { trv_lah(A) = local(x='x,B); B = sum(i=0, #A-1, A[i+1]*x^i/i!); B = subst(B+O(x^#A),x,x/(1-x+O(x^#A))); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* inverse lah transform */ { trv_i_lah(A) = local(x='x,B); B = sum(i=0, #A-1, A[i+1]*x^i/i!); B = subst(B+O(x^#A),x,x/(1+x+O(x^#A))); for(i=0,#A-1, A[i+1] = polcoeff(B,i)*i!); A } /* lyndon transform for matrices * this function ignores a_0,0 (A[1,1]) and assumes it to be 0 * the output sequence always has b_0,0=1 * function uses trm_i_euler */ { trm_lyndon(A) = A = -trm_i_euler(-A); A[1,1] = 1; A } /* inverse lyndon transform for matrices * this function ignores a_0,0 (A[1,1]) and assumes it to be 1 * the output sequence always has b_0,0=0 * function uses trm_euler */ { trm_i_lyndon(A) = A = -trm_euler(-A); A[1,1] = 0; A } /* lyndon product transform * apply lyndon word structure (or aperiodic cycles or * necklaces) to sequence of unlabeled objects * to assemble them by the product of their sizes * this function ignores a_0, a_1 (A[0], A[1]) and assumes them to be 0 * the output sequence always has b_0=0, b_1=1 * function uses trv_i_euler_p */ { trv_lyndon_p(A) = A = -trv_i_euler_p(-A); A[2] = 1; A } /* inverse lyndon product transform * this function ignores a_0, a_1 (A[0], A[1]) and assumes them to be 0, 1 * the output sequence always has b_0=b_1=0 * function uses trv_euler_p */ { trv_i_lyndon_p(A) = A = -trv_euler_p(-A); A[2] = 0; A } /* mask transform * From 6/9/2001 seqfan posting by Marc LeBrun */ { trv_mask(A) = forstep(i=#A-1,0,-1,A[i+1]=sum(j=0,i,(bitand(i,j)==j)*A[j+1])); A } /* inverse mask transform */ { trv_i_mask(A) = for(i=1,#A-1,for(j=0,i-1,A[i+1]-=(bitand(i,j)==j)*A[j+1])); A } /* moebius transform for matrices * this function ignores a_0,0 (A[1,1]) and assumes it to be 0 * the output sequence always has b_0,0=0 */ { trm_moebius(A) = local(m=matsize(A)[1],n=matsize(A)[2],j=1,ki,kj); for(i=0, (m-1)\2, for(j=j, (n-1)\2, ki=i+i; kj=j+j; while(ki1, A[j]-=C[j]/k); forstep(m=j+j, #A, j, C[m]-=C[j]); k++; v*=A[i]; ); ); concat(0,A) } /* inverse witt transform * slice of lyndon transform and related to somos transform * this function ignores a_0 (A[1]) and assumes it to be 0 * the output sequence always has b_0=0 */ { trv_i_witt(A) = local(B=vector(#A-1),C=vector(#B),k,v); for(i=1, #B, forstep(j=i, #B, i, C[j]=0); v=A[i+1]; k=1; forstep(j=i, #B, i, C[j]+=v; B[j]+=C[j]/k; forstep(m=j+j, #B, j, C[m]-=C[j]); k++; v*=A[i+1]; ); ); concat(0,B) } /* and convolution operation * this operation is not reversible */ { opv_conv_and(A, B) = local(C=vector(min(#A,#B))); forstep(i=#C-1,0,-1,C[i+1]=sum(j=0,i,bitand(A[j+1],B[i-j+1]))); C } /* gcd convolution operation * this operation is not reversible */ { opv_conv_gcd(A, B) = local(C=vector(min(#A,#B))); forstep(i=#C-1,0,-1,C[i+1]=sum(j=0,i,gcd(A[j+1],B[i-j+1]))); C } /* lcm convolution operation * this operation is not reversible */ { opv_conv_lcm(A, B) = local(C=vector(min(#A,#B))); forstep(i=#C-1,0,-1,C[i+1]=sum(j=0,i,lcm(A[j+1],B[i-j+1]))); C } /* or convolution operation * this operation is not reversible */ { opv_conv_or(A, B) = local(C=vector(min(#A,#B))); forstep(i=#C-1,0,-1,C[i+1]=sum(j=0,i,bitor(A[j+1],B[i-j+1]))); C } /* xor convolution operation */ { opv_conv_xor(A, B) = local(C=vector(min(#A,#B))); forstep(i=#C-1,0,-1,C[i+1]=sum(j=0,i,bitxor(A[j+1],B[i-j+1]))); C } /* xor division operation */ { opv_div_xor(A, B) = local(C=vector(min(#A,#B))); for(i=0,#C-1,C[i+1]=bitxor(A[1],B[i+1])-sum(j=0,i-1,bitxor(A[i-j+1],C[j+1]))); C } /* cycle "boy-girl" operation * apply cycle structure to sequence B of unlabeled objects * while sequence A governs when neighbors can be the same object types * sequence length is the shorter of the 2 * this function assumes a_0 (A[1]) == 1 and b_0 (B[1]) == 0 * ignoring their actual values * the output sequence always has c_0=1 * function uses opv_invert_bg, trv_i_euler and trv_i_moebius */ { opv_cycle_bg(A, B) = local(C = trv_i_euler(A)); C = dirmul(vecextract(B,"2.."),vecextract(C,"2..")); C = trv_i_moebius(trv_i_euler(opv_invert_bg(A,B)) - concat(0,C)); C[1] = 1; C + concat(0,dirmul(vecextract(A,"2.."),vecextract(B,"2.."))); } /* euler "type" operation * With A as "base" and B as "exponent" * apply set structure to sequence of unlabeled objects (B) allowing * a_n cases where there are n of the same type of object * sequence length is the shorter of the 2 * this function assumes a_0 (A[1]) == 1 and b_0 (B[1]) == 0 * ignoring their actual values * the output sequence always has c_0=1 * function uses trv_euler and trv_i_euler */ { opv_euler_type(A, B) = A = trv_i_euler(A); A = dirmul(vecextract(A,"2.."),vecextract(B,"2..")); trv_euler(concat(0,A)) } /* exponential "type" operation * With A as "base" and B as "exponent" * apply set structure to sequence of labeled objects (B) allowing * a_n cases where there are n of the same type of object * (as an L-species) * sequence length is the shorter of the 2 * this function assumes a_0 (A[1]) == 1 and b_0 (B[1]) == 0 * ignoring their actual values * the output sequence always has c_0=1 * function uses trv_exp, trv_log and opv_mul_ed */ { opv_exp_type(A, B) = trv_exp(opv_mul_ed(trv_log(A),B)) } /* invert "boy-girl" operation * apply list structure to sequence B of unlabeled objects * while sequence A governs when neighbors can be the same object types * sequence length is the shorter of the 2 * this function assumes a_0 (A[1]) == 1 and b_0 (B[1]) == 0 * ignoring their actual values * the output sequence always has c_0=1 * function uses trv_invert and trv_i_invert */ { opv_invert_bg(A, B) = A = trv_i_invert(A); A = dirmul(vecextract(A,"2.."),vecextract(B,"2..")); trv_invert(concat(0,A)) }