/*

  The tate collection

*/

#include "pari.h"

GEN cpnonsquare(GEN);
GEN cpadd(GEN,GEN);
GEN cpsub(GEN,GEN);
GEN cpmul(GEN,GEN);
GEN cpdiv(GEN,GEN);
GEN cpsqr(GEN);
GEN cpinv(GEN);
GEN cpneg(GEN);
GEN cpsqrt(GEN,GEN,long);
GEN cpchangew(GEN,GEN);
GEN cptoquad(GEN);
GEN cptomat(GEN);
GEN cppow(GEN,long);
GEN ellchangeofvars(GEN,GEN,GEN,long);
GEN elltate(GEN,GEN);

GEN cpnonsquare(GEN Gp) {
  /* this should be modified sometime to handle arbitrary sized primes, and to use gcarreparfait */
  long ltop = avma;
  long i,j,*L,p = itos(Gp);
  if (p == 2) err(talker,"cpgen does work at p=2");
  else if ((p % 4) == 3) return stoi(-1);
  else if ((p % 8) == 5) return gdeux;
  else {
    L = malloc(p*sizeof(long));
    for (i = 0; i < p; i++) L[i] = 0;
    for (i = 0; i < p; i++) L[(i*i) % p] = 1;
    for (i = 1; i < p; i++) {
      if (L[i] = 0) { j = i; break; }
      else if (L[p-i] == 0) { j = -i; break; }
    }
    free(L);
    return(stoi(j));
  }
}

GEN cpadd(GEN a, GEN b) {
  GEN outvec;

  if (typ(a) == t_VEC && lg(a) == 4 && typ(b) != t_VEC) {
    outvec = cgetg(4,t_VEC);
    outvec[1] = ladd((GEN) a[1],b);
    outvec[2] = lcopy((GEN) a[2]);
    outvec[3] = lcopy((GEN) a[3]);
  } else if (typ(a) != t_VEC && typ(b) == t_VEC && lg(b) == 4) {
    outvec = cgetg(4,t_VEC);
    outvec[1] = ladd((GEN) b[1],a);
    outvec[2] = lcopy((GEN) b[2]);
    outvec[3] = lcopy((GEN) b[3]);
  } else if (typ(a) == t_VEC && lg(a) == 4 && typ(b) == t_VEC && lg(b) == 4) {
    if (gegal((GEN) a[3],(GEN) b[3])) {
      outvec = cgetg(4,t_VEC);
      outvec[1] = ladd((GEN) a[1],(GEN) b[1]);
      outvec[2] = ladd((GEN) a[2],(GEN) b[2]);
      outvec[3] = lcopy((GEN) a[3]);
    } else err(talker,"nonmatching choices of nonsquare");
  } else {
    outvec = gadd(a,b);
  }

  return outvec;
}

GEN cpsub(GEN a, GEN b) {
  GEN outvec;

  if (typ(a) == t_VEC && lg(a) == 4 && typ(b) != t_VEC) {
    outvec = cgetg(4,t_VEC);
    outvec[1] = lsub((GEN) a[1],b);
    outvec[2] = lcopy((GEN) a[2]);
    outvec[3] = lcopy((GEN) a[3]);
  } else if (typ(a) != t_VEC && typ(b) == t_VEC && lg(b) == 4) {
    outvec = cgetg(4,t_VEC);
    outvec[1] = lsub(a,(GEN) b[1]);
    outvec[2] = lneg((GEN) b[2]);
    outvec[3] = lcopy((GEN) b[3]);
  } else if (typ(a) == t_VEC && lg(a) == 4 && typ(b) == t_VEC && lg(b) == 4) {
    if (gegal((GEN) a[3],(GEN) b[3])) {
      outvec = cgetg(4,t_VEC);
      outvec[1] = lsub((GEN) a[1],(GEN) b[1]);
      outvec[2] = lsub((GEN) a[2],(GEN) b[2]);
      outvec[3] = lcopy((GEN) a[3]);
    } else err(talker,"nonmatching choices of nonsquare");
  } else {
    outvec = gsub(a,b);
  }

  return outvec;
}

