/*
  Author:  Pate Williams (c) 1997

  Algorithm 7.5.1 (Reduction of an Elliptic Curve
  Modulo p). See "A Course in Computational Alge-
  braic Number Theory" by Henri Cohen page 407.
*/

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

#define BOUND 1000000l
#define DEBUG
#define POLY_SIZE 8192l

struct factor {long expon, prime;};

void zpoly_copy(long da, verylong *za, verylong *zb, long *db)
{
  long i;

  *db = da;
  for (i = 0; i <= da; i++) zcopy(za[i], &zb[i]);
}

void zpoly_mul(long m, long n, verylong *za, verylong *zb,
               verylong *zc, long *p)
{
  long i, j, k;
  verylong zd = 0, zai = 0, zbk = 0, zsum = 0, zterm = 0;

  *p = m + n;
  for (k = 0; k <= *p; k++) {
    zzero(&zsum);
    for (i = 0; i <= k; i++) {
      j = k - i;
      if (i > m) zzero(&zai);
      else zcopy(za[i], &zai);
      if (j > n) zzero(&zbk);
      else zcopy(zb[j], &zbk);
      zmul(zai, zbk, &zterm);
      zcopy(zsum, &zd);
      zadd(zterm, zd, &zsum);
    }
    zcopy(zsum, &zc[k]);
  }
  zfree(&zd);
  zfree(&zai);
  zfree(&zbk);
  zfree(&zsum);
  zfree(&zterm);
}

void zpoly_div(long m, long n, verylong *zu, verylong *zv,
               verylong *zq, verylong *zr, long *p, long *s)
{
  long j, jk, k, nk;
  verylong za = 0, zb = 0, zvn = 0;

  zcopy(zv[n], &zvn);
  for (j = 0; j <= m; j++)
    zcopy(zu[j], &zr[j]);
  if (m < n) {
    *p = 0, *s = m;
    zzero(&zq[0]);
  }
  else {
    *p = m - n, *s = n - 1;
    for (k = *p; k >= 0; k--) {
      nk = n + k;
      zsexp(zvn, k, &za);
      zmul(zr[nk], za, &zq[k]);
      for (j = nk - 1; j >= 0; j--) {
        jk = j - k;
        if (jk >= 0) {
          zmul(zvn, zr[j], &za);
          zmul(zr[nk], zv[jk], &zb);
          zsub(za, zb, &zr[j]);
        }
        else {
          zcopy(zr[j], &za);
          zmul(zvn, za, &zr[j]);
        }
      }
    }
    while (*p > 0 && zscompare(zq[*p], 0l) == 0) *p = *p - 1;
    while (*s > 0 && zscompare(zr[*s], 0l) == 0) *s = *s - 1;
  }
  zfree(&za);
  zfree(&zb);
  zfree(&zvn);
}

void zpoly_mod(verylong zp, verylong *za, long *da)
{
  long i;
  verylong zb = 0;

  for (i = 0; i <= *da; i++) {
    zmod(za[i], zp, &zb);
    zcopy(zb, &za[i]);
  }
  while (*da > 0 && zscompare(za[*da], 0l) == 0) *da = *da - 1;
  zfree(&zb);
}

