/* ===== 1. Auxiliary functions ===== */

quadredtrip(D) =
/* computes the set T_D of the triples (a, b, c) corresponding
 * to the reduced binary forms ax^2 + bxy + cy^2 of discriminant
 * D < 0, b >= 0 */
{
  my (TD = List());
  forstep (b = D % 2, floor(sqrt(-D/3)), 2,
    my (ac = (b^2-D) / 4);
    fordiv (ac, a,
      my(c = ac / a);
      if (a > c, break);
      if (b > a || gcd([a,b,c]) != 1, next);
      listput(TD, [a,b,c]);
      if (b != a && b && a != c, listput(TD, [a,-b,c]));
    );
  );
  listsort(TD); return(TD);
}

elljevaltrip(D, T) = my ([a,b] = T); abs( ellj((-b + sqrt(D)) / (2*a)) );

nfeval(nf, x) = abs(nfeltembed(nf, x, 1)); */

/* ===== 2. Theorem 1.1 ===== */

/* ===== 2.1. The case Q(tau) = Q(tau') ===== */

is0pm1(x) = x == 0 || x == 1 || x == -1;

modcurve(T1, T2, D) =
/* given two triples T1 and T2 of the respective sets T_{4D} and
 * T_D, checks if the corresponding CM-point (j(t1), j(t2))
 * belongs to the modular curve Y_0(2) */
{
  my ([a1,b1] = T1, [a2,b2] = T2, w = quadgen(D));
  my(u1 = -b1 / (2*a1), v1 = 1 / a1, t1 = u1 - v1 + 2*v1*w);
  my(u2 = -b2 / (2*a2), v2 = 1 / (2*a2), t2 = u2 - v2 + 2*v2*w);
  is0pm1(t1 - 2*t2)
     || is0pm1(2*t1 - t2)
     || is0pm1(t1 + 2 / t2)
     || is0pm1(t1 + 2 / (t2+1))
     || is0pm1(t1 + 1 / ((t2+1)/2 - 1));
}

conjcm(D) =
/* returns the full list of the conjugates of (j(t1), j(t2)),
 * where j(t1) and j(t2) are the dominant j-values of respective
 * discriminants 4D and D, keeping only one conjugate for each
 * pair of complex conjugates */
{
  my(TD1 = quadredtrip(4*D), TD2 = quadredtrip(D), C = List());
  for (i = 1, #TD1,
    my(j = 1, T = TD1[i], U);
    if (T[2] < 0, next);
    until(modcurve(T, U, D), U = TD2[j]; j++);
    listput(C, [elljevaltrip(4*D,T), elljevaltrip(D,U), T, U]);
  );
  return(C);
}

boundmn(C, i, j) =
/* given the set C of conjugates of (j(tau), j(tau')), computes
 * an interval for m / n given by c,d such that |m/n - c| < d;
 * with respect to the collinearity of the points C[1], C[i] and
 * C[j], where m and n are positive integers such that 1,
 * j(tau)^m and j(tau')^n are linearly dependent */
{
  my(M, b);
  my([a1,b1] = C[1], [ai,bi] = C[i], [aj,bj] = C[j]);
  if (bi > bj,
    M = (bj/bi + aj/a1 + bj/b1 + aj/ai) / (1 - bj/bi - aj/a1);
    b = bi;
 ,
    M = (aj/ai * (1 + bi/b1) + bi/bj + ai/a1) / (1 - bi/bj - ai/a1);
    b = bj;
  );
  return([log(b1/b)/log(a1/ai), M / ((1-M) * log(a1/ai))]);
}

