\\ Given the Igusa invariants this will try to find a small curve. igusa2curve(ivm) = { local(pol,p,dum,fmt); pol = sectix(ivm); if(pol == 0, print("The curve can not be defined over the rationals!"), if(pol == 1,, print("Found the big curve!"); fmt = default(format,,1); default(format,"g0.2"); pol = wangsmall(pol); default(format,fmt); if(ivm == igusainvpol(pol), pol, print("Oops: the igusa invariant changed!")))) } \\ Given the igusa invariants and a zero of mestre's conic, this gives \\ the sectix. sex(imv,lst) = { local(dum,tlst,mesmat); mesmat = mestmat(imv); dum = [1,x,t*(x-lst[2]/lst[1])+lst[3]/lst[1]]; dum = dum*mesmat*dum~/(x-lst[2]/lst[1]); tlst = [component(dum,2),-component(dum,1),t*(-component(dum,1)-lst[2]/lst[1]*component(dum,2))+lst[3]/lst[1]*component(dum,2)]; thecubic(imv,tlst) } \\ Given the Igusa invariants this gives a big sectix sectix(imv) = { local(dum); dum = mestratq(imv); if(type(dum) == "t_INT",dum,subst(hypred(clrden(sex(imv,dum[2]))),t,x)) } clrden(pol) = { local(dum); dum = denominator(component(pol,1)); for(k=2,7,dum = lcm(dum,denominator(component(pol,k)))); pol*dum } hypred(pol) = { local(dum); dum = component(pol,1); for(k=2,7,dum = gcd(dum,component(pol,k))); pol = pol/dum } thecubic(iv,sx) = { local(i200,i020,i002,i012,i013,i022,i030,i031,i040,i102,i103,i111,i112,i120,i121,i130,i201,i202,i210,i211,i220,i300,i301,i310,i400,im41,im50); i200 = iv[1]*iv[1]; i020 = iv[2]*iv[2]; i002 = iv[3]*iv[3]; i300 = i200*iv[1]; i210 = i200*iv[2]; i201 = i200*iv[3]; i400 = i300*iv[1]; i310 = i300*iv[2]; i301 = i300*iv[3]; i220 = i210*iv[2]; i211 = i210*iv[3]; i202 = i201*iv[3]; i102 = i002*iv[1]; i103 = i102*iv[3]; i120 = i020*iv[1]; i130 = i120*iv[2]; i121 = i120*iv[3]; i012 = i002*iv[2]; i013 = i012*iv[3]; i022 = i012*iv[2]; i030 = i020*iv[2]; i031 = i030*iv[3]; i040 = i030*iv[2]; i112 = i012*iv[1]; i111 = iv[1]*iv[2]*iv[3]; im41 = i040*iv[3]/iv[1]; im50 = i040*iv[2]/iv[1]; ((-1792*im41)/9 + (6464*im50)/135 - (8960*i013)/3 + (18752*i022)/9 - 38400*i030 - (12640*i031)/27 + (14144*i040)/405 - 640000*i102 + (352*i103)/3 + (800000*i111)/3 - (448*i112)/9 - (269440*i120)/9 + (64*i121)/15 + (496*i130)/3375 + (20800*i201)/3 - (472*i202)/225 - (12352*i210)/9 + (916*i211)/1125 - (428*i220)/5625 - (656*i300)/45 + (148*i301)/16875 - (452*i310)/253125 - i400/93750)*sx[3]^3 + sx[1]*(((-3375*i211)/16 + (2925*i220)/64 - (253125*i300)/8 + (12825*i301)/1024 - (1125*i310)/512 - (243*i400)/8192)*sx[1]^2 + sx[2]*(((-75*i130)/4 - (3375*i202)/4 + 225*i211 - (1155*i220)/128 - (253125*i300)/16 + (1485*i301)/128 - (1089*i310)/512 - (243*i400)/10240)*sx[1] + ((75*i130)/8 - (3375*i202)/8 - 33750*i210 + (2475*i211)/16 - (465*i220)/32 - (3375*i300)/2 + (225*i301)/64 - (213*i310)/320 - (81*i400)/12800)*sx[2] + (-40*i040 + 1800*i112 - 780*i121 + (331*i130)/4 + 540000*i201 - (1215*i202)/4 - 105750*i210 + (435*i211)/4 - (46*i220)/5 - (5175*i300)/2 + (21*i301)/10 - (163*i310)/400 - (27*i400)/8000)* sx[3]) + sx[3]*(((75*i130)/8 - (3375*i202)/8 - 33750*i210 + (2475*i211)/16 - (465*i220)/32 - (3375*i300)/2 + (225*i301)/64 - (213*i310)/320 - (81*i400)/12800)*sx[1] + ((160*i031)/3 + (26*i040)/9 + 2400*i103 - 1410*i112 - 36000*i120 + 276*i121 - (833*i130)/45 + 45000*i201 - (119*i202)/2 - 9600*i210 + (1411*i211)/60 - (2101*i220)/900 - 165*i300 + (31*i301)/100 - (23*i310)/375 - (9*i400)/20000)*sx[3])) + sx[2]*((-20*i040 - 300*i112 + 180*i121 - (68*i130)/3 - 45000*i201 - (105*i202)/2 + 6000*i210 + (85*i211)/4 - (71*i220)/30 + 75*i300 + (7*i301)/20 - (163*i310)/2400 - (9*i400)/16000)*sx[2]^2 + sx[3]*(((-1520*i031)/3 + (1148*i040)/9 - 1200*i103 + 1380*i112 - 108000*i120 - 420*i121 + (1717*i130)/45 + 180000*i201 - 49*i202 - 34800*i210 + (259*i211)/15 - (319*i220)/225 - 570*i300 + (31*i301)/100 - (23*i310)/375 - (9*i400)/20000)*sx[2] + ((-64*im50)/3 - 2880*i022 + (4544*i031)/3 - (8824*i040)/45 - 1200000*i111 + 312*i112 + 304000*i120 - (708*i121)/5 + (140*i130)/9 - 129600000*i200 + 84000*i201 - (86*i202)/5 - 16400*i210 + (461*i211)/75 - (557*i220)/1125 - 212*i300 + (34*i301)/375 - (41*i310)/2250 - (3*i400)/25000)*sx[3])) } \\ Given a symmetric matrix qmat, this will try to find a matrix P such that \\ P*Q*P~ is diagonal. \\ It uses the Gram-Schmidt orthogonalization process and gives a message if \\ it finds a vector such that v*Q*v~ = 0 diagsymm(qmat) = { local(dum,len,zro,bs); len = length(qmat); bs = vector(len,k,vector(len,j,k==j)); dum = vector(len,k,0); zro = 0; for(k=1,len, bs[k] = bs[k] - sum(j=1,k-1,(bs[k]*qmat*bs[j]~/dum[j])*bs[j]); if((dum[k] = bs[k]*qmat*bs[k]~) == 0,if(zro == 0,print("v = ",bs[k]," gives v*Q*v~ = 0");zro = bs[k],);dum[k]=1,)); if(zro == 0, matrix(len,len,k,l,bs[k][l]), zro) } \\ This gives the matrix for mestre's conic. Input are the igusa invariants \\ as given above. mestmat(imv) = { local(mx,my,mz); mx = 8/225*(20*imv[2]+imv[1])/imv[1]; my = 16/3375*(-600*imv[3]+imv[1]+80*imv[2])/imv[1]; mz = -64/253125*(-10800000*imv[1]-9*imv[1]^2-700*imv[2]*imv[1]+3600*imv[3]*imv[1]+12400*imv[2]^2-48000*imv[2]*imv[3])/imv[1]^2; [mx+6*my, 6*mx^2+2*my, 2*mz; 6*mx^2+2*my, 2*mz , 9*mx^3+4*mx*my+6*my^2; 2*mz, 9*mx^3+4*mx*my+6*my^2, 6*mx^2*my+2*my^2+3*mx*mz] } \\ Give this the igusa invariants mestratq(imv) = { local(dum,dum2,a,b,am,bm,mesmat,dmat); mesmat = mestmat(imv); dmat = diagsymm(mesmat); if(type(dmat) != "t_MAT",[mesmat,dmat], dum = dmat*mesmat*dmat~; dum = vector(3,k,factor(dum[k,k])); print("done factoring!"); dum2 = fmult(dum[3],frec(dum[1])); dum = fmult(dum[2],frec(dum[1])); am = fdenom(dum); bm = fdenom(dum2); dum = fmult(fmult(dum,am),am); dum2 = fmult(fmult(dum2,bm),bm); a = fsqtpart(dum); b = fsqtpart(dum2); am = f2r(fmult(am,frec(a))); bm = f2r(fmult(bm,frec(b))); a = fmult([-1,1],fmult(dum,frec(fmult(a,a)))); a = [f2r(a),a]; b = fmult([-1,1],fmult(dum2,frec(fmult(b,b)))); b = [f2r(b),b]; if(abs(a[1]) 10^50) & (type(lst[dum[1]][2]) == "t_INT"), fnum = fnum+1; if(factlist[fnum] == 0, print("don't wanna factor ",dum2); return(1), if(dum2 == f2r(factlist[fnum]), nb = [dum2,factlist[fnum]], print("factlist is not correct!");floor(i))), nb = [dum2,factor(dum2)]); dum = lst[dum[1]][3]; b = fsqrfrprt(nb); nb = f2r(fsqtpart(nb[2]))); if(abs(a[1]) 5, for(l=1,dum[k,2]\30, pol = redpdsc(dum[k,1],pol)),)); pol = nsmpol(pol); if(sizeofpol(pol) < 100, pol = smallhyp(pol,11,20),); pol) } nsmpol(pol) = { local(bsz,done,dum,dum2,inv,dir); bsz = sizeofpol(pol); dir = 1; done = 0; inv = 0; until(done, if(inv,print("low size ",bsz),); if((dum = sizeofpol(dum2 = hypred(hypsub(pol,1,round(dir),0,1)))) < bsz, pol = dum2;inv=0;bsz = dum; dir = dir*1.4, if(abs(dir)>1.5,dir = dir/2;inv=0, if((dum=sizeofpol(dum2=hypred(hypsub(pol,1,-sign(dir),0,1)))) < bsz, pol = dum2;inv=0;bsz = dum;dir=-1.4*sign(dir), pol = hypred(hypsub(pol,0,1,1,0)); dir = 1; if(inv,done=1,inv=1))))); pol } sizeofpol(pol) = log(sum(k=1,7,component(pol,k)^2)/7)/(log(10)*2) hypsub(pol,a,b,c,d) = { pol = subst(pol,x,x/y)*y^6; subst(subst(pol,x,a*x+b),y,c*x+d) } \\ This attempts to factor the discriminant of a sectix giving a simple \\ hyperelliptic curve. facdisc(dsc) = { local(dum,pow,dum2,pres); dum2 = factor(dsc,0); if(dum2[length(dum2~),1] < 10^10,[1,dum2,1], dsc = dum2[length(dum2~),1]; pow = 30/dum2[length(dum2~),2]; dum2 = matrix(length(dum2~)-1,2,k,l,dum2[k,l]); pres = default(realprecision,,1); default(realprecision,round(log(dsc)/(log(10)*pow)*1.2)); dum = round(dsc^(1/pow)); default(realprecision,pres); if(dum^pow == dsc, [1,dum2,dum], print("Discriminant not equal to small primes times a 30'th power"); [0,dsc])) } \\ This takes p^30 out of the discriminant of pol redpdsc(p,pol) = { local(dum); dum = gcd(Mod(1,p)*pol,Mod(1,p)*deriv(pol,x)); while(length(dum) > 2,dum = gcd(Mod(1,p)*dum,Mod(1,p)*deriv(dum,x))); if(length(dum) == 1, pol, subst(pol,x,p*x-lift(component(dum,1)/component(dum,2)))/p^6) } \\ This takes pol and finds a smaller pol such that the two y^2 = pol's are \\ isomorphic. This only works if the discriminant is small. Use wangsmall... \\ mmh is the maximum size of linear transformations that will be tried. \\ lbl is the number of best pol's it keeps and works on. \\ mmh = 5, and lbl = 1 works fine if the discriminant is reduced. smallhyp(pol,mmh,lbl) = { local(dum,k,mh,h); pol = hypred(pol); dum = sizeofpol(pol); blst = vector(lbl,k,[pol,dum*2,1]); mh=2; until(mh == mmh, mh = mh+1; print("Size of linear transforms now : ",mh); k=1; until(k == lbl+1, until((dum <= k) | (h >= mh), h = blst[k][3]; h = h+1; if(h > mh,, blst[k] = [blst[k][1],blst[k][2],h]; dum = smpol(blst[k][1],h))); if(dum <= k, k = dum, k = k+1))); blst[1][1] } smpol(pol,h) = { local(bsz,done,a,b,c,d,dum); done = length(blst)+1; bsz = blst[done-1][2]; for(b=-h,h, for(c = -(h-abs(b)),h-abs(b), for(d = -(h-abs(b)-abs(c)),h-abs(b)-abs(c), a = h-abs(b)-abs(c)-abs(d); if((c == 0 & d==0),, dum = type(simplify((a*x+b)/(c*x+d))); if(dum == "t_POL" || dum == "t_RFRAC", if((dum = sizeofpol(hypsub(pol,a,b,c,d))) < bsz, done = min(done,updblst(slowhypred(hypsub(pol,a,b,c,d)),dum)); bsz = blst[length(blst)][2],)))))); done; } \\ This takes out common factors in the coefficient of pol and also finds \\ the best substitutions of the form x*m or x/m. Best in the sense of \\ minimizing the size of the product of all abs(coefficients). slowhypred(pol) = { local(dum); dum = component(pol,1); for(k=2,7,dum = gcd(dum,component(pol,k))); pol = pol/dum; pol = clrden(subst(pol,x,x/newpowdiv(vector(7,k,component(pol,k))))); dum = newpowdiv(vector(7,k,component(pol,8-k)));\ pol = clrden(subst(pol,x,x*dum)/dum^6) } \\ If vct contains the coefficients of pol (degree 6) in either \\ ascending or descending order this will find the product of all \\ primes, p, to the power k, where k is the biggest integer less than or \\ equal to \\ min{(val(vct[7],p)-1)/3,(val(vct[6],p)-1)/2,(val(vct[5],p)-1),val(vct[4],p)} \\ This will ensure that if we subst x/ans and clear denominators the \\ corresponding pol will have smaller product of coefficients. newpowdiv(vct) = { local(dum,dum2,lst,maxk,kk); dum = vct[4]; for(k=5,7,dum = gcd(dum,vct[k])); dum2 = 1; dum = factor(dum); for(k=1,length(dum~), lst = vector(3,l,myval(vct[l+4],dum[k,1])); maxk = floor(vecmin(vector(3,l,(lst[l]-1)/l))); if(maxk == 0,, lst = concat(vector(3,l,myval(vct[l+1],dum[k,1])),lst); dum2 = dum2*dum[k,1]^(vecsort(vector(maxk,kk,-3*kk-vecmin(vector(6,l,lst[l]-l*kk))),,1)[1]))); dum2 } bigval = 1000; myval(a,p) = if(a==0,bigval,valuation(a,p)) updblst(pol,sz) = { local(len,dum,done); len = length(blst); for(k=1,len, if(sz < blst[len+1-k][2]+10^(-5),dum = len+1-k,)); if(sz > blst[dum][2]-10^(-5),len+1, for(k=1,len-dum, blst[len+1-k] = blst[len-k]); blst[dum] = [pol,sz,1]; if(dum == 1, print("new low size: ",blst[1][2]),); dum) }