\\ 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])<abs(b[1]),
      if(type(dum=bquadratq(a,b)) == "t_INT",,dum = dum[2]),
      if(type(dum=bquadratq(b,a)) == "t_INT",,dum = [dum[2][1],dum[2][3],dum[2][2]]));
    if(type(dum) == "t_INT",dum,
      [mesmat,[dum[1],dum[2]*am,dum[3]*bm]*dmat]))
}

\\ This multiplies two rational numbers given in factored form
fmult(a,b) =
{
  local(dum,dum2,outv,m);

  dum = concat(a~,b~);
  dum2 = vecsort(vector(length(dum),k,dum[1,k]),,1);
  outv = [dum[1,dum2[1]];dum[2,dum2[1]]];
  m = 1;
  for(k=2,length(dum2),
    if(outv[1,m] == dum[1,dum2[k]],
      outv[2,m] = outv[2,m]+dum[2,dum2[k]],
      m = m+1;
      outv = concat(outv,[dum[1,dum2[k]];dum[2,dum2[k]]])));
  outv~
}

\\ This takes the reciprocal of a rational number given in factored form
frec(a) = a*[1,0;0,-1]

\\ This takes the reciprocal of a rational number given in factored form
frec(a) = a*[1,0;0,-1]

\\ This finds the denominator of a rational number given in factored form
fdenom(a) =\
  for(k=1,length(a~),a[k,2] = if(a[k,2] < 0,-a[k,2],0));\
  a

\\ This finds the square root of the biggest square in a rational number
\\ given in factored form
fsqtpart(a) =
{
  for(k=1,length(a~),a[k,2] = (a[k,2] \ 2));
  a
}

\\ This find the rational number given a rational number in factored form
f2r(a) = prod(k=1,length(a~),a[k,1]^a[k,2])

\\ This finds the denominator of a rational number given in factored form
fdenom(a) =
{
  for(k=1,length(a~),a[k,2] = if(a[k,2] < 0,-a[k,2],0));
  a
}

\\The number of squareroots of a mod b to try
numsqr = 20;
fnum = 0;
factlist = vector(40,k,0);

\\ This finds a zero over Q of the quadratic form x^2-a*y^2-b*z^2 if it exists.
\\ It is assumed that abs(a) =< abs(b) and that a and b are squarefree.
\\ a and b should be vectors with the first entry the actual number and the 
\\ second the number in factored form. 
\\ This even works for really big a and b.
bquadratq(a,b) =
{
  local(dum,dum2,lst,nb,k);

  if(a[1] == 1,[1,[1,1,0]],
  if(abs(b[1])==1,if(a[1]==-1 && b[1]==-1,0,
    if(a[1] == 1,[1,[1,1,0]],[1,[1,0,1]])),
    if(type(dum = mysqrt(a[1],b[2])) != "t_INTMOD",0,
      dum = component(dum,2);
      if(abs(b[1]) < 10^35,
        if(abs(dum-b[1]) < abs(dum),dum = dum-b[1],);
        nb = [dum2 = ((dum^2-a[1])/b[1]),factor(dum2)];
        b = fsqrfrprt(nb);
        nb = f2r(fsqtpart(nb[2])),
        dum = vector(numsqr,k,dum+(k-floor(numsqr/2))*b[1]);
        lst = [];
        for(k=1,numsqr,
          dum2 = smallsqfrpt((dum[k]^2-a[1])/b[1]);
          if(abs(dum2[1]) < abs(b[1]),
            if(isprime(dum2[2]),dum2[2]=dum2[2]^(0.8),);
            lst = concat(lst,[[(dum[k]^2-a[1])/b[1],dum2[2],dum[k]]]),));
        dum = vecsort(vector(length(lst),k,lst[k][2]),,1);
        dum2 = lst[dum[1]][1];
        if((abs(lst[dum[1]][2]) > 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])<abs(b[1]),
        dum2 = bquadratq(a,b),
        if(type(dum2 = bquadratq(b,a))== "t_INT",,
          dum2 = [1,[dum2[2][1],dum2[2][3],dum2[2][2]]]));
      if(type(dum2) == "t_INT",dum2,
        if(dum2[2][3] == 0,dum2,
          dum2=[dum2[2][1]*dum+a[1]*dum2[2][2],dum2[2][2]*dum+dum2[2][1],b[1]*dum2[2][3]*nb];
          dum = gcd(gcd(dum2[1],dum2[2]),dum2[3]);
         [1,dum2/dum])))))
}

\\ This finds a squareroot of a mod m if m is squarefree. m should be
\\ given in factored form
mysqrt(a,m) =
{
  local(srt,nsqr);

  srt = Mod(0,1);
  nsqr = 0;
  for(k=1,length(m~),
    if((m[k,1] == -1) | (m[k,2] == 0),,
    if(kronecker(a,m[k,1]) == -1 && m[k,1] != 2,
      nsqr = 1,
      srt = chinese(srt,sqrt(Mod(a,m[k,1]))))));
  if(nsqr == 0,srt,0)
}

\\ This finds the squarefree part of a factored number, where the number is
\\ given as a two component vector, the first the number the second the
\\ factorization
fsqrfrprt(a) =
{
  a = a[2];
  for(k=1,length(a~),a[k,2] = (a[k,2] % 2));
  [f2r(a),a]
}

\\ This takes out only squares of small factors
smallsqfrpt(int) =
{
  local(dum);

  if(int == 1,[1,1],
    dum = factor(int,0);
    [prod(k=1,length(dum~),dum[k,1]^(dum[k,2] % 2)),dum[length(dum~),1]])
}

\\ This uses Wang's idea to reduce the discrimininant. 
\\ Then it smallhyp's the remaining bigpol.
wangsmall(pol) =
{
  local(dum,p,m);

  pol = nsmpol(pol);
  dum = facdisc(poldisc(pol));
  if(dum[1] == 0,pol,
  m = dum[3];
  if(m == 1,,
    pol = redpdsc(m,pol));
  dum = dum[2];
  for(k=1,length(dum~),
    if(dum[k,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)
}