checkbounds(C) =
/* computes all bounds on m / n for the set of conjugates of
 * (j(tau), j(tau')), where j(tau) and j(tau') are the dominant
 * j-values of respective discriminants 4* D and D, with respect
 * to the previous function boundmn, and checks if two bounds
 * are contradicting each other. In this case, the discriminant
 * D can be eliminated */
{
  my(B = List());
  for(i = 2, #C - 1,
    for(j = i + 1, #C, listput(B, boundmn(C, i, j)));
  );
  for(i = 1, #B - 1,
    my([ci,di] = B[i]);
    for(j = i + 1, #B,
      my([cj,dj] = B[j]);
      if (abs(ci-cj) > di+dj, return(1));
    );
  );
  return(0);
}

approx(alpha, c4,c6,c7) =
{
  my (N = ceil(log(c6)/log((1+sqrt(5))/2)));
  my (pq = contfracpnqn(contfrac(alpha), N));
  for (i = 2, #pq,
    my ([p,q] = pq[,i]);
    if (log(abs(alpha * q - p)) <=  c7 + q * c4, return (0));
    if (q > c6, break)
  );
  return(1);
}

/* C = conjcm(D) */
checklogform(C, D) =
{
  if (D == -23 || D == -31, return (1)); /* two unsolved cases */
  my(h, a, b, M, c1, c3, c4, c5, c6);
  my ([J1,j1, T1,t1] = C[1]);
  my ([J2,j2, T2,t2] = C[2]);
  my ([J3,j3, T3,t3] = C[3]);
  h = qfbclassno(D);
  a = J1 / J2;
  if (j2 > j3,
    b = j1/j2;
    M = (j3/j2 + J3/J1 + j3/j1 + J3/J2) / (1 - j3/j2 - J3/J1)
  ,
    b = j1/j3;
    M = (J3/J2 * (1+j2/j1) + j2/j3 + J2/J1) / (1 - j2/j3 - J2/J1);
  );
  c1 = (-M + (1-M) * log(b)) / ((1-M) * log(a));
  if (j2 > j3,
    c3 = 4 / (1 - j3/j2 - J3/J1);
    c4 = log( max(j3/j2, (J3/J2)^c1) )
  ,
    c3 = 4 / (1 - j2/j3 - J2/J1);
    c4 = log( vecmax([(J3/J2)^c1, j2/j3, (J2/J1)^c1]) );
  );
  c5 = (-2^(6*3+20)*10*19*(2*h)^5 * abs(D) * log(2*exp(1)*h)
       - log(c3/(1-M))) / c4;
  c6 = floor((1 + exp(-1)) * c5 * log(3 * (1 + exp(1)) * c5));

  if (1, \\ block for localbitprec, high accuracy
    localbitprec(floor(log(c6)/log(2)) + 64);
    a = log(elljevaltrip(4*D, T1) / elljevaltrip(4*D, T2));
    b = log(elljevaltrip(D, t1) /
            elljevaltrip(D, if (j2>j3, t2,t3)));
  );
  my (c7 = log(c3 / (a*(1 - M))));
  return (approx(b/a, c4,c6,c7));
}

checkDall() =
/* apply checkbounds(D) on all remaining discriminants |D| <
 * 1024, and prints the list of discriminants that cannot be
 * eliminated this way */
{
  forstep(D = -1023, -1, 8,
    if (qfbclassno(D) < 3, next);
    my (C = conjcm(D));
    if (!checkbounds(C) && !checklogform(C,D), error(D));
  );
}

/* ===== 2.2. The case Q(tau) <> Q(tau') ===== */

conjcmbis(D1, D2) =
/* returns the full list of the conjugates of (j(t1), j(t2)),
 * where j(t1) and j(t2) are the dominant j-values of respective
 * discriminants D1 and D2, D1 and D2 belonging to a common
 * entry of Table 4.3 with [L : Q] >= 4 */

{
  my(C, H1 = polclass(D1,0,'X), H2 = polclass(D2,0,'X));
  my (L = nfinit(polredbest(subst(H1, X, Y))));
  my (k0, x = nfroots(L, H1), x0 = apply(u -> nfeval(L,u), x));
  my (l0, y = nfroots(L, H2), y0 = apply(u -> nfeval(L,u), y));
  vecmax(x0, &k0);
  vecmax(y0, &l0);
  C = apply(s ->
    my(a = nfgaloisapply(L, s, x[k0]));
    my(b = nfgaloisapply(L, s, y[l0]));
    [nfeval(L,a), nfeval(L,b), a, b], nfgaloisconj(L));
  return([L, vecsort(C, (a,b)->sign(b[1]-a[1]))]);
}

checklogformbis(pair) =
{
  my ([D1,D2] = pair);
  my([L,C] = conjcmbis(D1,D2), h, a, b, M, c1, c3, c4, c5, c6, c7, c8 = 0);
  my([J1,j1, T1,t1] = C[1]);
  my([J2,j2, T2,t2] = C[2]);
  my([J3,j3, T3,t3] = C[3]);

  my( t, p0 = 0, p1 = 1, q0 = 1, q1 = 0, p_aux, q_aux, i = 2);
  h = qfbclassno(D1);
  a = J1 / J2;
  if(j2 > j3,
    b = j1/j2;
    M = (j3/j2 + J3/J1 + j3/j1 + J3/J2) / (1 - j3/j2 - J3/J1),
    b = j1 / j3;
    M = (J3/J2 * (1+j2/j1) + j2/j3 + J2/J1) / (1 - j2/j3 - J2/J1);
  );
  c1 = (-M + (1 - M) * log(b)) / ((1 - M) * log(a));
  if(j2 > j3,
    c3 = 4 / (1 - j3/j2 - J3/J1);
    c4 = log( max(j3 / j2, (J3/J2)^c1) )
  ,
    c3 = 4 / (1 - j2/j3 - J2/J1);
    c4 = log( vecmax([(J3/J2)^c1, j2/j3, (J2/J1)^c1]) );
  );
  c5 = (- 2748779069440 * h^5 * abs(D1*D2)^(1/2) * log(exp(1) * h)
       - log(c3 / (1 - M))) / c4;
  c6 = floor((1 + exp(-1)) * c5 * log(3 * (1 + exp(1)) * c5));

  if (1, \\ block for localbitprec, high accuracy
    localbitprec(floor(log(c6)/log(2)) + 64);
    a = log(nfeval(L,T1) / nfeval(L,T2));
    b = log(nfeval(L,t1) / nfeval(L, if(j2>j3, t2,t3)));
  );
  my (c7 = log(c3 / (a*(1-M))));
  return (approx(b/a, c4,c6,c7));
}

checklogformall() =
/* checks if checklogformbis(D1, D2) = 1, for each pair of
 * discriminants {D1, D2} given by a common entry of Table 4.3
 * with [L : Q] >= 4 */
{
  localprec(50);
  my(LD = [[-96,-192], [-96,-288], [-192,-288], [-180,-240], [-195,-520],
  [-195,-715], [-520,-715], [-120,-160], [-120,-280], [-120,-760], [-160,-280],
  [-160,-760], [-280,-760], [-340,-595], [-480,-960]]);
  for(k = 1, #LD,
    if (!checklogformbis(LD[k]), error(LD[k]));
  );
}

/* ===== 3. Theorem 1.2 ===== */

multind(L, x, y) =
/* checks if the algebraic numbers x and y of the number field L
 * are multiplicatively independent or not. x and y are two
 * algebraic numbers of modulus <> 1 */
{
  my (p, fa, f, k, l, u0, u, v);
  fa = idealfactor(L, idealadd(L,x,y));
  if (!fa,
    if (idealhnf(L,x) != 1 || idealhnf(L,y) != 1, return (1));
    error("unit case, not handled");
  );
  p = fa[1,1];
  f = idealval(L, y, p) / idealval(L, x, p);
  k = numerator(f);
  l = denominator(f);
  x = nfbasistoalg(L,x);
  y = nfbasistoalg(L,y); v = x^k / y^l;
  u0 = nfbasistoalg(L, nfrootsof1(L)[2]);
  u = u0;
  while(u <> 1 && u <> v, u *= u0);
  return(u == 1);
}

elimmultind(pair) =
/* checks if j(tau)^m * j(tau')^n is not rational, for any
 * singular moduli j(tau) and j(tau') of respective
 * discriminants D1 and D2, and for any integers m, n <> 0 */
{
  my ([D1,D2] = pair);
  my (H1 = polclass(D1,0,'y), H2 = polclass(D2,0,'y), L, x, y);
  print(pair);
  if (quaddisc(D1)==quaddisc(D2),
    L = subst(rnfequation(polredbest(H1), 'x^2-D1), 'x,'y)
  , L = H1);
  L = nfinit(polredbest(L));
  x = nfroots(L, subst(H1,'y,'x));
  y = nfroots(L, subst(H2,'y,'x));
  if (#y < #x, return(1));
  /* Gx[i] = sigma_i . x */
  my (G = nfgaloisconj(L));
  Gx = [ apply(t->nfgaloisapply(L,s,t),x) | s <- G ];
  for (i = if (D1<>D2, 1, 2), #x,
    my(k = 1);
    while (k <= #G,
      my (sx = Gx[k]);
      if (sx[1] != x[1] &&
          (D1<>D2 || [sx[i],sx[1]] <> [x[1],x[i]]), break);
      k++;
    );
    if (k > #G || !multind(L, x[1] / Gx[k][1],
                              nfgaloisapply(L, G[k], y[i]) / y[i]),
      return(0);
    );
  );
  return(1);
}

elimmultindall() =
/* checks if elimmultind(D1, D2) = 1, for each remaining pair
 * {D1, D2}; see Subsections 5.2 and 5.3. */
{
  my(LD2 = List(), LD);
  LD = List([[-71,-71], [-95,-95], [-96,-96],
  [-96,-192],  [-96,-288],  [-192,-192], [-192,-288], [-288,-288],
  [-180,-180], [-180,-240], [-240,-240], [-195,-195], [-195,-520],
  [-195,-715], [-520,-520], [-520,-715], [-715,-715], [-120,-120],
  [-120,-160], [-120,-280], [-120,-760], [-160,-160], [-160,-280],
  [-160,-760], [-280,-280], [-280,-760], [-760,-760], [-340,-340],
  [-340,-595], [-595,-595], [-480,-480], [-480,-960], [-960,-960]
  ]);
  for (D = -427, -1,
    if (D%4 <= 1 && qfbclassno(D) == 2, listput(LD2, D));
  );
  for (k = 1, #LD2-1,
    for (l = k+1, #LD2, listput(LD, [LD2[k], LD2[l]]));
  );
  forstep (D = -255, -1, 8,
    if (qfbclassno(D) >= 3, listput(LD, [D, 4*D]));
  );
  for (D = -4075, -1,
    if (D % 4 > 1, next);
    my (h = qfbclassno(D));
    if (h >= 3 && h <= 6,
      listput(LD, [D, D]);
      if (D % 8 == 1, listput(LD, [D, 4*D]));
    );
  );
  forstep (D = -303, -1, 4,
    if (qfbclassno(D) > 6, listput(LD, [D, 4*D]));
  );
  for (k = 1, #LD,
    if (!elimmultind(LD[k]), error(LD[k]));
  );
}