void zpoly_pow(long degreeA, long degreem, verylong zn,
               verylong zp, verylong *zA,
               verylong *zm, verylong *zs,
               long *ds)
{
  long dP, dq, dx = degreeA, i;
  verylong za = 0, zb = 0, zP[POLY_SIZE], zq[POLY_SIZE],
           zx[POLY_SIZE], zy[POLY_SIZE];

  for (i = 0; i < POLY_SIZE; i++)
    zP[i] = zq[i] = zx[i] = zy[i] = 0;
  *ds = 0;
  zcopy(zn, &za);
  zone(&zs[0]);
  for (i = 0; i <= dx; i++) zcopy(zA[i], &zx[i]);
  while (zscompare(za, 0l) > 0) {
    if (zodd(za)) {
      /* s = (s * x) % m; */
      zpoly_mul(*ds, dx, zs, zx, zP, &dP);
      zpoly_div(dP, degreem, zP, zm, zq, zs, &dq, ds);
      zpoly_mod(zp, zs, ds);
    }
    zcopy(za, &zb);
    zrshift(zb, 1l, &za);
    if (zscompare(za, 0l) > 0) {
      /* x = (x * x) % m; */
      for (i = 0; i <= dx; i++) zcopy(zx[i], &zy[i]);
      zpoly_mul(dx, dx, zx, zy, zP, &dP);
      zpoly_div(dP, degreem, zP, zm, zq, zx, &dq, &dx);
      zpoly_mod(zp, zx, &dx);
    }
  }
  zfree(&za);
  zfree(&zb);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zP[i]);
    zfree(&zq[i]);
    zfree(&zx[i]);
    zfree(&zy[i]);
  }
}

void zpoly_sub(long da, long db, verylong *za, verylong *zb,
               verylong *zc, long *dc)
{
  long i;
  verylong zz = 0;

  zzero(&zz);
  if (da >= db) {
    for (i = 0; i <= db; i++)
      zsub(za[i], zb[i], &zc[i]);
    for (i = db + 1; i <= da; i++)
      zcopy(za[i], &zc[i]);
    *dc = da;
  }
  else {
    for (i = 0; i <= da; i++)
      zsub(za[i], zb[i], &zc[i]);
    for (i = da + 1; i <= db; i++)
      zsub(zz, zb[i], &zc[i]);
    *dc = db;
  }
  zfree(&zz);
}

void zpoly_gcd(long degreeA, long degreeB, verylong zp,
               verylong *zA, verylong *zB, verylong *za,
               long *da)
{
  int nonzero = 0, zero;
  long db, dq, dr, i;
  verylong zc = 0;
  verylong zb[POLY_SIZE], zq[POLY_SIZE], zr[POLY_SIZE];

  for (i = 0; i < POLY_SIZE; i++)
    zb[i] = zq[i] = zr[i] = 0;
  if (degreeA > degreeB) {
    *da = degreeA;
    db = degreeB;
    for (i = 0; i <= *da; i++) zcopy(zA[i], &za[i]);
    for (i = 0; i <= db; i++) zcopy(zB[i], &zb[i]);
  }
  else {
    *da = degreeB;
    db = degreeA;
    for (i = 0; i <= *da; i++) zcopy(zB[i], &za[i]);
    for (i = 0; i <= db; i++) zcopy(zA[i], &zb[i]);
  }
  for (i = 0; i <= db && !nonzero; i++)
    nonzero = zscompare(zb[i], 0l) != 0;
  while (nonzero) {
    zpoly_div(*da, db, za, zb, zq, zr, &dq, &dr);
    for (i = 0; i <= dr; i++) {
      zcopy(zr[i], &zc);
      zmod(zc, zp, &zr[i]);
    }
    zero = 1;
    for (i = dr; i >= 0 && zero; i--) {
      zero = zscompare(zr[i], 0l) == 0;
      if (zero && dr > 0) dr--;
    }
    for (i = 0; i <= db; i++) zcopy(zb[i], &za[i]);
    *da = db;
    for (i = 0; i <= dr; i++) zcopy(zr[i], &zb[i]);
    db = dr;
    nonzero = 0;
    for (i = 0; i <= db && !nonzero; i++)
      nonzero = zscompare(zb[i], 0l) != 0;
  }
  zfree(&zc);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zb[i]);
    zfree(&zq[i]);
    zfree(&zr[i]);
  }
}

void zpoly_print(long da, verylong *za)
{
  long i;

  for (i = da; i >= 0; i--) {
    zwrite(za[i]);
    printf(" ");
  }
  printf("\n");
}