GEN cpmul(GEN a, GEN b){
  GEN outvec;
  GEN a1b1,a2b2s,a1b2,a2b1;

  if (typ(a) == t_VEC && lg(a) == 4 && typ(b) != t_VEC) {
    outvec = cgetg(4,t_VEC);
    outvec[1] = lmul((GEN) a[1],b);
    outvec[2] = lmul((GEN) a[2],b);
    outvec[3] = lcopy((GEN) a[3]);
  } else if (typ(a) != t_VEC && typ(b) == t_VEC && lg(b) == 4) {
    outvec = cgetg(4,t_VEC);
    outvec[1] = lmul(a,(GEN) b[1]);
    outvec[2] = lmul(a,(GEN) b[2]);
    outvec[3] = lcopy((GEN) b[3]);
  } else if (typ(a) == t_VEC && lg(a) == 4 && typ(b) == t_VEC && lg(b) == 4) {
    if (gegal((GEN) a[3],(GEN) b[3])) {
      a1b1 = gmul((GEN) a[1], (GEN) b[1]);
      a2b2s = gmul(gmul((GEN) a[2], (GEN) b[2]),(GEN) a[3]);
      a1b2 = gmul((GEN) a[1], (GEN) b[2]);
      a2b1 = gmul((GEN) a[2], (GEN) b[1]);
      outvec = cgetg(4,t_VEC);
      outvec[1] = ladd(a1b1,a2b2s);
      outvec[2] = ladd(a1b2,a2b1);
      outvec[3] = lcopy((GEN) a[3]);
    } else err(talker,"nonmatching choices of nonsquare");
  } else {
    outvec = gmul(a,b);
  }

  return outvec;
}

GEN cpdiv(GEN a, GEN b) {
  GEN outvec;
  GEN denom,ab1,mab2;
  GEN a1b1msa2b2,a2b1ma1b2;

  if (typ(a) == t_VEC && lg(a) == 4 && typ(b) != t_VEC) {
    outvec = cgetg(4,t_VEC);
    outvec[1] = ldiv((GEN) a[1],b);
    outvec[2] = ldiv((GEN) a[2],b);
    outvec[3] = lcopy((GEN) a[3]);
  } else if (typ(a) != t_VEC && typ(b) == t_VEC && lg(b) == 4) {
    denom = gsub(gsqr((GEN) b[1]),gmul((GEN) b[3],gsqr((GEN) b[2])));
    ab1 = gmul(a,(GEN) b[1]);
    mab2 = gneg(gmul(a,(GEN) b[2]));
    outvec = cgetg(4,t_VEC);
    outvec[1] = ldiv(ab1,denom);
    outvec[2] = ldiv(mab2,denom);
    outvec[3] = lcopy((GEN) b[3]);
  } else if (typ(a) == t_VEC && lg(a) == 4 && typ(b) == t_VEC && lg(b) == 4) {
    if (gegal((GEN) a[3],(GEN) b[3])) {
      denom = gsub(gsqr((GEN) b[1]),gmul((GEN) b[3],gsqr((GEN) b[2])));
      a1b1msa2b2 = gsub(gmul((GEN) a[1],(GEN) b[1]),gmul((GEN) b[3],gmul((GEN) a[2], (GEN) b[2])));
      a2b1ma1b2 = gsub(gmul((GEN) a[2],(GEN) b[1]),gmul((GEN) a[1],(GEN) b[2]));
      outvec = cgetg(4,t_VEC);
      outvec[1] = ldiv(a1b1msa2b2,denom);
      outvec[2] = ldiv(a2b1ma1b2,denom);
      outvec[3] = lcopy((GEN) a[3]);
    } else err(talker,"nonmatching choices of nonsquare");
  } else {
    outvec = gdiv(a,b);
  }

  return outvec;
}

