/*

  Complex Multiplication Tools for PARI-GP

*/

#define mod(a,b) ((a) < 0 ? ((((a) % (b)) + (b)) % (b)) : ((a) % (b)))
#define gcd(a,b) (abs(gcd00(a,b)))

#include "pari.h"
#include "LL.h"

typedef struct bqf {
  long A;
  long B;
  long C;
} BQF;

GEN elltautoz(GEN,GEN,long);
GEN elltautopoint(GEN,GEN,long);
GEN qfbclassgroup(GEN);
LList qfbclassgroupneg(long);
LList qfbclassgrouppos(long);
BQF qfbclassgroupposR(long,long,BQF);
GEN qfblift(GEN,GEN,GEN);
GEN ellHeegnerPoints(GEN,GEN,GEN,long);
long gcd00(long,long);

GEN elltautoz(GEN e, GEN t, long prec) {
  long ltop = avma;

  long M,j;
  GEN a,q,z,P,digits;

  if (typ(e) != t_VEC) err(talker,"first argument should be elliptic curve");
  if (lg(e) != 20) err(talker,"first argument should be elliptic curve");

  q = gexp(gmul(gmul(gdeux,gmul(mppi(prec),gi)),t),prec);
  digits = ggprecision(q);
  M = 10 * 2 * itos(gceil(gmul( gneg(digits), gdiv( glog(stoi(10),prec), glog(gnorm(q),prec) ) )));
  a = anell(e,M);

  for (j = 1; j <= M; j++) {
    a[j] = (long) gdiv((GEN) a[j],stoi(j));
  }
  P = gtopolyrev(a,0);
  z = gmul(q,gsubst(P,0,q));

  return gerepileupto(ltop,z);
}

GEN elltautopoint(GEN e, GEN t, long prec) {
  long ltop = avma;
  GEN z = elltautoz(e,t,prec);
  GEN point = pointell(e,z,prec);
  return gerepileupto(ltop,point);
}

GEN qfbclassgroup(GEN D) {
  long delta,currpos;
  LList CG;
  LLCell currcell;
  GEN outvec,z;
  BQF F;

  if (typ(D) != t_INT) err(talker,"input must be integer");
  delta = itos(D); /* should also check D is within range of long! */
  if (delta < 0) {
    CG = qfbclassgroupneg(delta);
  } else {
    CG = qfbclassgrouppos(delta);
  }

  outvec = cgetg(LLlength(CG)+1,t_VEC);
  currcell = LLfirst(CG);
  currpos = 1;
  while (currcell != NULL) {
    if (currpos > LLlength(CG)) err(talker,"serious conversion error");
    F = LLGetData(currcell,BQF);
    if (delta > 0) {
      z = cgetg(5,t_QFR);
      z[4] = (long) gzero;
    } else {
      z = cgetg(4,t_QFI);
    }
    z[1] = (long) stoi(F.A);
    z[2] = (long) stoi(F.B);
    z[3] = (long) stoi(F.C);
    outvec[currpos] = (long) z;
    currcell = LLnext(currcell);
    currpos++;
  }
  LLKillList(CG);
  return outvec;
}

LList qfbclassgroupneg(long delta) {
  LList L = LLNewList();
  long A,B,C,Abd;
  BQF F;

  if (delta >= 0) err(talker,"input must be negative integer");
  if ((mod(delta,4) == 2) || (mod(delta,4) == 3)) err(talker,"input must be square mod 4");

  Abd = (long) sqrt(-((float) delta)/3);
  for(A = 1; A <= Abd; A++) {
    for(B = 0; B <= A; B++) {
      if ((B*B - delta) % (4*A) == 0) {
	C = (B*B - delta) / (4*A);
	if (A <= C && gcd(gcd(A,B),gcd(B,C)) == 1) {
	  F.A = A; F.B = B; F.C = C;
	  LLInsertDataAtEnd(L,BQF,F);
	  if (!(B == 0 || B == A || A == C)) {
	    F.A = A; F.B = -B; F.C = C;
	    LLInsertDataAtEnd(L,BQF,F);
	  }
	}
      }
    }
  }

  return L;
}