void zpoly_ext_euclid(long dg, long dh, verylong zp, verylong *zg,
                      verylong *zh, verylong *zs,
                      verylong *zt, verylong *zd,
                      long *ds, long *dt, long *dd)
{
  long da, dq, dr, ds1 = 0, ds2 = 0, dt1 = 0, dt2 = 0, i;
  verylong za[POLY_SIZE], zb[POLY_SIZE];
  verylong zq[POLY_SIZE], zr[POLY_SIZE];
  verylong zs1[POLY_SIZE], zs2[POLY_SIZE];
  verylong zt1[POLY_SIZE], zt2[POLY_SIZE];

  if (dh == 0 && zscompare(zh[0], 0l) == 0) {
    zpoly_copy(dg, zg, zd, dd);
    *ds = *dt = 0;
    zone(&zs[0]);
    zzero(&zt[0]);
  }
  for (i = 0; i < POLY_SIZE; i++) {
    za[i] = zb[i] = zq[i] = zr[i] = 0;
    zs1[i] = zs2[i] = zt1[i] = zt2[i] = 0;
  }
  zone(&zs2[0]);
  zzero(&zs1[0]);
  zzero(&zt2[0]);
  zone(&zt1[0]);
  while (dh != 0 || zscompare(zh[0], 0l) != 0) {
    zpoly_div(dg, dh, zg, zh, zq, zr, &dq, &dr);
    zpoly_mod(zp, zq, &dq);
    zpoly_mod(zp, zr, &dr);
    zpoly_mul(dq, ds1, zq, zs1, za, &da);
    zpoly_sub(ds2, da, zs2, za, zs, ds);
    zpoly_mul(dq, dt1, zq, zt1, za, &da);
    zpoly_sub(dt2, da, zt2, za, zt, dt);
    zpoly_mod(zp, zs, ds);
    zpoly_mod(zp, zt, dt);
    zpoly_copy(dh, zh, zg, &dg);
    zpoly_copy(dr, zr, zh, &dh);
    zpoly_copy(ds1, zs1, zs2, &ds2);
    zpoly_copy(*ds, zs, zs1, &ds1);
    zpoly_copy(dt1, zt1, zt2, &dt2);
    zpoly_copy(*dt, zt, zt1, &dt1);
    #ifdef DEBUG
    printf("q  = "); zpoly_print(dq, zq);
    printf("r  = "); zpoly_print(dr, zr);
    printf("s  = "); zpoly_print(*ds, zs);
    printf("t  = "); zpoly_print(*dt, zt);
    printf("g  = "); zpoly_print(dg, zg);
    printf("h  = "); zpoly_print(dh, zh);
    printf("s2 = "); zpoly_print(ds2, zs2);
    printf("s1 = "); zpoly_print(ds1, zs1);
    printf("t2 = "); zpoly_print(dt2, zt2);
    printf("t1 = "); zpoly_print(dt1, zt1);
    #endif
  }
  zpoly_copy(dg, zg, zd, dd);
  zpoly_copy(ds2, zs2, zs, ds);
  zpoly_copy(dt2, zt2, zt, dt);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&za[i]);
    zfree(&zb[i]);
    zfree(&zq[i]);
    zfree(&zr[i]);
    zfree(&zs1[i]);
    zfree(&zs2[i]);
    zfree(&zt1[i]);
    zfree(&zt2[i]);
  }
}

