/*
  Author:  Pate Williams (c) 1997

  Algorithm 1.5.3 (Modified Cornacchia). See "A Course
  in Computational Algebraic Number Theory" by Henri
  Cohen page 36. Let p be a prime number and D be a
  negative number such that D = 0 or 1 modulo 4 and
  | D | < 4 * p. This algorithm either outputs an
  integer solution (x, y) to the Diophantine equation
  x * x + | D | * y * y = 4 * 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 modified_Cornacchia(verylong zD, verylong zp,
                        verylong *zx, verylong *zy)
{
  int value;
  long d, x;
  verylong za = 0, zb = 0, zc = 0, zd = 0, ze = 0;
  verylong zl = 0, zr = 0, zx0 = 0;

  if (zscompare(zp, 2l) == 0) {
    zsadd(zD, 8l, &za);
    if (zscompare(za, 0l) > 0) {
      if (square_test(za, zx)) {
        zone(zy);
        value = 1;
      }
      else value = 0;
    }
    else value = 0;
  }
  else {
    if (zjacobi(zD, zp) == - 1) value = 0;
    else {
      zsqrtmod(zD, zp, &zx0);
      if (zscompare(zx0, 0l) < 0) {
        zadd(zx0, zp, &za);
        zcopy(za, &zx0);
      }
      d = zsmod(zD, 2l);
      x = zsmod(zx0, 2l);
      if (d != x) {
        zsub(zp, zx0, &za);
        zcopy(za, &zx0);
      }
      zcopy(zp, &za);
      zcopy(zx0, &zb);
      zsqrt(zp, &zc, &zd);
      zlshift(zc, 1l, &zl);
      while (zcompare(zb, zl) > 0) {
        zmod(za, zb, &zr);
        zcopy(zb, &za);
        zcopy(zr, &zb);
      }
      zlshift(zp, 2l, &zc);
      zcopy(zD, &zd);
      zabs(&zd);
      zsq(zb, &za);
      zsub(zc, za, &ze);
      zdiv(ze, zd, &zc, &zr);
      if (zscompare(zr, 0l) == 0 && square_test(zc, zy)) {
        zcopy(zb, zx);
        value = 1;
      }
      else value = 0;
    }
  }
  zfree(&za);
  zfree(&zb);
  zfree(&zc);
  zfree(&zd);
  zfree(&ze);
  zfree(&zl);
  zfree(&zr);
  zfree(&zx0);
  return value;
}

int main(void)
{
  char answer[256];
  int pr;
  long m;
  verylong zD = 0, zP = 0, 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);
    zlshift(zp, 2l, &zP);
    do {
      printf("D = ");
      zread(&zD);
      if (zscompare(zD, 0l) < 0) {
        zcopy(zD, &zd);
        zabs(&zd);
        pr = zcompare(zd, zP) < 0;
        if (pr) {
          m = zsmod(zD, 4l);
          pr = m == 0 || m == 1;
          if (!pr) printf("*error*\nmust have D = 0 or 1 modulo 4\n");
        }
        else
          printf("*error*\nD must be | D | < 4 * p\n");
      }
      else {
        printf("D must be negative\n");
        pr = 0;
      }
    } while (!pr);
    if (modified_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;
}