LList qfbclassgrouppos(long delta) {
  LList L = LLNewList();
  LLCell cellC,cellD;
  long A,B,C,Bbd,Abd1,Abd2,sqrtdelta;
  BQF F,F0,currF;

  sqrtdelta = (long) sqrt((float) delta);
  if ((delta <= 0) || (sqrtdelta*sqrtdelta == delta)) err(talker,"input must be positive nonsquare integer");
  if ((mod(delta,4) == 2) || (mod(delta,4) == 3)) err(talker,"input must be square mod 4");

  /* list the reduced forms */
  Bbd = (long) sqrt((float) delta);
  for(B = 1; B <= Bbd; B++) {
    Abd1 = ((long) ((sqrt((float) delta) - B)/2)) + 1;
    Abd2 = (long) ((sqrt((float) delta) + B)/2);
    for (A = Abd1; A <= Abd2; A++) {
      if ((B*B - delta) % (4*A) == 0) {
	C = (B*B - delta) / (4*A);
        if (gcd(gcd(A,B),gcd(B,C)) == 1) {
      	  F.A = A; F.B = B; F.C = C;
	  LLInsertDataAtEnd(L,BQF,F);
      	  F.A = -A; F.B = B; F.C = -C;
	  LLInsertDataAtEnd(L,BQF,F);
        }
      }
    }
  }


  /* remove the redundant ones by cycling through right neighbors */

  cellC = LLfirst(L);
  while (cellC != NULL) {
    long problem = 0;
    F0 = LLGetData(cellC,BQF);
    /* printf("[%i,%i,%i]",F0.A,F0.B,F0.C); */
    F = qfbclassgroupposR(delta,sqrtdelta,F0);
    while (!((F.A == F0.A) && (F.B == F0.B) && (F.C == F0.C))) {
      /* printf("[%i,%i,%i]",F.A,F.B,F.C); */
      cellD = LLnext(cellC);
      while (cellD != NULL) {
        currF = LLGetData(cellD,BQF);
        if ((currF.A == F.A) && (currF.B == F.B) && (currF.C == F.C)) { LLKillCell(L,cellD); break; }
        cellD = LLnext(cellD);
      }
      F = qfbclassgroupposR(delta,sqrtdelta,F);
    }
    /* printf("\n"); */
    cellC = LLnext(cellC);
  }


  return L;
}

BQF qfbclassgroupposR(long delta, long sqrtdelta, BQF F) {
  BQF RF;
  long B1,Bbd1,Bbd2;
  RF.A = F.C;

  Bbd1 = sqrtdelta - abs(2*RF.A) + 1;
  Bbd2 = sqrtdelta;
  for (B1 = Bbd1; B1 <= Bbd2; B1++) {
    if ((F.B+B1) % abs(2*RF.A) == 0) break;
  }
  RF.B = B1;

  RF.C = (B1*B1 - delta)/(4*RF.A);

  return RF;
}