void Recurse(int degreeA, verylong zp, verylong *zA,
             verylong *zroot, long *rootSize)
{
  long dd, degreeB, dq, dr, du = 1, i;
  verylong zD = 0, za = 0, zb = 0, zc = 0, ze = 0;
  verylong zn = 0;
  verylong zB[POLY_SIZE], zd[POLY_SIZE];
  verylong zq[POLY_SIZE], zr[POLY_SIZE];
  verylong zu[2];

  for (i = 0; i < POLY_SIZE; i++)
    zB[i] = zd[i] = zq[i] = zr[i] = 0;
  zu[0] = zu[1] = 0;
  if (degreeA != 0) {
    if (degreeA == 1) {
      if (zscompare(zA[1], 0l) != 0) {
         zinvmod(zA[1], zp, &za);
         zmul(zA[0], za, &zb);
         znegate(&zb);
         zmod(zb, zp, &zroot[*rootSize]);
      }
      *rootSize = *rootSize + 1;
    }
    else if (degreeA == 2) {
      zsq(zA[1], &za);
      zmul(zA[0], zA[2], &zb);
      zlshift(zb, 2l, &zc);
      zsub(za, zc, &zb);
      zmod(zb, zp, &zD);
      zsqrtmod(zD, zp, &ze);
      zlshift(zA[2], 1l, &za);
      zinvmod(za, zp, &zD);
      zsub(ze, zA[1], &za);
      zmulmod(za, zD, zp, &zroot[*rootSize]);
      *rootSize = *rootSize + 1;
      znegate(&zA[1]);
      znegate(&ze);
      zadd(zA[1], ze, &za);
      zmulmod(za, zD, zp, &zroot[*rootSize]);
      *rootSize = *rootSize + 1;
    }
    else {
      zsadd(zp, - 1l, &za);
      zrshift(za, 1l, &zn);
      do {
        zrandomb(zp, &za);
        zcopy(za, &zu[0]);
        zone(&zu[1]);
        zpoly_pow(du, degreeA, zn, zp, zu, zA, zd, &dd);
        zsadd(zd[0], - 1l, &za);
        zcopy(za, &zd[0]);
        zpoly_gcd(dd, degreeA, zp, zd, zA, zB, &degreeB);
      } while (degreeB == 0 || degreeB == degreeA);
      Recurse(degreeB, zp, zB, zroot, rootSize);
      zpoly_div(degreeA, degreeB, zA, zB, zq, zr, &dq, &dr);
      zpoly_mod(zp, zq, &dq);
      Recurse(dq, zp, zq, zroot, rootSize);
    }
  }
  zfree(&zD);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&ze);
  zfree(&zn);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zB[i]);
    zfree(&zd[i]);
    zfree(&zq[i]);
    zfree(&zr[i]);
  }
  zfree(&zu[0]);
  zfree(&zu[1]);
}

void FindRootsModuloAPrime(long degreeP, verylong zp,
                           verylong *zP, verylong *zroot,
                           long *rootSize)
{
  long degreeA, dy, i, j;
  verylong za =0, zt = 0;
  verylong zA[POLY_SIZE], zx[POLY_SIZE], zy[POLY_SIZE];

  for (i = 0; i < POLY_SIZE; i++)
    zA[i] = zx[i] = zy[i] = 0;
  *rootSize = 0;
  zzero(&zx[0]);
  zone(&zx[1]);
  zpoly_pow(1, degreeP, zp, zp, zx, zP, zy, &dy);
  zsadd(zy[1], - 1l, &za);
  zcopy(za, &zy[1]);
  zpoly_gcd(dy, degreeP, zp, zy, zP, zA, &degreeA);
  if (zscompare(zA[0], 0l) == 0) {
    zzero(&zroot[*rootSize]);
    *rootSize = *rootSize + 1;
    for (i = 0; i < degreeA; i++)
      zcopy(zA[i + 1], &zA[i]);
    degreeA--;
  }
  Recurse(degreeA, zp, zA, zroot, rootSize);
  /* sort the roots using selection sort */
  for (i = 0; i < *rootSize - 1; i++) {
    for (j = i + 1; j < *rootSize; j++) {
      if (zcompare(zroot[i], zroot[j]) > 0) {
        zcopy(zroot[i], &zt);
        zcopy(zroot[j], &zroot[i]);
        zcopy(zt, &zroot[j]);
      }
    }
  }
  zfree(&zt);
  for (i = 0; i < POLY_SIZE; i++) {
    zfree(&zA[i]);
    zfree(&zx[i]);
    zfree(&zy[i]);
  }
}