GEN cpsqr(GEN a) {
  GEN outvec;
  GEN a1sqr,sa2sqr,a1a2;

  if (typ(a) == t_VEC && lg(a) == 4) {
    a1sqr = gsqr((GEN) a[1]);
    sa2sqr = gmul((GEN) a[3],gsqr((GEN) a[2]));
    a1a2 = gmul((GEN) a[1],(GEN) a[2]);
    outvec = cgetg(4,t_VEC);
    outvec[1] = ladd(a1sqr,sa2sqr);
    outvec[2] = lmul(gdeux,a1a2);
    outvec[3] = lcopy((GEN) a[3]);
  } else {
    outvec = gsqr(a);
  }

  return outvec;
}

GEN cpinv(GEN a) {
  GEN outvec;
  GEN denom,nega2;

  if (typ(a) == t_VEC && lg(a) == 4) {
    denom = gsub(gsqr((GEN) a[1]),gmul((GEN) a[3],gsqr((GEN) a[2])));
    nega2 = gneg((GEN) a[2]);
    outvec = cgetg(4,t_VEC);
    outvec[1] = ldiv((GEN) a[1],denom);
    outvec[2] = ldiv(nega2,denom);
    outvec[3] = lcopy((GEN) a[3]);
  } else {
    outvec = ginv(a);
  }

  return outvec;
}

GEN cpneg(GEN a) {
  GEN outvec;

  if (typ(a) == t_VEC && lg(a) == 4) {
    outvec = cgetg(4,t_VEC);
    outvec[1] = lneg((GEN) a[1]);
    outvec[2] = lneg((GEN) a[2]);
    outvec[3] = lcopy((GEN) a[3]);
  } else {
    outvec = gneg(a);
  }

  return outvec;

}

GEN cpsqrt(GEN x, GEN d1, long prec) {
  GEN outvec,temp,d;

  if (typ(x) == t_PADIC) {
    if (kronecker((GEN) x[4],(GEN) x[2]) == -1) {
      if (d1 == NULL) d = cpnonsquare((GEN) x[2]);
      else d = d1;
      temp = gdiv(x,d);
      outvec = cgetg(4,t_VEC);
      outvec[1] = zero;
      outvec[2] = lsqrt(temp,0);
      outvec[3] = lcopy(d);
    } else {
      outvec = gsqrt(x,0);
    }
  } else {
    outvec = gsqrt(x,prec);
  }

  return outvec;
}

GEN cpchangew(GEN a, GEN newd) {
  GEN z,a1,a2,a3,z2,p,d,one1,temp;

  if (typ(a) == t_VEC && lg(a) == 4) {

    a1 = (GEN) a[1]; a2 = (GEN) a[2]; a3 = (GEN) a[3];
    z2 = NULL;

    if ((typ(a2) == t_PADIC) ? gcmp0((GEN) a2[4]) : gcmp0(a2)) {
      z2 = a2;
    } else if (typ(a2) == t_PADIC) {
      p = (GEN) a2[2];
      one1 = gadd(gun,gsub(a2,a2));
    } else if (typ(a1) == t_PADIC) {
      p = (GEN) a1[2];
      one1 = gadd(gun,gsub(a1,a1));
    } else err(talker,"need specification of p in arguments to cpchangew");

    if (z2 == NULL) {
      if (newd == NULL) d = cpnonsquare(p);
      else if (typ(newd) == t_INT) d = newd;
      else err(talker,"second argument to cpchangew should be integer");
      
      temp = gsqrt(gmul(gdiv((GEN) a[3],d),one1),0);
      z2 = gmul(a2,temp);
    }

    if (newd == NULL) d = cpnonsquare(p);
    else if (typ(newd) == t_INT) d = newd;
    else err(talker,"second argument to cpchangew should be integer");

    temp = gsqrt(gmul(gdiv((GEN) a[3],d),one1),0);
    z = cgetg(4,t_VEC);
    z[1] = lcopy(a1);
    z[2] = lcopy(z2);
    z[3] = lcopy(d);

  } else {
    z = gcopy(a);
  }

  return z;
}

