/*
  Author:  Pate Williams (c) 1997

  Heuristic Algorithm 5.4.10 (h(D) Using Baby-Step
  Giant-Step). See "A Course in Computational
  Algebraic Number Theory" by Henri Cohen page 251.
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "lip.h"

struct zform { verylong za, zb, zc; };

int compare_zform(const void *f1, const void *f2)

{
  struct zform *zform1 = (struct zform *) f1;
  struct zform *zform2 = (struct zform *) f2;
  int i = zcompare(zform1->za, zform2->za), j;

  if (i != 0) return i;
  j = zcompare(zform1->zb, zform2->zb);
  if (j != 0) return j;
  return zcompare(zform1->zc, zform2->zc);
}

void copy_zform(struct zform zform1, struct zform *zform2)
{
  zcopy(zform1.za, &zform2->za);
  zcopy(zform1.zb, &zform2->zb);
  zcopy(zform1.zc, &zform2->zc);
}

void free_zform(struct zform *zform1)
{
  zfree(&zform1->za);
  zfree(&zform1->zb);
  zfree(&zform1->zc);
}

void get_zform(verylong *zD, struct zform *zform1)
{
  verylong zd = 0, ze = 0, zf = 0;

  printf("a = "); zread(&zform1->za);
  printf("b = "); zread(&zform1->zb);
  printf("c = "); zread(&zform1->zc);
  zmul(zform1->za, zform1->zc, &zd);
  zlshift(zd, 2l, &ze);
  zsq(zform1->zb, &zf);
  zsub(zf, ze, zD);
  printf("D = "); zwriteln(*zD);
  zfree(&zd);
  zfree(&ze);
  zfree(&zf);
}

void init_zform(struct zform *zform1)
{
  zform1->za = zform1->zb = zform1->zc = 0;
}

void put_zform(struct zform zform1)
{
  verylong zD = 0, zd = 0, ze = 0, zf = 0;

  zmul(zform1.za, zform1.zc, &zd);
  zlshift(zd, 2l, &ze);
  zsq(zform1.zb, &zf);
  zsub(zf, ze, &zD);
  printf("a = "); zwriteln(zform1.za);
  printf("b = "); zwriteln(zform1.zb);
  printf("c = "); zwriteln(zform1.zc);
  printf("D = "); zwriteln(zD);
  zfree(&zD);
  zfree(&zd);
  zfree(&ze);
  zfree(&zf);
}

struct zform *search_zform(long n, struct zform key,
                           struct zform *array)
{
  long i;
  struct zform *address;

  for (i = 0; i < n; i++) {
    address = &array[i];
    if (compare_zform(&key, address) == 0) return address;
  }
  return 0;
}

void reduce(verylong *za, verylong *zb, verylong *zc)
{
  verylong zd = 0, ze = 0, zf = 0, zq = 0, zr = 0;

  zcopy(*za, &zd);
  znegate(&zd);
  if (zcompare(zd, *zb) < 0 && zcompare(*zb, *za) <= 0)
    goto L3;
  L2:
  zlshift(*za, 1l, &zd);
  zdiv(*zb, zd, &zq, &zr);
  if (zscompare(zr, 0l) < 0) {
    if (zscompare(zd, 0l) < 0) znegate(&zd);
    zadd(zr, zd, &ze);
    znegate(&zd);
    zcopy(ze, &zr);
  }
  if (zcompare(zr, *za) > 0) {
    zsub(zr, zd, &ze);
    zcopy(ze, &zr);
    zsadd(zq, 1l, &ze);
    zcopy(ze, &zq);
  }
  zadd(*zb, zr, &zd);
  zmul(zd, zq, &ze);
  zrshift(ze, 1l, &zf);
  zsub(*zc, zf, &zd);
  zcopy(zd, zc);
  zcopy(zr, zb);
  L3:
  if (zcompare(*za, *zc) > 0) {
    znegate(zb);
    zcopy(*za, &zd);
    zcopy(*zc, za);
    zcopy(zd, zc);
    goto L2;
  }
  if (zcompare(*za, *zc) == 0 && zscompare(*zb, 0l) < 0)
    znegate(zb);
  zfree(&zd);
  zfree(&ze);
  zfree(&zf);
  zfree(&zq);
  zfree(&zr);
}

void PARTEUCL(verylong zL, verylong za,
              verylong zb, verylong *zv2,
              verylong *zv3, verylong *zd,
              verylong *zz)
{
  verylong zq = 0, zr = 0, zs = 0, zt = 0, zv = 0;
  verylong zt2 = 0, zt3 = 0;

  zzero(&zv);
  zcopy(za, zd);
  zone(zv2);
  zcopy(zb, zv3);
  zzero(zz);
  L2:
  zcopy(*zv3, &zr);
  zabs(&zr);
  if (zcompare(zr, zL) > 0) goto L3;
  if (zodd(*zz)) {
    znegate(zv2);
    znegate(zv3);
  }
  zfree(&zq);
  zfree(&zr);
  zfree(&zs);
  zfree(&zt);
  zfree(&zv);
  zfree(&zt2);
  zfree(&zt3);
  return;
  L3:
  zdiv(*zd, *zv3, &zq, &zt3);
  if (zscompare(zt3, 0l) < 0) {
    zadd(zt3, zr, &zt);
    zcopy(zt, &zt3);
  }
  zmul(zq, *zv2, &zs);
  zsub(zv, zs, &zt2);
  zcopy(*zv2, &zv);
  zcopy(*zv3, zd);
  zcopy(zt2, zv2);
  zcopy(zt3, zv3);
  zsadd(*zz, 1l, &zt);
  zcopy(zt, zz);
  goto L2;
}

void NUDUPL(verylong zL, struct zform zform1,
            struct zform *zform2)
{
  verylong zA = 0, zB = 0, zC = 0, za = 0, zb = 0;
  verylong zc = 0, zd = 0, ze = 0, zg = 0, zr = 0;
  verylong zs = 0, zt = 0, zu = 0, zv = 0, zz = 0;
  verylong zC1 = 0, za2 = 0, zb2 = 0, zc2 = 0;
  verylong zd1 = 0, zv2 = 0, zv3 = 0;

  zcopy(zform1.za, &za);
  zcopy(zform1.zb, &zb);
  zcopy(zform1.zc, &zc);
  zexteucl(zb, &zu, za, &zv, &zd1);
  zdiv(za, zd1, &zA, &zr);
  zdiv(zb, zd1, &zB, &zr);
  zcopy(zc, &zr);
  znegate(&zr);
  zmulmod(zr, zu, zA, &zC);
  zsub(zA, zC, &zC1);
  if (zcompare(zC1, zC) < 0) {
    zcopy(zC1, &zC);
    znegate(&zC);
  }
  PARTEUCL(zL, zA, zC, &zv2, &zv3, &zd, &zz);
  if (zscompare(zz, 0l) == 0) {
    zmul(zB, zv3, &zr);
    zadd(zr, zc, &zs);
    zdiv(zs, zd, &zg, &zr);
    zsq(zd, &za2);
    zsq(zv3, &zc2);
    zadd(zd, zv3, &zr);
    zsq(zr, &zs);
    zadd(zb, zs, &zt);
    zsub(zt, za2, &zr);
    zsub(zr, zc2, &zb2);
    zmul(zg, zd1, &zr);
    zadd(zc2, zr, &zs);
    zcopy(zs, &zc2);
  }
  else {
    zmul(zc, zv, &zr);
    zmul(zB, zd, &zs);
    zadd(zr, zs, &zt);
    zdiv(zt, zA, &ze, &zr);
    zmul(ze, zv2, &zr);
    zsub(zr, zB, &zs);
    zdiv(zs, zv, &zg, &zt);
    zmul(zv, zg, &zs);
    zadd(zr, zs, &zb2);
    if (zscompare(zd1, 1l) > 0) {
      zmul(zd1, zb2, &zr);
      zcopy(zr, &zb2);
      zmul(zd1, zv, &zr);
      zcopy(zr, &zv);
      zmul(zd1, zv2, &zr);
      zcopy(zr, &zv2);
    }
    zsq(zd, &za2);
    zsq(zv3, &zc2);
    zadd(zd, zv3, &zr);
    zsq(zr, &zs);
    zadd(zb2, zs, &zr);
    zsub(zr, za2, &zs);
    zsub(zs, zc2, &zb2);
    zmul(ze, zv, &zr);
    zadd(za2, zr, &zs);
    zcopy(zs, &za2);
    zmul(zg, zv2, &zr);
    zadd(zc2, zr, &zs);
    zcopy(zs, &zc2);
  }
  reduce(&za2, &zb2, &zc2);
  zcopy(za2, &zform2->za);
  zcopy(zb2, &zform2->zb);
  zcopy(zc2, &zform2->zc);
  zfree(&zA);
  zfree(&zB);
  zfree(&zC);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&ze);
  zfree(&zg);
  zfree(&zr);
  zfree(&zs);
  zfree(&zt);
  zfree(&zu);
  zfree(&zv);
  zfree(&zz);
  zfree(&zC1);
  zfree(&za2);
  zfree(&zb2);
  zfree(&zc2);
  zfree(&zd1);
  zfree(&zv2);
  zfree(&zv3);
}

void composition(struct zform zform1, struct zform zform2,
                 struct zform *zform3)
{
  verylong za1 = 0, za2 = 0, zb1 = 0, zb2 = 0;
  verylong zc1 = 0, zc2 = 0, zd1 = 0, zv1 = 0;
  verylong zv2 = 0, zx1 = 0, zx2 = 0, zy1 = 0;
  verylong zy2 = 0;
  verylong za = 0, zb = 0, zc = 0, zd = 0, zn = 0;
  verylong zr = 0, zs = 0, zu = 0, zv = 0;

  zcopy(zform1.za, &za1);
  zcopy(zform1.zb, &zb1);
  zcopy(zform1.zc, &zc1);
  zcopy(zform2.za, &za2);
  zcopy(zform2.zb, &zb2);
  zcopy(zform2.zc, &zc2);
  if (zcompare(za1, za2) > 0) {
    /* excange forms 1 and 2 */
    zcopy(zform2.za, &za1);
    zcopy(zform2.zb, &zb1);
    zcopy(zform2.zc, &zc1);
    zcopy(zform1.za, &za2);
    zcopy(zform1.zb, &zb2);
    zcopy(zform1.zc, &zc2);
  }
  zadd(zb1, zb2, &za);
  zrshift(za, 1l, &zs);
  zsub(zb2, zs, &zn);
  zmod(za2, za1, &zr);
  if (zscompare(zr, 0l) == 0) {
    zzero(&zy1);
    zcopy(za1, &zd);
  }
  else {
    zexteucl(za2, &zu, za1, &zv, &zd);
    zcopy(zu, &zy1);
  }
  zmod(zs, zd, &zr);
  if (zscompare(zr, 0l) == 0) {
    zone(&zy2);
    znegate(&zy2);
    zzero(&zx2);
    zcopy(zd, &zd1);
  }
  else {
    zexteucl(zs, &zu, zd, &zv, &zd1);
    zcopy(zu, &zx2);
    zcopy(zv, &zy2);
    znegate(&zy2);
  }
  zdiv(za1, zd1, &zv1, &zr);
  zdiv(za2, zd1, &zv2, &zr);
  zmul(zy1, zy2, &za);
  zmul(za, zn, &zb);
  zmul(zx2, zc2, &zc);
  zsubmod(zb, zc, zv1, &zr);
  zmul(zv2, zr, &za);
  zlshift(za, 1l, &zb);
  zmul(zv1, zv2, &zform3->za);
  zadd(zb2, zb, &zform3->zb);
  zmul(zc2, zd1, &zb);
  zadd(zb2, za, &zc);
  zmul(zc, zr, &zd);
  zadd(zb, zd, &za);
  zdiv(za, zv1, &zform3->zc, &zr);
  reduce(&zform3->za, &zform3->zb, &zform3->zc);
  zfree(&za1);
  zfree(&za2);
  zfree(&zb1);
  zfree(&zb2);
  zfree(&zc1);
  zfree(&zc2);
  zfree(&zd1);
  zfree(&zv1);
  zfree(&zv2);
  zfree(&zx1);
  zfree(&zx2);
  zfree(&zy1);
  zfree(&zy2);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&zn);
  zfree(&zr);
  zfree(&zs);
  zfree(&zu);
  zfree(&zv);
}