GEN qfblift(GEN F, GEN GN, GEN Gs) {
  long A,B,C,delta,N,s;
  long A0,B0,C0,a,b,c,d,p;
  long A1,B1,C1,A2,B2,C2;
  GEN t1,t2,L,M;
  long lbot,ltop=avma;

  if (typ(F) != t_QFI) err(talker,"first argument should be imaginary quadratic form");
  A = itos((GEN) F[1]); B = itos((GEN) F[2]); C = itos((GEN) F[3]);
  delta = B*B-4*A*C;
  if ((delta >= 0) || (A <= 0)) err(talker,"first argument should be positive definite quadratic form");
  N = itos(GN); s = itos(Gs);
  if (mod(s*s - delta,N) != 0) err(talker,"third argument must be square root of discriminant (mod N)");

  /* find equivalent form with first coeff prime to N */
  A0 = A; a = 1; c = 0;
  while ((gcd(A0,N) != 1) || (gcd(a,c) != 1)) {
    if (++a == N) { a = 0; if (++c == N) err(talker,"could not find equivalent form with first coeff prime to N"); }
    A0 = A*a*a + B*a*c + C*c*c;
  }
  gbezout(stoi(a),stoi(c),&t1,&t2);
  d = itos(t1); b = -itos(t2);
  B0 = 2*A*a*b + B*(a*d + b*c) + 2*C*c*d;
  C0 = A*b*b + B*b*d + C*d*d;

  /* find equivalent form satisfying congruence condtions */
  p = 0; while (mod(2*A0*p+B0+s,N) != 0) p++;
  A1 = A0*p*p + B0*p + C0; B1 = -(2*A0*p + B0); C1 = A0;

  /* minimize the first coefficient naively */
  M = cgetg(3,t_MAT); t1 = cgetg(3,t_COL); t2 = cgetg(3,t_COL);
  M[1] = (long) t1; M[2] = (long) t2;
  t1[1] = (long) gmul(gdeux,stoi(A1)); t2[1] = (long) gmul(stoi(B1),GN);
  t1[2] = (long) t2[1]; t2[2] = (long) gmul(gmul(gdeux,stoi(C1)),gmul(GN,GN));
  L = qfminim0(M,gzero,gun,0,0);
  a = itos((GEN) (((GEN) ((GEN) L[3])[1])[1])); c = itos((GEN) ((GEN) ((GEN) L[3])[1])[2]);
  if (gcd(a,N*c) != 1) { 
    printf("Warning: problems minimizing in qfblift\n");
    a=1; c=0; 
  }
  gbezout(stoi(a),stoi(N*c),&t1,&t2);
  d = itos(t1); b = -itos(t2);
  A2 = A1*a*a + B1*N*a*c + C1*N*N*c*c;
  B2 = 2*A1*a*b + B1*(a*d + N*b*c) + 2*C1*N*c*d;
  C2 = A1*b*b + B1*b*d + C1*d*d;
  
  /* return result */
  lbot = avma;
  gerepileupto(ltop,(GEN) lbot);
  L = cgetg(4,t_QFI);
  L[1] = (long) stoi(A2); L[2] =  (long) stoi(B2); L[3] = (long) stoi(C2);
  return L;
}

GEN ellHeegnerPoints(GEN E, GEN GD, GEN Gs0, long prec) {
  long D,N,s,h,j;
  GEN CG,GN,pts,F,outvec,Gs=Gs0;
  long ltop = avma;

  if ((typ(E) != t_VEC) || lg(E) != 20) err(talker,"first argument must be elliptic curve");
  if (typ(GD) != t_INT) err(talker,"second argument must be integer");
  D = itos(GD);
  if ((D >= 0) || (mod(D,4) == 2) || (mod(D,4) == 3)) err(talker,"second argument must be quadratic imaginary discriminant");
  N = itos((GEN) globalreduction(E)[1]);
  GN = stoi(N);

  if (Gs == NULL) {
    s = 0;
    while ((mod(D-s*s,N) != 0) && (s < N)) s++;
    if (s == N) err(talker,"second argument is not a square modulo the conductor of the elliptic curve");
    Gs = stoi(s);
  } else {
    if (typ(Gs) != t_INT) err(talker,"third argument must be integer");
  }
  if (mod(D-s*s,N) != 0) err(talker,"third argument must be an orientation");


  CG = qfbclassgroup(GD);
  h = lg(CG) - 1;
  pts = cgetg(lg(CG),t_VEC);
  for (j = 1; j <= h; j++) {
    F = qfblift((GEN) CG[j],GN,Gs);
    pts[j] = (long) gdiv(gsub(gsqrt(GD,prec),(GEN) F[2]),gmul(gdeux,(GEN) F[1]));
  }
  outvec = cgetg(lg(CG),t_VEC);
  for (j = 1; j <= h; j++) {
    outvec[j] = (long) elltautopoint(E,(GEN) pts[j],prec);
  }
  return gerepileupto(ltop,outvec);
}


long gcd00(long a, long b) {
  if (b == 0) return a;
  else return gcd00(b,mod(a,b));
}