/* two versions of guessgf, the latter much better */
guessgf(v,B=10)=local(l,m,p,q);l=length(v);m=matrix(l,B,x,y,if(x<y,0,v[x-y+1]));q=qflll(m);p=sum(k=1,B,x^(k-1)*q[k,1]);q=Pol(Pol(vector(l,n,v[l-n+1]))*p+O(x^(B+1)));q/p
ggf(v)=local(l,m,p,q,B);l=length(v);B=floor(l/2);if(B<4,return(0));m=matrix(B,B,x,y,v[x-y+B+1]);q=qflll(m,4)[1];if(length(q)==0,return(0));p=sum(k=1,B,x^(k-1)*q[k,1]);q=Pol(Pol(vector(l,n,v[l-n+1]))*p+O(x^(B+1)));if(polcoeff(p,0)<0,q=-q;p=-p);q=q/p;p=Ser(q+O(x^(l+1)));for(m=1,l,if(polcoeff(p,m-1)!=v[m],return(0)));q

/* Zeckendorf representation, Fibonomial coefficients */
Zec(n)=local(k,v,m);k=0;while(fibonacci(k)<=n,k=k+1);m=k-1;v=vector(m-1);v[1]=1;n=n-fibonacci(k-1);while(n>0,k=0;while(fibonacci(k)<=n,k=k+1);v[m-k+2]=1;n=n-fibonacci(k-1));v
lucas(n)=if(n<1,2*(n>=0),fibonacci(n-1)+fibonacci(n+1))
FF(k,n)=if(k<0||k<n,0,prod(i=k-n+1,k,fibonacci(i))/prod(i=1,n,fibonacci(i)))
PP(k,n)=if(k<0||k<n,0,prod(i=k-n+1,k,pell(i))/prod(i=1,n,pell(i)))
pell(n)=if(n<2,n>0,2*pell(n-1)+pell(n-2))

/* cyclic resultant */
CRes(pol,n)=polresultant(pol,x^n-1,x)
/* find numbers n such that n^2,(n+1)^2,...,(n+k-1)^2 is square */

A000792(n) = floor(3^(n-4-(n-4)\3*2)*2^(-n%3))
A001317(n) = local(k,s);s=0;for(k=0,n,s=2*s+(binomial(n,k)%2));s
A002516(n) = local(t);if(n<1,0,if(n%2==0,2*A002516(n/2),t=(n-1)/2;3*t+1/2-(t-5/2)*(-1)^t))
A003239(n)=if(n<1,n>=0,sumdiv(n,k,eulerphi(n/k)*binomial(2*k,k))/(2*n))
A003458(n) = local(m,i);m=0;i=n+1;while(m<=n,i=i+1;m=factor(binomial(i,n))[1,1]);i
A004718(n) = if(n<1,0, if(n%2==0, -A004718(n/2), A004718((n-1)/2)+1))
A007601(n) = local(m,t);m=n;for(p=0,2,t=2*p+3*ceil(log3(n/2^p));if(t<m,m=t));m
A025480(n)=if(n%2==0,n/2,A025480((n-1)/2))

A072894(n,N=1000)=local(v,s);v=vector(N);v[1]=1;v[2]=n;for(k=3,N,s=v[k-1]+v[k-2];v[k]=if(s%2,(s+1)/2,s/2);if(v[k]==v[k-1],return(v[k])));v[N]
fillA086450(m)=A086450=vector(m);A086450[1]=1;s=1;for(n=2,m,A086450[n]=if(n%2==0,A086450[n/2],s=s+A086450[(n+1)/2];s))
A081882(n)=local(c,s);if(n<0,0,s=n-1/4;c=0;while(frac(s)>0,s=s*ceil(s);c++);c)
A083368(n)=local(l1,l2,z1,z2);if(n==1,return(1));z1=Zec(n-1);z2=Zec(n);l1=length(z1);l2=length(z2);if(l2>l1,return(l1));for(k=1,l1,if(z1[k]<>z2[k],return(l1-k)));-1

xD(n, pol) = for(k=1,n,pol=x*deriv(pol,x));pol
ntothe(e) = xD(e,1/(1-x))