void zsquare_root(verylong za, verylong zp, verylong *zr)
{
  long i, s = 0, t;
  verylong zb = 0, zc = 0, zd = 0, ze = 0, zi = 0;
  verylong zq = 0, zs = 0, zt = 0, zx = 0, zy = 0;

  if (zcompare(za, zp) < 0) zcopy(za, &ze);
  else
    zmod(za, zp, &ze);
  if (zscompare(ze, 0l) == 0)
    zzero(zr);
  else if (zjacobi(ze, zp) == - 1)
    zzero(zr);
  else {
    do
      zrandomb(zp, &zb);
    while (zscompare(zb, 0l) == 0 || zjacobi(zb, zp) != - 1);
    zsadd(zp, - 1l, &zq);
    zcopy(zq, &zy);
    do {
      zrshift(zy, 1l, &zt);
      s++;
      t = zodd(zt);
      zcopy(zt, &zy);
    } while (!t);
    zinvmod(ze, zp, &zi);
    zexpmod(zb, zt, zp, &zc);
    zsadd(zt, 1l, &zs);
    zrshift(zs, 1l, &zx);
    zexpmod(ze, zx, zp, zr);
    for (i = 1; i < s; i++) {
      zsq(*zr, &zs);
      zmul(zs, zi, &zx);
      zsexpmod(zx, pow(2, s - i - 1), zp, &zd);
      if (zcompare(zd, zq) == 0) {
        zmulmod(*zr, zc, zp, &zx);
        zcopy(zx, zr);
      }
      zmulmod(zc, zc, zp, &zx);
      zcopy(zx, &zc);
    }
  }
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&ze);
  zfree(&zi);
  zfree(&zq);
  zfree(&zs);
  zfree(&zt);
  zfree(&zx);
}