GEN cptoquad(GEN a) {
  GEN z,mp;

  if (typ(a) == t_VEC && lg(a) == 4) {
    z = cgetg(4,t_QUAD);
    mp = cgetg(5,t_POL);
    setsigne(mp,1); setlgef(mp,5); setvarn(mp,0);
    mp[2] = lneg((GEN) a[3]); mp[3] = zero; mp[4] = un;
    z[1] = (long) mp; z[2] = lcopy((GEN) a[1]); z[3] = lcopy((GEN) a[2]);
  } else {
    z = gcopy(a);
  }

  return z;
}

GEN cptomat(GEN a) {
  GEN z,c1,c2;

  if (typ(a) == t_VEC && lg(a) == 4) {
    z = cgetg(3,t_MAT);
    c1 = cgetg(3,t_COL); c2 = cgetg(3,t_COL);
    c1[1] = lcopy((GEN) a[1]); c2[1] = lmul((GEN) a[2], (GEN) a[3]);
    c1[2] = lcopy((GEN) a[2]); c2[2] = lcopy((GEN) a[1]);
    z[1] = (long) c1; z[2] = (long) c2;
  } else {
    z = gcopy(a);
  }

  return z;
}

GEN cppow(GEN ina, long inmu) {
  long t,mu;
  GEN a,PROD,ato2toi;

  if (inmu < 0) {
    a = cpinv(ina);
    mu = -inmu;
  } else {
    a = ina;
    mu = inmu;
  }

  ato2toi = a;
  PROD = gun;

  while (mu != 0) {
    if (t = mu % 2 == 1) PROD = cpmul(PROD,ato2toi);
    mu = (mu-t) / 2;
    ato2toi = cpmul(ato2toi,ato2toi);
  }

  return PROD;
}

GEN ellchangeofvars(GEN E1, GEN E2, GEN d, long prec) {
  long ltop = avma,lbot;
  GEN u,usqr,r,s,t,outvec;

  usqr = gdiv(gmul((GEN) E1[11],(GEN) E2[10]),gmul((GEN) E1[10],(GEN) E2[11]));
  u = cpsqrt(usqr,d,prec);
  s = cpdiv(cpsub(cpmul(u,(GEN) E2[1]),(GEN) E1[1]),gdeux);
  r = gsub(gdiv(gsub(gmul(usqr,(GEN) E2[2]),(GEN) E1[2]),stoi(3)),gdiv(gsub(gsqr((GEN) E1[1]),gmul(usqr,gsqr((GEN) E2[1]))),stoi(12)));
  t = cpadd(cpsub(cpdiv(cpsub(cpmul(cpmul(u,usqr),(GEN) E2[3]),(GEN) E1[3]),gdeux),cpdiv(cpsub(cpmul(cpmul(usqr,(GEN) E1[1]),(GEN) E2[2]),cpmul((GEN) E1[1],(GEN) E1[2])),stoi(6))),cpdiv(cpsub(cpmul((GEN) E1[1],cpmul((GEN) E1[1],(GEN) E1[1])),gmul(usqr,gmul((GEN) E1[1],gsqr((GEN) E2[1])))),stoi(24)));

  outvec = cgetg(5,t_VEC);
  outvec[1] = lcopy(u); outvec[2] = lcopy(r); outvec[3] = lcopy(s); outvec[4] = lcopy(t);
  return gerepileupto(ltop,outvec);
}