int trial_division(struct factor *f, verylong zn, long *number)
{
  int pr = 0;
  long e, m, p;
  verylong za = 0, zb = 0;

  *number = 0;
  zcopy(zn, &za);
  zabs(&za);
  zpstart2();
  do {
    p = zpnext();
    if (zsmod(za, p) == 0) {
      e = 0;
      do {
        e++;
        zsdiv(za, p, &zb);
        zcopy(zb, &za);
        m = zsmod(za, p);
        pr = zprobprime(za, 5l);
        if (zscompare(za, 1l) == 0) pr = 1;
        if (zscompare(za, p) == 0) pr = 0;
      } while (!pr && m == 0);
      f[*number].expon = e;
      f[*number].prime = p;
      *number = *number + 1;
    }
  } while (!pr && p <= BOUND);
  if (pr && zscompare(za, 1l) != 0 && zscompare(za, LONG_MAX) <= 0) {
    f[*number].expon = 1;
    f[*number].prime = ztoint(za);
    *number = *number + 1;
  }
  zfree(&za);
  zfree(&zb);
  return pr;
}

long gcd(long x, long y)
{
  long g;

  if (x < 0) x = - x;
  if (y < 0) y = - y;
  g = y;
  while (x > 0)
    g = x, x = y % x, y = g;
  return g;
}