PisotE(a,b,n)=if(n<2,if(n<1,if(n<0,0,a),b),floor(PisotE(a,b,n-1)^2/PisotE(a,b,n-2)+1/2))
PisotL(a,b,n)=if(n<2,if(n<1,if(n<0,0,a),b),ceil(PisotL(a,b,n-1)^2/PisotL(a,b,n-2)))
PisotP(a,b,n)=if(n<2,if(n<1,if(n<0,0,a),b),ceil(PisotP(a,b,n-1)^2/PisotP(a,b,n-2)-1/2))
PisotT(a,b,n)=if(n<2,if(n<1,if(n<0,0,a),b),floor(PisotT(a,b,n-1)^2/PisotT(a,b,n-2)))

/* no 3-term arithmetic progression, two methods */
NAPold(sv,N)=local(v,f,m,n,k,l,sl);sl=length(sv);v=vector(N);for(k=1,sl,v[k]=sv[k]);for(n=sl+1,N,f=2;m=v[n-1]+1;while(1,forstep(k=n-1,1,-1,if(v[k]<(m+1)/2,f=1;break);for(l=1,k-1,if(m-v[k]==v[k]-v[l],f=0;break));if(f<2,break));if(!f,m=m+1;f=2);if(f==1,break));v[n]=m);v
NAP(sv,N)=local(v,vv,m,k,l,sl,vvl);sl=length(sv);vvl=min(N*N,10^5);v=vector(N);vv=vector(vvl);for(k=1,sl,v[k]=sv[k];for(l=1,k-1,vv[2*v[k]-v[l]]=1));m=v[sl]+1;for(k=sl+1,N,while(m<=vvl&&vv[m],m=m+1);if(m>vvl,return(v));for(l=1,k-1,sl=2*m-v[l];if(sl<=vvl,vv[sl]=1));vv[m]=1;v[k]=m);v

/* triple sum free sequence starting with sv, and max value m */
trisumfree(sv,m)=local(nv,v,vl,nvl,c,t);vl=length(sv);nv=vector(max(m,3*sv[vl]));v=listcreate(m+5);for(k=1,vl,listput(v,sv[k]));for(i=1,vl,for(j=1,vl,for(k=1,vl,t=sv[i]+sv[j]+sv[k];if(t,nv[t]=1))));nvl=3*sv[vl];c=sv[vl]+1;while(1,while(c<=m&&nv[c],c++);if(c>m,break);vv=Vec(v);listput(v,c);for(k=1,vl,for(l=1,vl,t=c+vv[k]+vv[l];if(t<=m,nv[t]=1)));for(k=1,vl,t=c+c+vv[k];if(t<=m,nv[t]=1));if(3*c<=m,nv[3*c]=1);nvl=3*c;c=c+1);Vec(v)

/* number of unique elements in n X n multiplication table */
multab(N)=local(v,cv,s,t);v=vector(N);cv=vector(N*N);v[1]=1;cv[1]=1;for(k=2,N,s=0;for(l=1,k,t=k*l;if(cv[t]==0,s++);cv[t]++);v[k]=v[k-1]+s);v