GEN elltate(GEN E, GEN c) {
  long ltop = avma,ltop2,lbot2;
  long i,N;
  GEN q,Eq,X,Y,v,d,qtoi;
  GEN Gi,t1,t3,t5,sa,s1,s3,s5,P1,P3,P5;
  GEN T1,T2,T3,T4,sX,sY;
  GEN u,r,s,t,usqrinv,usqrsqrinv,pt;

  if (typ(E) != t_VEC || lg(E) != 20) err(talker,"first argument should be elliptic curve");
  q = (GEN) E[17];
  if (typ(q) != t_PADIC) err(talker,"first argument should be p-adic elliptic curve");
  N = precp(q);

  /* compute s1,s3,s5 */
  /* this may be a waste if gsubst doesn't use Horner's rule */
  ltop2 = avma;
  t1 = cgetg((N+3),t_VEC); t3 = cgetg((N+3),t_VEC); t5 = cgetg((N+3),t_VEC);
  t1[1] = zero; t3[1] = zero; t5[1] = zero;
  for (i = 2; i <= N+2; i++) {
    Gi = stoi(i-1);
    t1[i] = (long) sumdiv(Gi);
    t3[i] = (long) gsumdivk(Gi,3);
    t5[i] = (long) gsumdivk(Gi,5);
  }
  P1 = gtopolyrev(t1,0); P3 = gtopolyrev(t3,0); P5 = gtopolyrev(t5,0);
  sa = cgetg(4,t_VEC);
  sa[1] = lsubst(P1,0,q); sa[2] = lsubst(P3,0,q); sa[3] = lsubst(P5,0,q);
  sa = gerepileupto(ltop2,sa);
  s1 = (GEN) sa[1]; s3 = (GEN) sa[2]; s5 = (GEN) sa[3];

  /* determine the curve E_q */
  Eq = cgetg(6,t_VEC);
  Eq[1] = un; Eq[2] = zero; Eq[3] = zero;
  Eq[4] = lmul(stoi(-5),s3);
  Eq[5] = ldiv(gadd(gmul(stoi(5),s3),gmul(stoi(7),s5)),stoi(-12));
  Eq = initell(Eq,0);

  /* compute the point on E_q */
  ltop2 = avma;
  sa = cgetg(3,t_VEC);
  qtoi = gun;
  T1 = gun;
  T2 = c;
  T3 = cpsub(gun,T2);
  T4 = cpsqr(T3);
  sX = cpdiv(T2,T4);
  sY = cpdiv(cpsqr(T2),cpmul(T3,T4));
  for (i = 1; i <= N+2; i++) {

    qtoi = gmul(qtoi,q);
    T1 = qtoi;
    T2 = cpmul(T1,c);
    T3 = cpsub(gun,T2);
    T4 = cpsqr(T3);
    sX = cpadd(sX,cpdiv(T2,T4));
    sY = cpadd(sY,cpdiv(cpsqr(T2),cpmul(T3,T4)));

    T1 = ginv(T1);
    T2 = cpmul(T1,c);
    T3 = cpsub(gun,T2);
    T4 = cpsqr(T3);
    sX = cpadd(sX,cpdiv(T2,T4));
    sY = cpadd(sY,cpdiv(cpsqr(T2),cpmul(T3,T4)));
  }
  sa = cgetg(3,t_VEC);
  sa[1] = lcopy(sX); sa[2] = lcopy(sY);
  sa = gerepileupto(ltop2,sa);
  X = cpsub((GEN) sa[1],cpmul(gdeux,s1));
  Y = cpadd((GEN) sa[2],s1);

  /* change variables and return */
  if (typ(c) == t_VEC && lg(c) == 4) {
    d = (GEN) c[3];
  } else {
    d = NULL;
  }
  v = ellchangeofvars(Eq,E,d,0);
  u = (GEN) v[1]; r = (GEN) v[2]; s = (GEN) v[3]; t = (GEN) v[4];
  usqrinv = cpinv(cpsqr(u)); usqrsqrinv = cpsqr(usqrinv);
  T1 = cpmul(cpsub(X,r),usqrinv);
  T2 = cpmul(cpmul(u,cpadd(cpsub(Y,cpmul(s,X)),cpsub(cpmul(s,r),t))),usqrsqrinv);
  pt = cgetg(3,t_VEC);
  pt[1] = lcopy(T1); pt[2] = lcopy(T2);
  return gerepileupto(ltop,pt);
}