int reduction(long p, verylong *zA, verylong *zr,
              verylong *zs, verylong *zt,
              verylong *zu, long *c, long *f, long *nu)
{
  int found, j, kodaira;
  long dindex, dnumber, jindex, jnumber, i, k, vpd, vpj;
  long rootSize;
  struct factor dfactor[32], jfactor[32];
  verylong za = 0, zb = 0, zc = 0, zd = 0, zj = 0;
  verylong zp = 0;
  verylong za0 = 0, za1 = 0, za2 = 0, za3 = 0;
  verylong za4 = 0, za6 = 0, zb2 = 0, zb4 = 0;
  verylong zb6 = 0, zb8 = 0, zc4 = 0, zc6 = 0;
  verylong za2p = 0, za3p = 0, zdelta = 0;
  verylong zP[4] = {0, 0, 0, 0};
  verylong zroot[4] = {0, 0, 0, 0};

  *nu = 0;
  zintoz(p, &zp);
  zcopy(zA[0], &za0);
  zcopy(zA[1], &za1);
  zcopy(zA[2], &za2);
  zcopy(zA[3], &za3);
  zcopy(zA[4], &za4);
  zcopy(zA[6], &za6);
  zsq(za1, &za);
  zlshift(za2, 2l, &zb);
  zadd(za, zb, &zb2);
  zmul(za1, za3, &za);
  zlshift(za4, 1l, &zb);
  zadd(za, zb, &zb4);
  zsq(za3, &za);
  zlshift(za6, 2l, &zb);
  zadd(za, zb, &zb6);
  zsq(za1, &za);
  zmul(za, za6, &zb);
  zmul(za2, za6, &zc);
  zlshift(zc, 2l, &zd);
  zadd(zb, zd, &za);
  zmul(za1, za3, &zb);
  zmul(zb, za4, &zc);
  zsub(za, zc, &zb);
  zsq(za3, &za);
  zmul(za, za2, &zc);
  zadd(zb, zc, &za);
  zsq(za4, &zb);
  zsub(za, zb, &zb8);
  zsq(zb2, &za);
  zsmul(zb4, 24l, &zb);
  zsub(za, zb, &zc4);
  zsexp(zb2, 3l, &za);
  znegate(&za);
  zsmul(zb2, 36l, &zb);
  zmul(zb, zb4, &zc);
  zadd(za, zc, &zb);
  zsmul(zb6, 216l, &za);
  zsub(zb, za, &zc6);
  zsq(zb2, &za);
  zmul(za, zb8, &zb);
  znegate(&zb);
  zsq(zb4, &za);
  zmul(za, zb4, &zd);
  zlshift(zd, 3l, &zc);
  znegate(&zc);
  zadd(zb, zc, &za);
  zsq(zb6, &zb);
  zsmul(zb, 27l, &zc);
  znegate(&zc);
  zadd(za, zc, &zb);
  zmul(zb2, zb4, &za);
  zmul(za, zb6, &zc);
  zsmul(zc, 9l, &zd);
  zadd(zb, zd, &zdelta);
  zsq(zc4, &za);
  zmul(za, zc4, &zb);
  zdiv(zb, zdelta, &zj, &zc);
  if (!trial_division(dfactor, zdelta, &dnumber)) {
    fprintf(stdout, "*error*\ncould not factor delta\n");
    exit(1);
  }
  if (!trial_division(jfactor, zb, &jnumber)) {
    fprintf(stdout, "*error*\ncould not factor j\n");
    exit(1);
  }
  #ifdef DEBUG
  printf("b2 = "); zwriteln(zb2);
  printf("b4 = "); zwriteln(zb4);
  printf("b6 = "); zwriteln(zb6);
  printf("b8 = "); zwriteln(zb8);
  printf("c4 = "); zwriteln(zc4);
  printf("c6 = "); zwriteln(zc6);
  printf("n  = "); zwriteln(zb);
  printf("d  = "); zwriteln(zdelta);
  printf("j  = "); zwriteln(zj);
  printf("factorization of n\n");
  for (i = 0; i < jnumber; i++)
    printf("%10ld %10ld\n", jfactor[i].prime, jfactor[i].expon);
  printf("factorization of d\n");
  for (i = 0; i < dnumber; i++)
    printf("%10ld %10ld\n", dfactor[i].prime, dfactor[i].expon);
  #endif
  dindex = - 1;
  found = 0;
  for (i = 0; i < dnumber && !found; i++) {
    found = p == dfactor[i].prime;
    if (found) dindex = i;
  }
  jindex = - 1;
  found = 0;
  for (i = 0; i < jnumber && !found; i++) {
    found = p == jfactor[i].prime;
    if (found) jindex = i;
  }
  if (dindex != - 1) vpd = dfactor[dindex].expon;
  else vpd = 0;
  if (jindex != - 1) vpj = jfactor[jindex].expon;
  else vpj = 0;
  vpj -= vpd;
  if (vpj < 0) k = vpd + vpj; else k = vpd;
  if (k < 12) {
    zone(zu);
    zzero(zr);
    zzero(zs);
    zzero(zt);
  }
  else {
    zsexp(zp, k / 12, zu);
    if (zodd(za1)) {
      zsub(*zu, za1, &za);
      zrshift(za, 1l, zs);
    }
    else {
      zrshift(za1, 1l, zs);
      znegate(zs);
    }
    zmul(*zs, za1, &za);
    zsub(za2, za, &zb);
    zsq(*zs, &za);
    zsub(zb, za, &za2p);
    switch (zsmod(za2p, 3l)) {
      case 0 :
        zsdiv(za2p, 3l, zr);
        znegate(zr);
        break;
      case 1 :
        zsq(*zu, &za);
        zsub(za, za2p, &zb);
        zsdiv(zb, 3l, zr);
        break;
      case 2 :
        zsq(*zu, &za);
        znegate(&za);
        zsub(za, za2p, &zb);
        zsdiv(zb, 3l, zr);
        break;
    }
    zmul(*zr, za1, &za);
    zadd(za3, za, &za3p);
    if (zodd(za3p)) {
      zsq(*zu, &za);
      zmul(za, *zu, &zb);
      zsub(zb, za3p, zt);
    }
    else {
      zrshift(za3p, 1l, zt);
      znegate(zt);
    }
    k %= 12;
    zcopy(*zu, &za);
    if (zscompare(*zu, 0l) < 0) znegate(&za);
    zsexp(za, 12l, &zb);
    zdiv(zdelta, zb, &zc, &zd);
    zcopy(zc, &zdelta);
    zsexp(za, 4l, &zb);
    zdiv(zc4, zb, &zc, &zd);
    zcopy(zc, &zc4);
    zsexp(za, 6l, &zb);
    zdiv(zc6, zb, &zc, &zd);
    zcopy(zc, &zc6);
  }
  if (vpj < 0) {
    *nu = - vpj;
    if (k == 0) {
      *f = 1;
      zcopy(zc6, &za);
      znegate(&za);
      j = zjacobi(za, zp);
      if (j == 1) *c = 1;
      else if (j == - 1) *c = gcd(2, *nu);
      kodaira = 1;
    }
    else if (k == 6) {
      *f = 2;
      if (*nu & 1) {
        zsexp(zp, 9 + *nu, &za);
        zmul(zdelta, zc6, &zb);
        zdiv(zb, za, &zc, &zd);
        *c = 3 + zjacobi(zc, zp);
      }
      else {
        zsexp(zp, 6 + *nu, &za);
        zmul(zdelta, zc6, &zb);
        zdiv(zb, za, &zc, &zd);
        *c = 3 + zjacobi(zc, zp);
      }
      kodaira = 6;
    }
  }
  else {
    if (k == 0) *f = 0; else *f = 2;
    zsmul(zc6, 6l, &za);
    znegate(&za);
    zsq(zp, &zb);
    zdiv(za, zb, &zc, &zd);
    j = zjacobi(zc, zp);
    zone(&zP[3]);
    zzero(&zP[2]);
    zsmul(zc4, 3l, &za);
    znegate(&za);
    zrshift(za, 2l, &zb);
    zsq(zp, &za);
    zdiv(zb, za, &zP[1], &zc);
    zsexp(zp, 3l, &za);
    zrshift(zc6, 2l, &zb);
    znegate(&zb);
    zdiv(zb, za, &zP[0], &zc);
    FindRootsModuloAPrime(3, zp, zP, zroot, &rootSize);
    switch (k) {
      case 0 :
        *c = 1;
        kodaira = 0;
        break;
      case 2 :
        *c = 1;
        kodaira = 2;
        break;
      case 3 :
        *c = 2;
        kodaira = 3;
        break;
      case 4 :
        *c = 2 + j;
        kodaira = 4;
        break;
      case 6 :
        *c = 1 + rootSize;
        kodaira = 5;
        break;
      case 8 :
        zsmul(zc6, 6l, &za);
        znegate(&za);
        zsexp(zp, 4l, &zb);
        zdiv(za, zb, &zc, &zd);
        *c = 2 + zjacobi(zc, zp);
        kodaira = 9;
        break;
      case 9 :
        *c = 2;
        kodaira = 8;
        break;
      case 10 :
        *c = 1;
        kodaira = 7;
        break;
      default :
        fprintf(stderr, "*error*\nunknown k = %ld\n", k);
        exit(1);
        break;
    }
  }
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&zj);
  zfree(&zp);
  zfree(&za0);
  zfree(&za1);
  zfree(&za2);
  zfree(&za3);
  zfree(&za4);
  zfree(&za6);
  zfree(&zb2);
  zfree(&zb4);
  zfree(&zb6);
  zfree(&zb8);
  zfree(&zc4);
  zfree(&zc6);
  zfree(&za2p);
  zfree(&za3p);
  zfree(&zdelta);
  for (k = 0; k < 4; k++) {
    zfree(&zP[k]);
    zfree(&zroot[k]);
  }
  return kodaira;
}

int main(void)
{
  int r;
  long a[7] = {12, 8, 15, 11, 20, 0, 10}, c, f, i, nu, p = 5;
  verylong za[7] = {0, 0, 0, 0, 0, 0, 0};
  verylong zr = 0, zs = 0, zt = 0, zu = 0;

  for (i = 0; i < 7; i++) zintoz(a[i], &za[i]);
  r = reduction(p, za, &zr, &zs, &zt, &zu,
                &c, &f, &nu);
  printf("r = "); zwriteln(zr);
  printf("s = "); zwriteln(zs);
  printf("t = "); zwriteln(zt);
  printf("u = "); zwriteln(zu);
  printf("c = %ld\n", c);
  printf("f = %ld\n", f);
  printf("v = %ld\n", nu);
  printf("k = %d\n", r);
  zfree(&zr);
  zfree(&zs);
  zfree(&zt);
  zfree(&zu);
  for (i = 0; i < 7; i++) zfree(&za[i]);
  return 0;
}