void zsquare_roots(verylong za, verylong zq,
                   verylong *zx1, verylong *zx2,
                   verylong *zy1, verylong *zy2)
{
  long i, n, r;
  verylong zb = 0, zc = 0, zd = 0, ze = 0, zg = 0;
  verylong zn = 0, zp = 0, zr = 0, zs = 0, zx = 0;
  verylong zy = 0, zz = 0;

  zintoz(4l, &zp);
  zcopy(za, &zn);
  n = zsmod(zn, 4l);
  for (i = 1; i < 4; i++)
    if ((i * i) % 4 == n) r = i;
  zintoz(r, &zr);
  if (zcompare(za, zq) < 0)
    zcopy(za, &zn);
  else
    zmod(za, zq, &zn);
  zsquare_root(zn, zq, &zs);
  zexteucl(zp, &zc, zq, &zd, &zg);
  if (zscompare(zg, 1l) == 0) {
    zmul(zp, zq, &zn);
    zmul(zr, zq, &zb);
    zmul(zs, zc, &ze);
    zmul(ze, zp, &zz);
    zaddmod(zb, zz, zn, &zx);
    zmul(zr, zq, &zb);
    zmul(zs, zc, &ze);
    zmul(ze, zp, &zz);
    zsubmod(zb, zz, zn, &zy);
    zmod(zx, zn, zx1);
    znegate(&zx);
    zmod(zx, zn, zx2);
    zmod(zy, zn, zy1);
    znegate(&zy);
    zmod(zy, zn, zy2);
  }
  else {
    zzero(zx1);
    zzero(zx2);
    zzero(zy1);
    zzero(zy2);
  }
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&ze);
  zfree(&zg);
  zfree(&zn);
  zfree(&zp);
  zfree(&zr);
  zfree(&zs);
  zfree(&zx);
  zfree(&zy);
  zfree(&zz);
}