S1(n,k)=if(n<1,0,n!*polcoeff(binomial(x,n),k))
S2(n,k)=(1/k!)*sum(i=0,k,(-1)^(k-i)*binomial(k,i)*i^n)
Bell(n)=sum(k=1,n,S2(n,k))
polyB(k,n)=if(n>=0,(-1)^n*sum(m=0,n,(-1)^m*m!*S2(n,m)/(m+1)^k),0)
BernB(n,x)=sum(i=0,n,binomial(n,i)*bernfrac(i)*x^(n-i))
Balg(v)=sum(m=0,length(v)-1,(-1)^m*m!*S2(n,m)*v[m+1])
Eu(n)=sum(m=0,n,(-1)^m*m!*S2(n+1,m+1)*(-1)^floor(m/4)*2^-floor(m/2)*((m+1)%4!=0))
EulerE(n,v='x)=2/(n+1)*(BernB(n+1,x)-2^(n+1)*BernB(n+1,x/2))
Eulerian(n,k)=sum(j=0,k,(-1)^j*(k-j)^n*binomial(n+1,j))
trinomial(n,k)=polcoeff((1+x+x^2)^n,k)
Zsig(AA,B,n)=local(d,f,D);d=divisors(AA^n-B^n);forstep(D=length(d),1,-1,f=0;for(m=1,n-1,if(gcd(d[D],AA^m-B^m)>1,f=1;break));if(!f,return(d[D])));1

bdA(n)=1/n!*sum(k=1,n,(-1)^(n-k)*S1(n,k)*A(k))
buA(n)=sum(k=1,n,(-1)^(n-k)*S2(n,k)*k!*A(k))

listdcgfA(al,pol,m)=local(l);for(n=0,m,l=if(!n,0,ceil(log2(n)));print1(polcoeff(sum(k=0,l,al^k*subst(pol,x,x^2^k))+O(x^(2^l+1)),n)","))
listdcgfB(spol,al,pol,m)=local(l);for(n=0,m,l=if(!n,0,ceil(log2(n)));print1(polcoeff(1/(1-x)*(spol+sum(k=0,l,al^k*subst(pol,x,x^2^k))),n)","))
listdcgfC(spol,al,pol,m)=local(l);for(n=0,m,l=if(!n,0,ceil(log2(n)));print1(polcoeff(1/(1-x)^2*(spol+sum(k=0,l,al^k*subst(pol,x,x^2^k))),n)","))
filldcgfA(al,pol,m)=local(l,v);v=vector(m+1);for(n=0,m,l=if(!n,0,ceil(log2(n)));v[n+1]=polcoeff(sum(k=0,l,al^k*subst(pol,x,x^2^k)),n));v
filldcgfB(spol,al,pol,m)=local(l,v);v=vector(m+1);for(n=0,m,l=if(!n,0,ceil(log2(n)));v[n+1]=polcoeff(1/(1-x)*(spol+sum(k=0,l,al^k*subst(pol,x,x^2^k))),n));v
filldcgfC(spol,al,pol,m)=local(l,v);v=vector(m+1);for(n=0,m,l=if(!n,0,ceil(log2(n)));v[n+1]=polcoeff(1/(1-x)^2*(spol+sum(k=0,l,al^k*subst(pol,x,x^2^k))),n));v
listdcgfA4(al,pol,m)=local(l);for(n=0,m,l=if(!n,0,ceil(log4(n)));print1(polcoeff(sum(k=0,l,al^k*subst(pol,x,x^4^k)),n)","))
listdcgfG(c,d,pol,m)=local(g,l);l=if(m<=0,0,ceil(log2(m)));g=pol+sum(k=0,l,subst(pol,x,x^(2^(k+1)))*prod(kk=0,k,c+d*x^(2^kk)));for(n=0,m,print1(polcoeff(g,n)","))
listdcgfQ(c,d,pol,m)=local(g,l);l=if(m<=0,0,ceil(log2(m)));g=sum(k=0,l,subst(pol,x,x^(2^(k)))*prod(kk=0,k,c+d*x^(2^kk)));for(n=0,m,print1(polcoeff(g,n)","))
/* ternary morphism */
termorph(n, m) = local(s1, s2); s1=[0];for(n=1,10,s2=vector(2*#s1);for(k=1,#s1,s2[2*k-1]=m[s1[k]+1,1];s2[2*k]=m[s1[k]+1,2]);s1=s2);s2[n+1]

fillA035928(m)=local(l,b,c);c=1;A035928=vector(m);for(n=1,8*4^floor(log2(m)),if(c<=m,b=binary(n);l=length(b);if(sum(i=1,l,abs(component(b,i)-component(b,l+1-i)))==l,A035928[c]=n;c=c+1;print(c)),break))

\\ partly borrowed from Michael Somos' tricks.gp
derivn(F,n,v='x)=while(n>0,F=deriv(F,v);n-=n);F
/* Chebyshev polynomial function ( T_n(x) ) of first kind */
ChebT(n,v='x)=subst(poltchebi(n),'x,v)
/* Chebyshev polynomial function ( U_n(x) ) of second kind */
ChebU(n,v='x)=if(n==-1,0,subst(poltchebi(n+1)'/(n+1),'x,v))
ChebS(n,v='x)=if(n==-1,0,subst(poltchebi(n+1)'/(n+1),'x,v/2))
/* Laguerre polynomial function ( L_n(x) ) */
LagL(n,v='x)=sum(k=0,n,(-1)^k*n!/(k!^2*(n-k)!)*v^k)
LagL3(n,al,v='x)=Pol(1/n!*exp(v)*v^-al*derivn(exp(-v)*v^(n+al),n,v))
LagG(n,al,v='x)=sum(k=0,n,binomial(n+al,n-k)*(-v)^k/k!)
/* Hermite polynomial function ( H_n(x) ) */
HermH(n,v='x)=sum(k=0,n\2,(-1)^k*n!/(k!*(n-2*k)!)*(v*2)^(n-2*k))
Hermh(n,v='x)=round((-I/sqrt(2))^n*HermH(n,I*x/sqrt(2)))
HermHH(n,mu,v='x)=(-1)^(n\2)*2^n*(n\2)!*LagG(n\2,if(n%2==0,mu-1/2,mu+1/2),v*v)
GegC(n,lam,x)=sum(k=0,n\2,(-1)^k*lam^(n-k)*(2*x)^(n-2*k)/k!/(n-2*k)!)
tprs(x)=local(t);t=type(x);(t=="t_POL"|t=="t_RFRAC"|t=="t_SER")
Bessel_J(n,x='x)=if(!tprs(x),besselj(n,x),if(n>=0,(x/2)^n/n!*besselj(n,x),(-1)^n*Bessel_J(-n,x)))

trunc(x)=if(x>=0,floor(x),ceil(x))
log2(x)=log(x)/log(2)
log3(x)=log(x)/log(3)
log4(x)=log(x)/log(4)
flog2(x)=local(c);c=0;while(2^(c+1)<x+1,c=c+1);c /* integers only */
clog2(x)=local(c);c=0;while(2^c<x,c=c+1);c
v2(n)=if(n==0,0,valuation(n,2))
lod(n)=n/(2^v2(n))
lod2(n)=(lod(n)-1)/2
e0(n)=if(n==0,0,sum(k=0,floor(log2(n)),!bittest(n,k)))
e1(n)=if(n==0,0,sum(k=0,floor(log2(n)),bittest(n,k)))
u1(c,d,n)=if(n<1,0,if(n%2==0,c*u1(c,d,n/2),d*u1(c,d,(n-1)/2)+1))
u2(c,d,n)=if(n<1,0,if(n%2==0,c*u2(c,d,n/2)+1,d*u2(c,d,(n-1)/2)))

ternary(n)=local(v,t,l);l=floor(log3(n+0.1));v=vector(l+1);forstep(k=l,0,-1,t=floor((n+0.1)/3^k);v[l-k+1]=t;n=n-t*3^k);v
engel(x,n)=(-1)^(n+1)*engelu(x,n)
engelu(x,n)=if(n<=1,floor(1/x),trunc(1/(x-sum(k=1,n-1,1/engelu(k))))) 
addmod2(x,y)=local(bx,by,lx,ly);bx=binary(x);by=binary(y);if(x<y,bx=by;by=binary(x));lx=length(bx);ly=length(by);for(k=1,ly,bx[lx-k+1]=(bx[lx-k+1]!=by[ly-k+1]));sum(k=1,lx,bx[k]*2^(lx-k))

/*A076215=[0 ,1 ,77 ,92 ,70 ,195 ,143 ,3854 ,357 ,245 ,413 ,4088 ,2257 ,2222 ,652 ,679 ,278949 ,3366 ,1281 ,67963 ,1612 ,8555 ,1518 ,63412 ,1159158 ,2619 ,2725 ,13862 ,60973 ,3069 ,10790 ,3128 ,4620 ,5083 ,42918 ,3406 ,34430 ,96045 ,4017 ,214099 ,8756 ,10412 ,37829 ,7990 ,199160 ,19912 ,8140 ,1455222 ,214638 ,358387 ,358387 ,28167 ,180104 ,41847 ,259739 ,45793 ,12443 ,1597754 ,10857 ,8033 ,17750 ,4019959 ,8075 ,79068 ,9010 ,15975 ,10384 ,12787 ,3280695 ,18497 ,14345 ,16817 ,184601560 ,14548378 ,43273 ,61248 ,264718 ,16332 ,369349 ,26610 ,28106 ,20802705 ,119277 ,17267 ,1634193 ,328322 ,4480969 ,441227 ,33659 ,326190 ,25185 ,35211 ,45318 ,244306 ,257412 ,34034 ,92240 ,70943 ,55930 ,242198 ,940383 ,445341 ,64960 ,227423 ,32500 ,84932 ,48604 ,29575 ,47377 ,28026001 ,1178332 ,2267337 ,76738 ,1009866]

A001032=[1,2,11,23,24,26,33,47,49,50,59,73,74,88,96,97,107,121,122,146,169,177,184,191,193,194,218,239,241,242,249,289,297,299,311,312,313,337,338,347,352,361,362,376,383,393,407,409,431,443,458,479,481,491,506,529,537,539,554,568,577,578,599,600,625,649,673,674,698,722,753,767,793,794,841,856,864,866,887,897,913,914,961,971,983,1009,1019,1041,1048,1058,1067,1079,1081,1129,1144,1153,1163,1175,1199,1201,1202,1225,1226,1248,1249,1257,1274,1298,1307,1321,1346,1369,1391]

A007475=[1,3,18,7,1,25,7,539,25,7,22,442,225,192,13,15,26914,244,50,5552,30,553,7,4493,83342,83,65,775,3807,64,556,20,106,132,2277,15,1788,5063,27,11320,280,358,1805,210,9985,802,183,71752,10123,16806,1081,7989,1663,11476,1778,255,68680,172,25,443,167065,3,2927,25,301,38,117,126031,322,131,197,6665193,516232,1122,1678,8617,65,12116,413,454,688013,3480,3,51958,9977,140563,13310,479,9548,183,500,805,6884,7090,379,2120,1472,1009,6387,26533,12240,1210,5873,223,1752,694,106,612,774564,31758,61127,1352,26379]

fillA081504(m)=local(c,f,t);A081504=vector(m);c=1;for(n=2,100000,f=0;for(i=1,n-1,t=2^n+2^i+1;if(isprime(t),f=1;break));if(!f,print1(n",");A[c]=n;c=c+1;if(c>=m,break)))

weight(q)=local(l);l=length(q~);sum(k=1,l,if(q[k,][1],q[k,][1]^2,-1000))
mw=10^11;v=vector(20,n,Lucas(n-1)^3);oo=0
/*for(k=-5,5,for(l=-10,10,if(k!=l,m=matrix(length(v)-oo,5,x,y,if(y==1,(-1)^x,if(y==2,1,if(y==3,if(2*x+k<0,0,fibonacci(2*x+k)),if(y==4,fibonacci(if(2*x+l<0,0,2*x+l)),v[x])))));q=qflll(m,4)[1];if(length(q)&&weight(q)<=mw,mw=weight(q);print(k","l"; "q)))))
for(k=-4,3,print(k);for(l=-10,10,if(k!=l,for(kk=-10,10,for(ll=-10,10,if(kk!=ll,m=matrix(length(v)-oo,6,x,y,if(y==1,if(1*x+kk<0,0,(-1)^x*fibonacci(1*x+kk)),if(y==2,if(1*x+ll<0,0,(-1)^(x)*fibonacci(1*x+ll)),if(y==3,if(3*x+k<0,0,fibonacci(3*x+k)),if(y==4,fibonacci(if(3*x+l<0,0,3*x+l)),if(y==5,1,v[x]))))));q=qflll(m,4)[1];if(length(q)&&weight(q)<=mw,mw=weight(q);print(kk","ll","k","l"; "q)))))))) */