/*
  Author:  Pate Williams (c) 1997

  Algorithm 1.5.2 (Cornacchia). See "A Course in
  Computational Algebraic Number Theory" by Henri
  Cohen pages 34 - 35. Let p be a prime number
  and 0 < d < p, this algorithm either outputs
  an integer solution (x, y) to the Diophantine
  equation x * x + d * y * y = p, or says that
  such a solution does not exist.
*/

#include <ctype.h>
#include <stdio.h>
#include "lip.h"

int square_test(verylong zn, verylong *zq)
{
  int square = 1;
  long q11[11], q63[63], q64[64], q65[65];
  long k, r, t;
  verylong zd = 0;

  for (k = 0; k < 11; k++) q11[k] = 0;
  for (k = 0; k < 6; k++) q11[(k * k) % 11] = 1;
  for (k = 0; k < 63; k++) q63[k] = 0;
  for (k = 0; k < 32; k++) q63[(k * k) % 63] = 1;
  for (k = 0; k < 64; k++) q64[k] = 0;
  for (k = 0; k < 32; k++) q64[(k * k) % 64] = 1;
  for (k = 0; k < 65; k++) q65[k] = 0;
  for (k = 0; k < 33; k++) q65[(k * k) % 65] = 1;
  t = zsmod(zn, 64l);
  r = zsmod(zn, 45045);
  if (q64[t] == 0) square = 0;
  if (square && q63[r % 63] == 0) square = 0;
  if (square && q65[r % 65] == 0) square = 0;
  if (square && q11[r % 11] == 0) square = 0;
  if (square) {
    zsqrt(zn, zq, &zd);
    if (zscompare(zd, 0l) != 0) square = 0;
  }
  zfree(&zd);
  return square;
}

int Cornacchia(verylong zd, verylong zp,
               verylong *zx, verylong *zy)
{
  int value;
  verylong zD = 0, za = 0, zb = 0, zc = 0, zl = 0;
  verylong zr = 0, zx0 = 0;

  zcopy(zd, &zD);
  znegate(&zD);
  if (zjacobi(zD, zp) != - 1) {
    zadd(zD, zp, &za);
    zsqrtmod(za, zp, &zx0);
    zrshift(zp, 1l, &za);
    zcopy(zx0, &zc);
    while (zcompare(zx0, za) <= 0) {
      zadd(zx0, zp, &zb);
      zcopy(zb, &zx0);
    }
    if (zcompare(zx0, zp) >= 0) {
      zcopy(zc, &zx0);
      znegate(&zx0);
      while (zcompare(zx0, za) <= 0) {
        zadd(zx0, zp, &zb);
        zcopy(zb, &zx0);
      }
    }
    zcopy(zp, &za);
    zcopy(zx0, &zb);
    zsqrt(zp, &zl, &zc);
    while (zcompare(zb, zl) > 0) {
      zmod(za, zb, &zr);
      zcopy(zb, &za);
      zcopy(zr, &zb);
    }
    zsq(zb, &za);
    zsub(zp, za, &zD);
    zdiv(zD, zd, &zc, &zr);
    if (zscompare(zr, 0l) == 0 && square_test(zc, &za)) {
      zcopy(zb, zx);
      zcopy(za, zy);
      value = 1;
    }
    else value = 0;
  }
  else value = 0;
  zfree(&zD);
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zl);
  zfree(&zr);
  zfree(&zx0);
  return value;
}

int main(void)
{
  char answer[256];
  int pr;
  verylong zd = 0, zp = 0, zx = 0, zy = 0;

  do {
    do {
      printf("p = ");
      zread(&zp);
      pr = zprobprime(zp, 5l);
      if (!pr) printf("*error*\np must be a prime\n");
    } while (!pr);
    do {
      printf("d = ");
      zread(&zd);
      pr = zscompare(zd, 0l) > 0 && zcompare(zd, zp) < 0;
      if (!pr) printf("*error*\nd must be 0 < d < p\n");
    } while (!pr);
    if (Cornacchia(zd, zp, &zx, &zy)) {
      printf("x = "); zwriteln(zx);
      printf("y = "); zwriteln(zy);
    }
    else printf("Diophantine equation has no solution\n");
    printf("another equation (n or y)? ");
    scanf("%s", answer);
  } while (tolower(answer[0]) == 'y');
  zfree(&zd);
  zfree(&zp);
  zfree(&zx);
  zfree(&zy);
  return 0;
}