\\---------------------------------------------------------------------- \\ GP scripts \\ \\ Experimental Number Theory \\ Fernando Rodriguez Villegas \\ Oxford University Press 2007 \\ \\---------------------------------------------------------------------- recognize(x)=local(m);m=contfracpnqn(contfrac(x)); m[1,1]/m[2,1] recognize2(x)=bestappr(x, 10^(precision(x) - 5)) count(f,p)=sum(a=0,p-1,subst(f,x,Mod(a,p))==0) bqfn(a,b)= ceil(1/2+b/shift(a,1))-1 bqfdisc(q)=q[2]^2-shift(q[1]*q[3],2) bqfeval(q,v)=q[1]*v[1]^2+q[2]*v[1]*v[2]+q[3]*v[2]^2 floorsqrt(x)=sqrtint(floor(x)) memb(g,v)=for(k=1,length(v),if(g==v[k],return(k)));0 sum2squares2(n)=4*sumdiv(n,d,kronecker(-4,d)) sum4squares2(n)=8*sumdiv(n,d,if(d%4,d)) sershift(u)=(u-polcoeff(u,0))/x partnum2(n) = Vec(1/etaps(n)-1) partnum3(n) = Vec(1/eta(x + O(x^(n+1))) - 1) harmp1(n,p,m)=sum(k=1,n,1/k,O(p^m)) logeps(d,p,m = 10) = log(quadunit(d) + O(p^m)) polgraeffer(f,r)=polresultant(subst(f,x,t),x-t^r,t) polcyclomax(n)=vecmax(abs(Vec(polcyclo(n)))) \\---------------------------------------------------------------------- testmod(v,bd)= { local(n,mn,S); mn=length(v); S=[0,mn]; for(D=-bd,bd, if(!isfundamental(D),next); n=dist(v,D); if(n < mn, mn=n; S=[D,mn])); S } dist(v,D)= { sum(j=1,length(v), abs(v[j][2]-kronecker(D,v[j][1]))) } sgn(f, p)= { if (poldisc(f)%p, return(0)); f = factormod(f,p)[,1]; (-1)^sum(j=1,length(f), poldegree(f[j]) + 1) } sgn1(p)= { if(6%p == 0, return(0)); centerlift(Mod(1,p)*sum(n=0,p-1, (6*n)!/n!/(2*n)!/(3*n)!/432^n)) } sgn2(p) = { local(M, S); if (6%p == 0, return (0)); S = M = Mod(1,p); for (n=1, p-1, S += (M *= (6*n-1)*(6*n-5) / (n^2*36))); centerlift(S) } gosper(v)= { local(a,b,j,l,c); a=contfrac(v[1]); b=contfrac(v[2]); j=1; l=min(length(a),length(b)); while(j <= l && a[j] == b[j],j++); c=vector(j-1,k,a[k]); if(j <= l, c=concat(c,if(a[j] < b[j],a[j],b[j])+1)); c=contfracpnqn(c); c[1,1]/c[2,1] } recognize1(x,n)= { local(e); if(n==0,n=precision(x)); e=.5/10^n; gosper([x-e,x+e]) } recognizemod(c,m)= { local(z,N); N=floor(sqrtint(floor(m/2))); z=bestappr(c/m,N); c-m*z } ellunit(D)= { local(z); z=(1+sqrt(D))/2; norm(eta(z/2,1)/eta(z,1))/sqrt(2) } bernoulli(n,x)= { local(h,s,c); h=0;s=0;c=-1; for(k=1,n+1, c*=1-(n+2)/k; s+=x^n; x++; h+=c*s/k); h } antidiff(f)= { local(x,h,s,y,c,n); x=variable(f); h=0;s=0;c=-1;y=x; n=poldegree(f); f=intformal(f); for(k=1,n+1, c*=1-(n+2)/k; s+=subst(f,x,y); y++; h+=c*s/k); h } polsum(f,a,b)= { local(h,s,c,n); h=0;s=0;c=-1;b++; n=poldegree(f); f=intformal(f); for(k=1,n+1, c*=1-(n+2)/k; s+=subst(f,x,b)-subst(f,x,a); a++;b++; h+=c*s/k); h } sumsquares(n, k) = { local(bd = sqrtint(n)); forvec(x = vector(k,i, [-bd,bd]), if (n == sum(i=1, k, x[i]^2), print(x))) } sum2squares(n) = { local(x1, x2, v); v=[]; if (n<0 || n%4 == 3, return([])); for (x1=0, floorsqrt(n/2), if (issquare(n - x1^2, &x2), v=concat(v,[[x1,x2]]))); v } sum2squares1(n) = { local(x1, x2); if (n<0 || n%4 == 3, return([])); for (x1=0, floorsqrt(n/2), if (issquare(n - x1^2, &x2), return([x1,x2]))); } sum3squares1(n)= { local(x_1,x_2,x_3,n_1); if(n<0 || n%8==7,return([])); for(x_1=0,floorsqrt(n/3), n_1=n-x_1^2; for(x_2=x_1,floorsqrt(n_1/2), if(issquare(n_1-x_2^2,&x_3), return([x_1,x_2,x_3])))) } sum4squares1(n)= { local(x_1,x_2,x_3,x_4,n_1,n_2); if(n<0,return([])); for(x_1=0,floorsqrt(n/4), n_1=n-x_1^2; for(x_2=x_1,floorsqrt(n_1/3), n_2=n_1-x_2^2; for(x_3=x_2,floorsqrt(n_2/2), if(issquare(n_2-x_3^2,&x_4), return([x_1,x_2,x_3,x_4]))))) } sum2squares3(n)= { local(x_1,x_2,v); ct=0; if(n<0 || n%4==3,return(0)); for(x_1=0,floorsqrt(n/2), if(issquare(n-x_1^2,&x_2), ct+=((x_1>0)+1)*((x_2>0)+1)*((x_1 2,, e = 3); bd=10^e; e=bd/100; forprime(p=2,bd, if(polfacttype(f,p)==t,a=a++); b++; if(b%e,,print(a/b*1.))); print(a/b*1.) } poltypedensity1(f,t,N,e)= { local(a,b,bd); if(e > 2,,e = 3); bd=10^e; e=bd/100; forprime(p=2,bd, if(polfacttype(f,p)==t,a=a++); b++; if(b%e,,print1(bestappr(a/b,N)," "))); } poltypedensity2(f,galv,t,N,e)= { local(a,b,bd); if(e>2,,e=3);bd=10^e;e=bd/100; forprime(p=2,bd, if(polallfrob(p,galv,f)==t,a=a++); b++; if(b%e,,print1(bestappr(a/b,N)," "))); } polisgalois(f, bd = 100 * poldegree(f))= { local(v,k); forprime(p=2,bd, v = polfacttype(f,p); for (k=2, length(v), if (v[k] != v[k-1], return(0)))); 1 } polgal8(bd)= { local(v,f); forvec(u = concat([[[1,1], [0,bd]], vector(6,k,[-bd,bd]), [[0,0]]]), f=Pol(u); gal8test(f+1); gal8test(f-1)) } gal8test(f)= { if(polisirreducible(f) && polisgalois(f), print(f," ",factor(nfdisc(f)))) } galisabelian(galv,f)= { local(j,k,c); for(j=1,length(galv)-1, for(k=j+1,length(galv), c=subst(galv[j],x,galv[k])-subst(galv[k],x,galv[j]); c=c%f; if(c,return(0)))); 1 } galorder(g,f)= { local(k,gg); k=1;gg=g; while(k<1000 && gg!=x, gg=subst(gg,x,g)%f;k++); k } bqfred(q)= { local(top,bot, a=q[1], b=q[2], c=q[3]); top=[-b,1]; bot=[shift(a,1),0]; while(1, applyTpow(); if (c >= a, break); applyS()); if(a == c && b > 0, applyS()); bot[1] = (bot[1]-b*bot[2])/(2*a); [[a,b,c], bot] } applyTpow() = { local(aux, n = ceil(1/2+b/shift(a,1))-1); if (!n, return); aux=a*n; c+=(aux-b)*n; b-=shift(aux,1); top+=n*bot; } applyS() = { local(aux); aux=c; c=a; a=aux; b=-b; aux=top; top=-bot; bot=aux; } fermat1(p)= { local(b); b=lift(sqrt(Mod(-1,p))); b=bqfred([p,2*b,(b^2+1)/p])[2]; vecsort(abs(b)) } checkfermat1(n)= { local(p,v); until(p%4==1,p=nextprime(n++)); v=fermat1(p); [p,v,p-(v[1]^2+v[2]^2)] } bqf(D)= { local(b0 = D%2, fv,a,c,zv); if(D >= 0 || D%4 > 1, return([])); fv = [[1,b0,(b0^2-D)/4]]; forstep(b = b0, floorsqrt(-D/3), 2, zv=divisors((b^2-D)/4); n=length(zv); forstep(j=(n+1)\2, 2, -1, a=zv[j]; if (a < b, break); c=zv[n-j+1]; if(gcd(gcd(a,b),c) != 1, next); fv = concat(fv,[[a,b,c]]); if(b && a != b && a != c, fv = concat(fv,[[a,-b,c]])))); fv } bqfold(D)= { local(fv,a,b0,j,c,z,zv); if(D < 0 && (D%4 == 0 || D%4 == 1), b0=D%2; fv=[[1,b0,(b0^2-D)/4]]; forstep(b = b0,floorsqrt(-D/3),2, z=(b^2-D)/4; zv=divisors(z); n=length(zv); j=(n+1)\2; a=zv[j]; c=zv[n-j+1]; while(a >= b && j > 1, if(gcd(gcd(a,b),c)==1, fv=concat(fv,[[a,b,c]]); if(b != 0 && a != b && a != c, fv=concat(fv,[[a,-b,c]]))); j--; a=zv[j]; c=zv[n-j+1])); fv, []) } checkbqf(x)= { for(d=x,x+1000, v=bqf(-d);if(v,if(length(v)-qfbclassno(-d),print(-d)))) } bqfcomp(q1,q2)= { local(aux,m,a,b,c,dd,n); dd=-bqfdisc(q1); if(bqfdisc(q2) != -dd,return([])); aux=bezout(q1[1],q2[1]); m=aux[2]; aux=bezout(aux[3],(q1[2]+q2[2])\2); m=m*aux[1];n=aux[2]; a=q2[1]\aux[3]; b=q2[2]+(m*(q1[2]-q2[2])\2-n*q2[3])*a*2; a=a*q1[1]\aux[3]; b=b%(2*a); c=(b*b+dd)\(4*a); bqfsred([a,b,c]) } bqfpow(q,n)= { local(t,dd); if(n == 0, dd=-bqfdisc(q); return([1,dd%2,(dd+dd%2)/4])); if(n == 1, return (q)); t = bqfpow(bqfcomp(q,q), n\2); if(n % 2, bqfcomp(q,t), t); } bqfprimef(p,D)= { local(b = centerlift(sqrt(Mod(D,p)))); if((b-D)%2,b=b+p); [p,b,(b^2-D)/4/p] } checkcubic(bd)= { forprime(p=2, bd, if(kronecker(-31,p)!=1, next); if((polfacttype(x^3+x+1,p)==[1,1,1]) != bqfistrivial(p,-31), print(p))) } bqfistrivial(p,D)= { local(q); q=bqfprimef(p,D); if(bqfsred(q)[1]==1,1,0) } bqftheta(Q,bd)= { local(b); b=sqrtint(bd); sum(n=-b,b, sum(m=-b,b, q^bqfeval(Q,[m,n]))+O(q^bd)) } trinomial(n)= { local(pol = 1, pol0 = 1+x+x^2); concat(1,vector(n, k, polcoeff(pol*=pol0, k))); } trinomial1(n)= { local(c = 1); 1 + sum(k=0, n\2-1, c *= (n-2*k)*(n-2*k-1)/(k+1)^2); } serlindep(v,l,flag)= { local(M); M=matrix(l,length(v),j,k,polcoeff(v[k],j-1)); M=if(flag, qflll(M,4)[1], matker(M)); if(M,M[,1]) } serdiffeq(u,r,d,l,flag = 1)= { local(w,S); if(l, u+=O(x^l), l=length(u)); if(l < (d+1)*(r+1), error("Not enough precision in series.")); w=vector(r+1); w[1]=u; for(k=1,r, w[k+1]=x*deriv(w[k],x)); v=w; for(j=1,d, v = concat(v, w*=x)); S=serlindep(v,l,flag); if(S, vector(r+1,k, sum(j=0,d, S[j*(r+1)+k]*x^j))); } seqlinrec(a,r,d,l,flag = 1)= { local(w,v,S); if(!l, l=length(a)); if(l < (d+1)*(r+1)+d, error("Not enough terms of the sequence.")); w=vector(d+1); w[1]=Polrev(a); for(j=1,d, w[j+1] = w[j]\x); v=w; for(k=1,r, v = concat(v, w=x*deriv(w,x))); S=serlindep(v,l-d,flag); if(S, S=vector(d+1,j, sum(k=0,r, S[k*(d+1)+j]*x^k))); } trinomial2(n)= { local(v); v=vector(n+1); v[1] = v[2] = 1; for(k=2,n, v[k+1] = ((2*k-1)*v[k]+3*(k-1)*v[k-1]) / k); v } trinomial3(n)= { local(a, a0 = 0, a1 = 1); for(k=1,n, a = ((2*k-1)*a1+3*(k-1)*a0) / k; a0=a1; a1=a); a1 } recsolve(N,Q,a)= { local(w,q, d = length(Q)-1); w = concat(a, vector(N+d+1 - length(a))); for(n=0,N, q = subst(Q,x,n); w[n+d+1] = -sum(k=0,d-1, q[k+1]*w[n+k+1])/q[d+1]); w } seralgdep(u,r,d,l,flag = 1)= { local(w,v); if(l, u+=O(x^l), l=length(u)); if(l < (d+1)*(r+1), error("Not enough precision on input.")); w=vector(d+1); w[1]=1; for(k=1,d, w[k+1] = u*w[k]); v=vector((d+1)*(r+1),k, w[(k-1)\(r+1)+1] * x^((k-1)%(r+1))); v=serlindep(v,l,flag); if(v, v/=content(v); v=vector(d+1,k, sum(j=0,r,v[(k-1)*(r+1)+j+1]*x^j)); sum(k=0,d,v[k+1]*y^k)); } hensel(f,w,n)= { local(g,d,c); d=poldegree(w); c=-1/truncate(subst(subst(deriv(f,u),x,0),u,w+O(x))); for(k=d+1,d+n, g=subst(f,u,w+O(x^(k+1))); g=polcoeff(g,k); w+=c*g*x^k); w+O(x^(d+n+1)) } hensel1(f,w,n)= { local(g,d,c,r); d=poldegree(w); c=subst(subst(deriv(f,u),x,0),u,w+O(x^(d+1))); r=valuation(c,x); c=-1/polcoeff(c,r); for(k=d+1,d+n, g=subst(f,u,w+O(x^(k+1))); g=polcoeff(g,k+r); w+=c*g*x^k); w+O(x^(d+n+1)) } newtsqrt(D, n = 5, a = 1, b = 1)= { local(aux); for(k=1,n, aux=a^2+D*b^2; b*=2*a; a=aux); a/b } newtsqrt1(D, n = 5, a = 1)= { local(m); m=1; for(k=1,n, m = shift(m,1); a = truncate((D+O(x^m))/a + a)/2); a+O(x^m) } sercontfrac(z)= { local(v,c); c=truncate(z+O(x)); v=[c]; z-=c; while(z, z=1/z; c=truncate(z+O(x)); z-=c; v=concat(v,c)); v } contfracper(ci,cp, m = 1)= { local(v,v0,v1); v0=[0,1];v1=[1,0]; for(n=1,length(ci), v=ci[n]*v1+v0; v0=v1;v1=v); for(k=1,m, for(n=1,length(cp), v=cp[n]*v1+v0; v0=v1;v1=v)); v1[1]/v1[2] } trinasympt(n)= { local(A,B,g,c); c=1; A=(-2+x)/sqrt(1-x+O(x^(n+2))); B=(-1+x)/sqrt(1-2*x+O(x^(n+2))); for(k=1,n, g=c+z*x^k+O(x^(k+2)); g=3*g+A*subst(g,x,x/(1-x))+B*subst(g,x,x/(1-2*x)); g=polcoeff(g,valuation(g,x)); c=c-polcoeff(g,0)/polcoeff(g,1)*x^k); c+O(x^(n+1)) } sumasympt(C,z)= { local(w,n); w=vector(length(C),k, abs(polcoeff(C,k-1)) * z^(k-1)); n=vecsort(w,,1)[1]; [ subst(truncate(C+O(x^n)),x,z), n ] } trinapprox(n)= { local(s); s=sumasympt(C,1./n); [1/sqrt(n)*sqrt(3/4/Pi)*s[1],s[2]] } seqispol(v,d)= { local(p,l); p=v[1]; l=length(v); for (k = 1, d+1, for(j=1,l-1, v[j]=v[j+1]-v[j]); v[l] = 0; if (!v, return(p)); l--; p += binomial(x,k)*v[1]); } seqisratnl(v,d, flag = 1)= { local(A); if(flag, A=concat(matrix(length(v),d+1,j,k,numerator(v[j])*(j-1)^(k-1)), matrix(length(v),d+1,j,k,-denominator(v[j])*(j-1)^(k-1))); A=qflll(A,4)[1], A=concat(matrix(length(v),d+1,j,k,v[j]*(j-1)^(k-1)), matrix(length(v),d+1,j,k,-(j-1)^(k-1))); A=matker(A)); if(A, A=A[,1]; sum(n=0,d,A[n+d+2]*x^n)/sum(n=0,d,A[n+1]*x^n), 0) } serrecognize(f)= { local(m); m=contfracpnqn(sercontfrac(f)); m[1,1]/m[2,1] } contfracrq(D)= { local(v); v = contfrac(sqrt(D)); serrecognize(Ser(v) + O(x^(length(v)-1)) ) } somos6(n)= { local(a); a=vector(n+1); for(k=1,6,a[k]=1); for(k=6,n, a[k+1]=(a[k]*a[k-4]+a[k-1]*a[k-3]+a[k-2]^2)/a[k-5]); a } lineq(mv,m,flag)= { local(k,j,sm,sj,s,sv, S = [], N = length(mv)); k = j = 1; sm = sj = vector(m\mv[1]+1); while(k, /* We are done when stack is empty */ /* Create new object s,j */ s = sm[k]+mv[j]; if (s > m, /* Object too large, backtrack */ until(j <= N, j = sj[k]+1; k--); /* k-- = pop last */ next); /* Object suitable, push it at the top of the stack */ k++; sm[k]=s; sj[k]=j; /* here s <= m, j <= N */ if (s < m, next); /* Not a solution to linear equation, go on */ /* Actual solution, add to solution vector */ if (flag, sv = vector(k-1,l, mv[sj[k-l+1]]) , /* else */ sv = vector(N); for(l=2,k, sv[sj[l]]++)); S = concat(S,[sv])); S } part(m)= { local(k,j,sm,sj,s, S = []); k = j = 1; sm = sj = vector(m+1); while(k, s = sm[k]+j; if (s > m, until(j <= m, j = sj[k]+1; k--); next); k++; sm[k]=s; sj[k]=j; if (s == m, S = concat(S, [vector(k-1,l, sj[k-l+1])]))); S } partnum(m)= { local(k,j,sm,sj,s, S = vector(m)); k = j = 1; sm = sj = vector(m+1); while(k, s = sm[k]+j; if (s > m, until(j <= m, j = sj[k]+1; k--); next); k++; sm[k]=s; sj[k]=j; S[s]++); S } Psi(n, q) = { local(a, b, c); a=sqrt(2/3)*Pi/q; b=n-1/24; c=sqrt(b); (sqrt(q)/(2*sqrt(2)*b*Pi))*(a*cosh(a*c)-(sinh(a*c)/c)) } G(h, q) = if(q<3,0,sum(k=1,q-1,k*(frac(h*k/q)-1/2))) L(n, q) = { if(q==1,1, sum(h=1,q-1, if(gcd(h,q)>1,0,cos((G(h,q)-2*h*n)*Pi/q)))) } partnum4(n) = round(sum(q=1,5 + 0.24*sqrt(n),L(n,q)*Psi(n,q))) etaps(n)= { local(z,t1,t2,q1,q2,q3); z=1+O(x^(n+1));t1=1;t2=1; q1=x; q2=x^2; q3=x^3; until(poldegree(t2)>n, t1*=-q1;t2*=-q2; z+=t1+t2; q1*=q3;q2*=q3); z } partdual(v)= { local(w); w=vector(v[1]); for(j=1,length(v), for(i=1,v[j], w[i]++)); w } partdual1(v) = { local(last, w = vector(v[1])); w[1] = last = length(v); for (i=2, v[1], while (v[last] < i, last--); w[i] = last); w; } partdim(v)= { local(w,ct,dim); w=partdual(v); dim=1;ct=1; for(j=1,v[1], for(k=1,w[j], dim=dim*ct/(v[k]+w[j]-j-k+1); ct++)); dim } partz(v)= { local(w); w=vector(v[1]); for(j=1,length(v), w[v[j]]++); prod(j=1,length(w), j^w[j]*w[j]!) } \\ Remove rim hook of length m from partition v. The input is vr=[v,r] \\ where r is a number mod 2, which keeps track of the sign for the \\ recursive calculation of character values. partrmrim(vr,m)= { local(Solutionv,newv,v,neww,w,row,rimlen,col,r,currlen); v=vr[1]; w=partdual(v); Solutionv=[]; r=0; newv=v; /* Col represents the column where the current rim begins Row represents the row where the current rim ends */ col=1; row=w[col]; /* Start off having removed the last row */ newv[row] = 0; rimlen = v[row]; if(rimlen == m, Solutionv=concat(Solutionv,[[newv,vr[2]]])); /* could say "if rim==m, Sol=this, Sol=[]" and take out init of Soln */ while(row > 1 || (rimlen >= m && v[row] != col), if ((rimlen < m) || (v[row] == col), /* Move up a row to increase the length if the length is too small, or if moving over a column will delete the rim */ row--; newv[row] = v[row+1]-1; rimlen = rimlen + v[row]-newv[row], /* else: move over a column to decrease the length */ neww=partdual(newv); if (length(neww) 0, neww[col] = w[col], neww=concat(neww,w[col])); newv = partdual(neww); col++); if(rimlen == m, Solutionv=concat(Solutionv,[[newv,(vr[2]+w[col]-row)%2]]))); Solutionv } partcharcheck(n)= { local(w,v,M); w=part(n); v=vector(length(w),k,partz(w[k])); M=matrix(length(w),length(w),j,k,partchar(w[j],w[k])); matid(length(w))==M*matdiagonal(vector(length(v),k,1/v[k]))*M~ } partchar(v,u)= { local(V,W); V=[[v,0]]; for(i=1,length(u), W=[]; for(j=1,length(V), W=concat(W,partrmrim(V[j],u[i]))); V=W); sum(k=1,length(W),if(W[k][2],-1,1)) } partchartable(n)= { local(w); w=part(n); matrix(length(w),length(w),j,k,partchar(w[j],w[k])) } partnosol(ccv)= { local(pv,m,n,nf); n=sum(k=1,length(ccv[1]),ccv[1][k]); pv=part(n);m=length(ccv);nf=n!; prod(k=1,m,nf/partz(ccv[k]))/nf* sum(j=1,length(pv), prod(k=1,m, partchar(pv[j],ccv[k]))/ partdim(pv[j])^(m-2)) } parttglhom(rv,n)= { /* Run a copy of part to get all partitions w of m <= n. */ local(k,j,sm,sj,s,ct,w); k = j = 1; sm = sj = vector(n+1); /* ct is the counting polynomial; set the intial value to account for the trivial homomomorphism into S_k. */ ct=sum(k=0,n,x^k)+O(x^(n+1)); while(k, s = sm[k]+j; if (s > n, until(j <= n, j = sj[k]+1; k--); next); k++; sm[k]=s; sj[k]=j; w=vector(k-1,i,sj[k-i+1]); /* With partition w compute the contribution to the total ct from non-trivial conjugacy classes of orders rv[i] i=1,2,... */ ct+=prod(i=1,length(rv), partpsum(rv[i],w))/partdim(w)*x^s/s!); ct } partpsum(i,v)= { local(s,k,n,z); /* n varies: n, n-i, n-2i,... with |v| = n as initial value; n!/z is the size of the corresponding conjugacy class C = (i,...,i,1,...1) with k i's and (current value of) n 1's */ n=sum(k=1,length(v),v[k]); k=0; z=1; s=0; while(i <= n, k++; z*=prod(j=1,i,n-j+1)/i/k; n-=i; s+=z*partchar(v,concat(vector(k,j,i),vector(n,j,1)))); s } parttglsubgps(rv,n)= { /* Run a copy of part to get all partitions w of m <= n. */ local(k,j,sm,sj,s,ct,w); k = j = 1; sm = sj = vector(n+1); /* ct is the counting polynomial; set the intial value to account for the trivial homomomorphism into S_k. */ ct=sum(k=0,n,x^k/k!)+O(x^(n+1)); while(k, s = sm[k]+j; if (s > n, until(j <= n, j = sj[k]+1; k--); next); k++; sm[k]=s; sj[k]=j; w=vector(k-1,i,sj[k-i+1]); /* With partition w compute the contribution to the total ct from non-trivial conjugacy classes of orders rv[i] i=1,2,... */ ct+=prod(i=1,length(rv), partpsum(rv[i],w))/partdim(w)*x^s/s!^2); x*deriv(ct,x)/ct } eulerphilist(n)= { local(p,pv,N,M,s,j,k,m,sph,sj,sm); pv=[];p=1;a=1; forprime(p=2,n+1,pv=concat(pv,p)); M=ceil(log(n)/log(2))+2; S=vector(n\2,k,[]); N=length(pv); k=1;j=1; sph = vector(M,l,1); sm=sph; sj=vector(M); while(k, /* We are done when stack is empty */ /* Create new object s,m,j */ s=sph[k]*(pv[j]-if(j==sj[k],0,1)); m=sm[k]*pv[j]; if(s > n, /* Object too large, backtrack */ until(j <= N,j = sj[k]+1; k--); /* k-- = pop last */ next); /* Object suitable, push it at the top of the stack */ k++; sph[k] = s; sm[k] = m; sj[k] = j; /* here s <= n, j <= N */ if(s > 1,S[s/2]=concat(S[s/2],m))); vector(length(S),k,vecsort(S[k])) } polcyclolist(n)= { local(nuv,nuv1,spol,sj,sdeg,N,j,k,pol,S,s); if(n > 1, nuv1 = eulerphilist(n)); phv=[1,1];nuv=[1,2]; S=[]; for(k=1,length(nuv1), for(j=1,length(nuv1[k]), phv=concat(phv,2*k); nuv=concat(nuv,nuv1[k][j]))); k=1; j=1; N=length(phv); spol = vector(N,k,1); sdeg = vector(N); sj = sdeg; while(k, /* Create new object s, pol, j */ s = sdeg[k]+phv[j]; pol = spol[k]*polcyclo(nuv[j]); if (s > n, /* Object too large, backtrack */ until(j <= N, j = sj[k]+1; k--); /* k-- = pop last */ next); /* Object suitable, push it at the top of the stack */ k++; sdeg[k] = s; sj[k] = j; spol[k] = pol; /* here s <= n, j <= N */ if (s < n, next); /* Not a solution to linear equation, go on */ S=concat(S,pol)); S } gammap(x)= { local(n,a,p,m,d,t,s); p=component(x,1); m=padicprec(x,p); n=ceil(p*m/(p-1)); a=-x%p; x=(x+a)/p;a++; x=truncate(x)+O(p^(m+floor(n/(p-1)))); d=vector(p); d[1]=1+O(p^m);t=1; for(j=1,p-1,d[j+1]=d[j]/j); s=d[a]; for(k=1,n-1, d[1]=d[1]+d[p]; for(j=1,p-1,d[j+1]=(d[j]+p*k*d[j+1])/(k*p+j)); t*=(x+k-1)/k; s+=d[a]*t); s } gammapv(x)= { local(n,a,p,m,d,t,v); p=component(x,1); m=padicprec(x,p); n=ceil(p*m/(p-1)); x=truncate(x)+O(p^(m+floor(n/(p-1)))); d=vector(p);v=d; d[1]=1+O(p^m);v[1]=d[1];t=1; for(j=1,p-1,d[j+1]=d[j]/j); for(j=2,p,v[j]=d[j]); for(k=1,n-1, t*=(x+k-1)/k; d[1]=d[1]+d[p]; v[1]+=d[1]*t; for(j=1,p-1,d[j+1]=(d[j]+p*k*d[j+1])/(k*p+j)); for(j=2,p,v[j]+=d[j]*t)); v } clogp(x,p)= { local(k); x=abs(x); k=0; until(x==0,x=x\p;k++); k } psip(x)= { local(n,a,p,pp,m,np,d,t,t1,s,s1); p=component(x,1); m=padicprec(x,p); n=ceil(p*m/(p-1)); np=clogp(n,p); while(n-floor(n/p)-np < m, n++; np=clogp(n,p)); m+=np; a=-x%p; x=(x+a)/p;a++; x=truncate(x)+O(p^(m+floor(n/(p-1)))); d=vector(p); d[1]=1+O(p^m);t=1;t1=0; for(j=1,p-1,d[j+1]=d[j]/j); s=d[a];s1=0; for(k=1,n-1, d[1]=d[1]+d[p]; for(j=1,p-1,d[j+1]=(d[j]+p*k*d[j+1])/(k*p+j)); t1*=(x+k-1);t1+=t;t1/=k; t*=(x+k-1)/k; s+=d[a]*t; s1+=d[a]*t1); s1/s/p } harmp(n, p, m = 20)= { local(psip1, s = 0,pp = 1); m+=clogp(n,p); psip1=psip(1+O(p^m)); while(n, s += (psip(n+1+O(p^m)) - psip1) / pp; n\=p; pp*=p); s } diam(d,p, m = 10)= { -sum(r=1,d-1, kronecker(d,r)*psip(r/d+O(p^m))) /sqrt(d+O(p^m))/(1-1/p)/2 } diamtest(p,Dmax,m)= { for(d=2,Dmax, if(issquarefree(d), D=quaddisc(d); if(D 1 && m < M && polisirreducible(f), print(f, "\t", m))); } polmahler(f)= { local(rv); rv=abs(polroots(f)); polcoeff(f,poldegree(f))* prod(k=1,length(rv),if(rv[k]>1,rv[k],1)) } polgraeffe(f)= { local(v,d,hd); d=poldegree(f); hd=d\2; v=[sum(k=0,hd,polcoeff(f,2*k)*x^k), sum(k=0,hd,polcoeff(f,2*k+1)*x^k)]; f=v[1]^2-x*v[2]^2; if(d%2,-f,f) } poliscyclo(f)= { local(f1,k,n); f=f/x^valuation(f,x); n=poldegree(f); if(abs(polcoeff(f,0)*polcoeff(f,n)) > 1, 0, k=1; n=shift(n,2); until(f1 == f || k > n, f1=f; f=polgraeffe(f1); k=shift(k,1)); if(k > n,0,1)) } polcyclofactor(f)= { local(g,h,k,p,l,n); g=1; k=1; l=1; f=f/x^valuation(f,x); until(l > n , p=prime(k);l*=p-1;k++; h=gcd(f,polgraeffer(f,p)); if(poldegree(h),f/=h;g*=h); n=poldegree(f)); g } polwedge(f,l)= { local(n,hn,v,ff,m,fv,a0); n = poldegree(f); a0=polcoeff(f,0); m=n; a0 = if(n%2,-a0,a0); if(l, ,l=n); hn=min(l,n\2); v = concat(f,vector(binomial(n,hn)-1,r,polgraeffer(f,r+1))); fv=vector(l+1); fv[1]=x-1;fv[2]=f; if(l >= n,fv[n+1]=x-a0); for(k=2,min(l,n-1), m*=(n-k+1)/k; ff=exp((-1)^(k-1)*sum(r=1,m, polcoeff(v[r],n-k)*x^r/r,O(x^(m+1)))); ff=truncate(ff); ff=ff/polcoeff(ff,0); ff=subst(ff,x,1/x)*x^m; fv[k+1]=ff); fv } interlacing(A,B)= { local(w,k,s, r = length(A)); if(r != length(B), error("Sets must be of the same size")); w=concat(vector(r,k,[A[k],0]), vector(r,k,[B[k],1])); w=vecsort(w,1); for (k = 2, length(w), if (w[k][1]==w[k-1][1] || w[k][2]==w[k-1][2], return(0))); 1; } polcycloexp(v)= { local(u,w,m); w=[];m=0; for(k=1,length(v), if(v[k] > m, m=v[k]; u=[]; if(v[k]==1,u=[0], for(n=1,v[k]-1, if(gcd(n,v[k])==1,u=concat(u,n/v[k]))))); w=concat(w,u)); w } ssasymptview(k, bd = 100)= { local(s,th,vk,wk,ck); th=sum(n=1,sqrtint(bd),2*x^(n^2),1+O(x^(bd+1))); vk=Vec(th^k-1); ck=Pi^(k/2)/gamma(k/2+1); wk=vector(bd); s=1; for(n=1,bd, s+=vk[n]; wk[n]=s/ck/n^(k/2)-1); ploth(x=1,bd,wk[round(x)]) } poltypesf(f,a = 2, b = 1000)= { local(v,pv); w=[];pv=[]; forprime(p=a,b, v=polfacttype(f,p); if(v==[] || memb(v,w),, w=concat(w,[v]); pv=concat(pv,p))); [pv,w] } idoneus(bd)= { local(hv,h,j,S); S=[]; for(d=1,bd, if(d%4==0 || d%4==3, h = qfbclassno(-d); if(h == 1, S=concat(S,[[-d,1]]), if(twopower(h), hv=bqf(-d); j=2; while(bqfpow(hv[j],2) == hv[1] && j < length(hv),j++); if(j == length(hv), S=concat(S,[[-d,j]])))))); S }} twopower(n)= { local(z,r); until(r, z=divrem(n,2); n=z[1]; r=z[2]); (n==0) } bhfred(H)= { local(a,b,c,top,bot,n,aux); a=H[1];b=H[2];c=H[3]; top=[-b,1]; bot=[shift(a,1),0]; n=bqfn(a,real(b))+I*bqfn(a,imag(b)); if(n, c+=a*norm(n)-shift(trace(n*conj(b)),-1); b-=shift(a*n,1); top+=n*bot); while(c < a, aux=c; c=a; a=aux; b=-conj(b); aux=top; top=-bot; bot=aux; n=bqfn(a,real(b))+I*bqfn(a,imag(b)); if(n, c+=a*norm(n)-shift(trace(n*conj(b)),-1); b-=shift(a*n,1); top+=n*bot)); if((imag(b) < 0) || (imag(b) == 0 && real(b) < 0), b=-b; top*=I; bot*=-I); if((a == c && real(b) < 0), b=-conj(b); aux=top; top=-bot; bot=aux); H=[a,b,c]; bot[1]=(bot[1]-H[2]*bot[2])/2/H[1]; [H,bot] } sum2sqmod(p)= { local(u,v); u=0; v=-1; while(kronecker(v,p) != 1, v-=shift(u,1)+1; u++); [u,lift(sqrt(Mod(v,p)))] } sum4sq(p)= { local(H,v); v=sum2sqmod(p); H=bhfred([p,shift(v[1]+v[2]*I,1),(v[1]^2+v[2]^2+1)/p])[2]; vecsort(abs([real(H[1]),imag(H[1]),real(H[2]),imag(H[2])])) } boydcheck(v,N)= { local(a,aux); a0=8;a1=55; for(n=3,N, aux=floor(a1^2/a0+1); a0=a1; a1=aux; if(v[n]!=a1,print(n))); } cantorcheck(v,N)= { local(a,aux); a0=3;a1=10; for(n=3,N-1, aux=floor(a1^2/a0+1/2); a0=a1; a1=aux; if(v[n+1]!=a1,print(n))); } randfib(m,b)= { local(f0,f,s,v); v=vector(m+1); s=0;v[1]=1;f0=1; for(n=1,m, f=1+b*(-1)^random(2)/f0; s=s+log(abs(f)); v[n+1]=s/(n+1);f0=f); v } syracuse(a)= { local(v); v=[a]; while(a>1,a=(1+a)*a-floor(a/2)*(2*a+1);v=concat(v,a)); v } nondecseq(n,a,b)= { local(k,j,sj,S); S=[]; k=1;j=a; sj=vector(n+1); while(k, if(k > n, until(j <= b, j=sj[k]+1; sj[k]=0; k--); next); k++; sj[k]=j; S=concat(S,[vector(k-1,l,sj[l+1])])); S } landaux(n)= { local(v); if(n<0,0, n++;v=vector(n,i,1); for(j=2,n, for(i=1,j-1, v[j]=max(v[j],lcm(i,v[j-i])))); v) } ser2prod(f,m)= { if(m==0,m=length(f)+1); f=f+O(x^m); if(subst(truncate(f),x,0)!=1, error("Series must be of the form 1 + ..."), f=-log(f); sum(n=1,length(f), -sumdiv(n,d,polcoeff(f,d)*moebius(n/d)*d)/n*x^n)) } \\----------------------------------------------------------------------