void zform_pow(long e, struct zform x,
               struct zform *y,
               verylong zL)
{
  struct zform s, t, u;

  init_zform(&s);
  init_zform(&t);
  init_zform(&u);
  if (e == 1) copy_zform(x, y);
  else {
    copy_zform(x, &t);
    zintoz(1l, &s.za);
    zcopy(x.zb, &s.zb);
    znegate(&s.zb);
    zmul(x.za, x.zc, &s.zc);
    reduce(&s.za, &s.zb, &s.zc);
    while (e > 0) {
      if (e & 1) {
        composition(s, t, &u);
        copy_zform(u, &s);
      }
      e >>= 1;
      if (e > 0) {
        NUDUPL(zL, t, &u);
        copy_zform(u, &t);
      }
    }
    copy_zform(s, y);
  }
  free_zform(&s);
  free_zform(&t);
  free_zform(&u);
}

int baby_step(verylong zD)
{
  double B, B1, C, C1, product = 1.0;
  double D = fabs(zdoub(zD)), h, sP, Q, Q1;
  int error, flag;
  long b = 10, c, e, i, j, n, p, q, q1, r;
  long L = pow(D, 0.25), P = max(262144L, L);
  struct zform *f, g, t, unit, x[128],y, z;
  verylong zE = 0, zL = 0, za = 0, zb = 0, zc = 0;
  verylong zp = 0, zq = 0, zr = 0;
  verylong zx2 = 0, zy1 = 0, zy2 = 0;

  zintoz(L, &zL);
  init_zform(&g);
  init_zform(&t);
  init_zform(&unit);
  for (r = 0; r < 128; r++)
    init_zform(&x[r]);
  init_zform(&y);
  init_zform(&z);
  zpstart2();
  do {
    p = zpnext();
    zintoz(p, &zp);
    j = zjacobi(zD, zp);
    if (j != 0) product /= 1.0 - j / (double) p;
  } while (p <= P);
  Q = floor(sqrt(D) * product / (3.0 * M_PI));
  sP = 1.0 / (2.0 * sqrt(P));
  B = floor(Q * (1.0 + sP));
  C =  ceil(Q * (1.0 - sP));
  e = 1, c = 0, B1 = B, C1 = C, Q1 = Q;
  zpstart();
  L3:
  do zintoz(zpnext(), &zp); while (zjacobi(zD, zp) != 1);
  zlshift(zp, 2l, &zq);
  zmod(zD, zq, &zE);
  zsquare_roots(zE, zp, &zb, &zx2, &zy1, &zy2);
  if (zscompare(zb, 1l) == 0) zcopy(zx2, &zb);
  zcopy(zp, &g.za);
  zcopy(zb, &g.zb);
  zsq(zb, &za);
  zsub(za, zD, &zc);
  zdiv(zc, zq, &g.zc, &zr);
  c++, q = ceil(sqrt((B1 - C1) / 2.0));
  if (q == 0) q1 = 1; else q1 = q;
  zone(&x[0].za);
  zcopy(g.zb, &x[0].zb);
  znegate(&x[0].zb);
  zmul(g.za, g.zc, &x[0].zc);
  copy_zform(x[0], &unit);
  reduce(&unit.za, &unit.zb, &unit.zc);
  copy_zform(unit, &x[0]);
  zform_pow(e, g, &x[1], zL);
  flag = 0;
  for (r = 2; r < q; r++) {
    composition(x[1], x[r - 1], &x[r]);
    if (compare_zform(&x[r], &unit) == 0) {
      flag = 1;
      n = r;
      break;
    }
  }
  if (flag == 1) goto L7;
  /* sort the forms using a crude algorithm */
  for (i = 0; i < q1 - 1; i++) {
    for (j = i + 1; j < q1; j++) {
      if (compare_zform(&x[i], &x[j]) > 0) {
        copy_zform(x[i], &t);
        copy_zform(x[j], &x[i]);
        copy_zform(t, &x[j]);
      }
    }
  }
  composition(x[1], x[q1 - 1], &y);
  NUDUPL(zL, y, &z);
  copy_zform(z, &y);
  zform_pow(Q1, x[1], &z, zL);
  n = Q1;
  L5:
  f = search_zform(q1, z, x);
  if (f != 0) {
    error = 0;
    n -= f - x;
    goto L7;
  }
  znegate(&z.zb);
  f = search_zform(q1, z, x);
  if (f != 0) {
    error = 0;
    n += f - x;
    goto L7;
  }
  znegate(&z.zb);
  composition(y, z, &t);
  copy_zform(t, &z);
  n += 2 * q;
  error = 0;
  if ((double) n <= B1) goto L5;
  else error = - 1;
  L7:
  if (!error) {
    zpstart2();
    p = zpnext();
    L7a:
    if (n > 1) {
      if (n % p == 0) {
        zform_pow(n / p, x[1], &t, zL);
        if (compare_zform(&t, &unit) == 0) {
          n /= p;
          goto L7a;
        }
      }
    }
    L8:
    e *= n;
    if (e > B - C) h = e * floor(B / e);
    else  if (c >= b) error = - 1;
    else {
      B1 = floor(B1 / n);
      C1 = ceil(C1 / n);
      goto L3;
    }
  }
  if (error) h = error;
  for (r = 0; r < 128; r++)
    free_zform(&x[r]);
  free_zform(&g);
  free_zform(&t);
  free_zform(&unit);
  free_zform(&y);
  free_zform(&z);
  zfree(&zE);
  zfree(&zL);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zp);
  zfree(&zq);
  zfree(&zr);
  zfree(&zx2);
  zfree(&zy1);
  zfree(&zy2);
  return h;
}

int main(void)
{
  long D[5] = {- 27l, - 51, - 83, - 155, - 227}, i;
  verylong zD = 0;

  for (i = 0; i < 5; i++) {
    zintoz(D[i], &zD);
    printf("D = %4ld h(D) = %d\n", D[i], baby_step(zD));
  }
  zfree(&zD);
  return 0;
}