/* Exercise mpz_*_kronecker_*() and mpz_jacobi() functions.

Copyright 1999-2004, 2013 Free Software Foundation, Inc.

This file is part of the GNU MP Library test suite.

The GNU MP Library test suite is free software; you can redistribute it

and/or modify it under the terms of the GNU General Public License as

published by the Free Software Foundation; either version 3 of the License,

or (at your option) any later version.

The GNU MP Library test suite is distributed in the hope that it will be

useful, but WITHOUT ANY WARRANTY; without even the implied warranty of

MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General

Public License for more details.

You should have received a copy of the GNU General Public License along with

the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */

/* With no arguments the various Kronecker/Jacobi symbol routines are

   checked against some test data and a lot of derived data.

   To check the test data against PARI-GP, run

	   t-jac -p | gp -q

   It takes a while because the output from "t-jac -p" is big.

   Enhancements:

   More big test cases than those given by check_squares_zi would be good.  */

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include "gmp.h"

#include "gmp-impl.h"

#include "tests.h"

#ifdef _LONG_LONG_LIMB

#define LL(l,ll)  ll

#else

#define LL(l,ll)  l

#endif

int option_pari = 0;

unsigned long

mpz_mod4 (mpz_srcptr z)

{

  mpz_t          m;

  unsigned long  ret;

  mpz_init (m);

  mpz_fdiv_r_2exp (m, z, 2);

  ret = mpz_get_ui (m);

  mpz_clear (m);

  return ret;

}

int

mpz_fits_ulimb_p (mpz_srcptr z)

{

  return (SIZ(z) == 1 || SIZ(z) == 0);

}

mp_limb_t

mpz_get_ulimb (mpz_srcptr z)

{

  if (SIZ(z) == 0)

    return 0;

  else

    return PTR(z)[0];

}

void

try_base (mp_limb_t a, mp_limb_t b, int answer)

{

  int  got;

  if ((b & 1) == 0 || b == 1 || a > b)

    return;

  got = mpn_jacobi_base (a, b, 0);

  if (got != answer)

    {

      printf (LL("mpn_jacobi_base (%lu, %lu) is %d should be %d\n",

		 "mpn_jacobi_base (%llu, %llu) is %d should be %d\n"),

	      a, b, got, answer);

      abort ();

    }

}

void

try_zi_ui (mpz_srcptr a, unsigned long b, int answer)

{

  int  got;

  got = mpz_kronecker_ui (a, b);

  if (got != answer)

    {

      printf ("mpz_kronecker_ui (");

      mpz_out_str (stdout, 10, a);

      printf (", %lu) is %d should be %d\n", b, got, answer);

      abort ();

    }

}

void

try_zi_si (mpz_srcptr a, long b, int answer)

{

  int  got;

  got = mpz_kronecker_si (a, b);

  if (got != answer)

    {

      printf ("mpz_kronecker_si (");

      mpz_out_str (stdout, 10, a);

      printf (", %ld) is %d should be %d\n", b, got, answer);

      abort ();

    }

}

void

try_ui_zi (unsigned long a, mpz_srcptr b, int answer)

{

  int  got;

  got = mpz_ui_kronecker (a, b);

  if (got != answer)

    {

      printf ("mpz_ui_kronecker (%lu, ", a);

      mpz_out_str (stdout, 10, b);

      printf (") is %d should be %d\n", got, answer);

      abort ();

    }

}

void

try_si_zi (long a, mpz_srcptr b, int answer)

{

  int  got;

  got = mpz_si_kronecker (a, b);

  if (got != answer)

    {

      printf ("mpz_si_kronecker (%ld, ", a);

      mpz_out_str (stdout, 10, b);

      printf (") is %d should be %d\n", got, answer);

      abort ();

    }

}

/* Don't bother checking mpz_jacobi, since it only differs for b even, and

   we don't have an actual expected answer for it.  tests/devel/try.c does

   some checks though.  */

void

try_zi_zi (mpz_srcptr a, mpz_srcptr b, int answer)

{

  int  got;

  got = mpz_kronecker (a, b);

  if (got != answer)

    {

      printf ("mpz_kronecker (");

      mpz_out_str (stdout, 10, a);

      printf (", ");

      mpz_out_str (stdout, 10, b);

      printf (") is %d should be %d\n", got, answer);

      abort ();

    }

}

void

try_pari (mpz_srcptr a, mpz_srcptr b, int answer)

{

  printf ("try(");

  mpz_out_str (stdout, 10, a);

  printf (",");

  mpz_out_str (stdout, 10, b);

  printf (",%d)\n", answer);

}

void

try_each (mpz_srcptr a, mpz_srcptr b, int answer)

{

#if 0

  fprintf(stderr, "asize = %d, bsize = %d\n",

	  mpz_sizeinbase (a, 2), mpz_sizeinbase (b, 2));

#endif

  if (option_pari)

    {

      try_pari (a, b, answer);

      return;

    }

  if (mpz_fits_ulimb_p (a) && mpz_fits_ulimb_p (b))

    try_base (mpz_get_ulimb (a), mpz_get_ulimb (b), answer);

  if (mpz_fits_ulong_p (b))

    try_zi_ui (a, mpz_get_ui (b), answer);

  if (mpz_fits_slong_p (b))

    try_zi_si (a, mpz_get_si (b), answer);

  if (mpz_fits_ulong_p (a))

    try_ui_zi (mpz_get_ui (a), b, answer);

  if (mpz_fits_sint_p (a))

    try_si_zi (mpz_get_si (a), b, answer);

  try_zi_zi (a, b, answer);

}

/* Try (a/b) and (a/-b). */

void

try_pn (mpz_srcptr a, mpz_srcptr b_orig, int answer)

{

  mpz_t  b;

  mpz_init_set (b, b_orig);

  try_each (a, b, answer);

  mpz_neg (b, b);

  if (mpz_sgn (a) < 0)

    answer = -answer;

  try_each (a, b, answer);

  mpz_clear (b);

}

/* Try (a+k*p/b) for various k, using the fact (a/b) is periodic in a with

   period p.  For b>0, p=b if b!=2mod4 or p=4*b if b==2mod4. */

void

try_periodic_num (mpz_srcptr a_orig, mpz_srcptr b, int answer)

{

  mpz_t  a, a_period;

  int    i;

  if (mpz_sgn (b) <= 0)

    return;

  mpz_init_set (a, a_orig);

  mpz_init_set (a_period, b);

  if (mpz_mod4 (b) == 2)

    mpz_mul_ui (a_period, a_period, 4);

  /* don't bother with these tests if they're only going to produce

     even/even */

  if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (a_period))

    goto done;

  for (i = 0; i < 6; i++)

    {

      mpz_add (a, a, a_period);

      try_pn (a, b, answer);

    }

  mpz_set (a, a_orig);

  for (i = 0; i < 6; i++)

    {

      mpz_sub (a, a, a_period);

      try_pn (a, b, answer);

    }

 done:

  mpz_clear (a);

  mpz_clear (a_period);

}

/* Try (a/b+k*p) for various k, using the fact (a/b) is periodic in b of

   period p.

			       period p

	   a==0,1mod4             a

	   a==2mod4              4*a

	   a==3mod4 and b odd    4*a

	   a==3mod4 and b even   8*a

   In Henri Cohen's book the period is given as 4*a for all a==2,3mod4, but

   a counterexample would seem to be (3/2)=-1 which with (3/14)=+1 doesn't

   have period 4*a (but rather 8*a with (3/26)=-1).  Maybe the plain 4*a is

   to be read as applying to a plain Jacobi symbol with b odd, rather than

   the Kronecker extension to b even. */

void

try_periodic_den (mpz_srcptr a, mpz_srcptr b_orig, int answer)

{

  mpz_t  b, b_period;

  int    i;

  if (mpz_sgn (a) == 0 || mpz_sgn (b_orig) == 0)

    return;

  mpz_init_set (b, b_orig);

  mpz_init_set (b_period, a);

  if (mpz_mod4 (a) == 3 && mpz_even_p (b))

    mpz_mul_ui (b_period, b_period, 8L);

  else if (mpz_mod4 (a) >= 2)

    mpz_mul_ui (b_period, b_period, 4L);

  /* don't bother with these tests if they're only going to produce

     even/even */

  if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (b_period))

    goto done;

  for (i = 0; i < 6; i++)

    {

      mpz_add (b, b, b_period);

      try_pn (a, b, answer);

    }

  mpz_set (b, b_orig);

  for (i = 0; i < 6; i++)

    {

      mpz_sub (b, b, b_period);

      try_pn (a, b, answer);

    }

 done:

  mpz_clear (b);

  mpz_clear (b_period);

}

static const unsigned long  ktable[] = {

  0, 1, 2, 3, 4, 5, 6, 7,

  GMP_NUMB_BITS-1, GMP_NUMB_BITS, GMP_NUMB_BITS+1,

  2*GMP_NUMB_BITS-1, 2*GMP_NUMB_BITS, 2*GMP_NUMB_BITS+1,

  3*GMP_NUMB_BITS-1, 3*GMP_NUMB_BITS, 3*GMP_NUMB_BITS+1

};

/* Try (a/b*2^k) for various k. */

void

try_2den (mpz_srcptr a, mpz_srcptr b_orig, int answer)

{

  mpz_t  b;

  int    kindex;

  int    answer_a2, answer_k;

  unsigned long k;

  /* don't bother when b==0 */

  if (mpz_sgn (b_orig) == 0)

    return;

  mpz_init_set (b, b_orig);

  /* (a/2) is 0 if a even, 1 if a==1 or 7 mod 8, -1 if a==3 or 5 mod 8 */

  answer_a2 = (mpz_even_p (a) ? 0

	       : (((SIZ(a) >= 0 ? PTR(a)[0] : -PTR(a)[0]) + 2) & 7) < 4 ? 1

	       : -1);

  for (kindex = 0; kindex < numberof (ktable); kindex++)

    {

      k = ktable[kindex];

      /* answer_k = answer*(answer_a2^k) */

      answer_k = (answer_a2 == 0 && k != 0 ? 0

		  : (k & 1) == 1 && answer_a2 == -1 ? -answer

		  : answer);

      mpz_mul_2exp (b, b_orig, k);

      try_pn (a, b, answer_k);

    }

  mpz_clear (b);

}

/* Try (a*2^k/b) for various k.  If it happens mpz_ui_kronecker() gets (2/b)

   wrong it will show up as wrong answers demanded. */

void

try_2num (mpz_srcptr a_orig, mpz_srcptr b, int answer)

{

  mpz_t  a;

  int    kindex;

  int    answer_2b, answer_k;

  unsigned long  k;

  /* don't bother when a==0 */

  if (mpz_sgn (a_orig) == 0)

    return;

  mpz_init (a);

  /* (2/b) is 0 if b even, 1 if b==1 or 7 mod 8, -1 if b==3 or 5 mod 8 */

  answer_2b = (mpz_even_p (b) ? 0

	       : (((SIZ(b) >= 0 ? PTR(b)[0] : -PTR(b)[0]) + 2) & 7) < 4 ? 1

	       : -1);

  for (kindex = 0; kindex < numberof (ktable); kindex++)

    {

      k = ktable[kindex];

      /* answer_k = answer*(answer_2b^k) */

      answer_k = (answer_2b == 0 && k != 0 ? 0

		  : (k & 1) == 1 && answer_2b == -1 ? -answer

		  : answer);

	mpz_mul_2exp (a, a_orig, k);

      try_pn (a, b, answer_k);

    }

  mpz_clear (a);

}

/* The try_2num() and try_2den() routines don't in turn call

   try_periodic_num() and try_periodic_den() because it hugely increases the

   number of tests performed, without obviously increasing coverage.

   Useful extra derived cases can be added here. */

void

try_all (mpz_t a, mpz_t b, int answer)

{

  try_pn (a, b, answer);

  try_periodic_num (a, b, answer);

  try_periodic_den (a, b, answer);

  try_2num (a, b, answer);

  try_2den (a, b, answer);

}

void

check_data (void)

{

  static const struct {

    const char  *a;

    const char  *b;

    int         answer;

  } data[] = {

    /* Note that the various derived checks in try_all() reduce the cases

       that need to be given here.  */

    /* some zeros */

    {  "0",  "0", 0 },

    {  "0",  "2", 0 },

    {  "0",  "6", 0 },

    {  "5",  "0", 0 },

    { "24", "60", 0 },

    /* (a/1) = 1, any a

       In particular note (0/1)=1 so that (a/b)=(a mod b/b). */

    { "0", "1", 1 },

    { "1", "1", 1 },

    { "2", "1", 1 },

    { "3", "1", 1 },

    { "4", "1", 1 },

    { "5", "1", 1 },

    /* (0/b) = 0, b != 1 */

    { "0",  "3", 0 },

    { "0",  "5", 0 },

    { "0",  "7", 0 },

    { "0",  "9", 0 },

    { "0", "11", 0 },

    { "0", "13", 0 },

    { "0", "15", 0 },

    /* (1/b) = 1 */

    { "1",  "1", 1 },

    { "1",  "3", 1 },

    { "1",  "5", 1 },

    { "1",  "7", 1 },

    { "1",  "9", 1 },

    { "1", "11", 1 },

    /* (-1/b) = (-1)^((b-1)/2) which is -1 for b==3 mod 4 */

    { "-1",  "1",  1 },

    { "-1",  "3", -1 },

    { "-1",  "5",  1 },

    { "-1",  "7", -1 },

    { "-1",  "9",  1 },

    { "-1", "11", -1 },

    { "-1", "13",  1 },

    { "-1", "15", -1 },

    { "-1", "17",  1 },

    { "-1", "19", -1 },

    /* (2/b) = (-1)^((b^2-1)/8) which is -1 for b==3,5 mod 8.

       try_2num() will exercise multiple powers of 2 in the numerator.  */

    { "2",  "1",  1 },

    { "2",  "3", -1 },

    { "2",  "5", -1 },

    { "2",  "7",  1 },

    { "2",  "9",  1 },

    { "2", "11", -1 },

    { "2", "13", -1 },

    { "2", "15",  1 },

    { "2", "17",  1 },

    /* (-2/b) = (-1)^((b^2-1)/8)*(-1)^((b-1)/2) which is -1 for b==5,7mod8.

       try_2num() will exercise multiple powers of 2 in the numerator, which

       will test that the shift in mpz_si_kronecker() uses unsigned not

       signed.  */

    { "-2",  "1",  1 },

    { "-2",  "3",  1 },

    { "-2",  "5", -1 },

    { "-2",  "7", -1 },

    { "-2",  "9",  1 },

    { "-2", "11",  1 },

    { "-2", "13", -1 },

    { "-2", "15", -1 },

    { "-2", "17",  1 },

    /* (a/2)=(2/a).

       try_2den() will exercise multiple powers of 2 in the denominator. */

    {  "3",  "2", -1 },

    {  "5",  "2", -1 },

    {  "7",  "2",  1 },

    {  "9",  "2",  1 },

    {  "11", "2", -1 },

    /* Harriet Griffin, "Elementary Theory of Numbers", page 155, various

       examples.  */

    {   "2", "135",  1 },

    { "135",  "19", -1 },

    {   "2",  "19", -1 },

    {  "19", "135",  1 },

    { "173", "135",  1 },

    {  "38", "135",  1 },

    { "135", "173",  1 },

    { "173",   "5", -1 },

    {   "3",   "5", -1 },

    {   "5", "173", -1 },

    { "173",   "3", -1 },

    {   "2",   "3", -1 },

    {   "3", "173", -1 },

    { "253",  "21",  1 },

    {   "1",  "21",  1 },

    {  "21", "253",  1 },

    {  "21",  "11", -1 },

    {  "-1",  "11", -1 },

    /* Griffin page 147 */

    {  "-1",  "17",  1 },

    {   "2",  "17",  1 },

    {  "-2",  "17",  1 },

    {  "-1",  "89",  1 },

    {   "2",  "89",  1 },

    /* Griffin page 148 */

    {  "89",  "11",  1 },

    {   "1",  "11",  1 },

    {  "89",   "3", -1 },

    {   "2",   "3", -1 },

    {   "3",  "89", -1 },

    {  "11",  "89",  1 },

    {  "33",  "89", -1 },

    /* H. Davenport, "The Higher Arithmetic", page 65, the quadratic

       residues and non-residues mod 19.  */

    {  "1", "19",  1 },

    {  "4", "19",  1 },

    {  "5", "19",  1 },

    {  "6", "19",  1 },

    {  "7", "19",  1 },

    {  "9", "19",  1 },

    { "11", "19",  1 },

    { "16", "19",  1 },

    { "17", "19",  1 },

    {  "2", "19", -1 },

    {  "3", "19", -1 },

    {  "8", "19", -1 },

    { "10", "19", -1 },

    { "12", "19", -1 },

    { "13", "19", -1 },

    { "14", "19", -1 },

    { "15", "19", -1 },

    { "18", "19", -1 },

    /* Residues and non-residues mod 13 */

    {  "0",  "13",  0 },

    {  "1",  "13",  1 },

    {  "2",  "13", -1 },

    {  "3",  "13",  1 },

    {  "4",  "13",  1 },

    {  "5",  "13", -1 },

    {  "6",  "13", -1 },

    {  "7",  "13", -1 },

    {  "8",  "13", -1 },

    {  "9",  "13",  1 },

    { "10",  "13",  1 },

    { "11",  "13", -1 },

    { "12",  "13",  1 },

    /* various */

    {  "5",   "7", -1 },

    { "15",  "17",  1 },

    { "67",  "89",  1 },

    /* special values inducing a==b==1 at the end of jac_or_kron() */

    { "0x10000000000000000000000000000000000000000000000001",

      "0x10000000000000000000000000000000000000000000000003", 1 },

    /* Test for previous bugs in jacobi_2. */

    { "0x43900000000", "0x42400000439", -1 }, /* 32-bit limbs */

    { "0x4390000000000000000", "0x4240000000000000439", -1 }, /* 64-bit limbs */

    { "198158408161039063", "198158360916398807", -1 },

    /* Some tests involving large quotients in the continued fraction

       expansion. */

    { "37200210845139167613356125645445281805",

      "451716845976689892447895811408978421929", -1 },

    { "67674091930576781943923596701346271058970643542491743605048620644676477275152701774960868941561652032482173612421015",

      "4902678867794567120224500687210807069172039735", 0 },

    { "2666617146103764067061017961903284334497474492754652499788571378062969111250584288683585223600172138551198546085281683283672592", "2666617146103764067061017961903284334497474492754652499788571378062969111250584288683585223600172138551198546085281683290481773", 1 },

    /* Exercises the case asize == 1, btwos > 0 in mpz_jacobi. */

    { "804609", "421248363205206617296534688032638102314410556521742428832362659824", 1 } ,

    { "4190209", "2239744742177804210557442048984321017460028974602978995388383905961079286530650825925074203175536427000", 1 },

    /* Exercises the case asize == 1, btwos = 63 in mpz_jacobi

       (relevant when GMP_LIMB_BITS == 64). */

    { "17311973299000934401", "1675975991242824637446753124775689449936871337036614677577044717424700351103148799107651171694863695242089956242888229458836426332300124417011114380886016", 1 },

    { "3220569220116583677", "41859917623035396746", -1 },

    /* Other test cases that triggered bugs during development. */

    { "37200210845139167613356125645445281805", "340116213441272389607827434472642576514", -1 },

    { "74400421690278335226712251290890563610", "451716845976689892447895811408978421929", -1 },

  };

  int    i;

  mpz_t  a, b;

  mpz_init (a);

  mpz_init (b);

  for (i = 0; i < numberof (data); i++)

    {

      mpz_set_str_or_abort (a, data[i].a, 0);

      mpz_set_str_or_abort (b, data[i].b, 0);

      try_all (a, b, data[i].answer);

    }

  mpz_clear (a);

  mpz_clear (b);

}

/* (a^2/b)=1 if gcd(a,b)=1, or (a^2/b)=0 if gcd(a,b)!=1.

   This includes when a=0 or b=0. */

void

check_squares_zi (void)

{

  gmp_randstate_ptr rands = RANDS;

  mpz_t  a, b, g;

  int    i, answer;

  mp_size_t size_range, an, bn;

  mpz_t bs;

  mpz_init (bs);

  mpz_init (a);

  mpz_init (b);

  mpz_init (g);

  for (i = 0; i < 50; i++)

    {

      mpz_urandomb (bs, rands, 32);

      size_range = mpz_get_ui (bs) % 10 + i/8 + 2;

      mpz_urandomb (bs, rands, size_range);

      an = mpz_get_ui (bs);

      mpz_rrandomb (a, rands, an);

      mpz_urandomb (bs, rands, size_range);

      bn = mpz_get_ui (bs);

      mpz_rrandomb (b, rands, bn);

      mpz_gcd (g, a, b);

      if (mpz_cmp_ui (g, 1L) == 0)

	answer = 1;

      else

	answer = 0;

      mpz_mul (a, a, a);

      try_all (a, b, answer);

    }

  mpz_clear (bs);

  mpz_clear (a);

  mpz_clear (b);

  mpz_clear (g);

}

/* Check the handling of asize==0, make sure it isn't affected by the low

   limb. */

void

check_a_zero (void)

{

  mpz_t  a, b;

  mpz_init_set_ui (a, 0);

  mpz_init (b);

  mpz_set_ui (b, 1L);

  PTR(a)[0] = 0;

  try_all (a, b, 1);   /* (0/1)=1 */

  PTR(a)[0] = 1;

  try_all (a, b, 1);   /* (0/1)=1 */

  mpz_set_si (b, -1L);

  PTR(a)[0] = 0;

  try_all (a, b, 1);   /* (0/-1)=1 */

  PTR(a)[0] = 1;

  try_all (a, b, 1);   /* (0/-1)=1 */

  mpz_set_ui (b, 0);

  PTR(a)[0] = 0;

  try_all (a, b, 0);   /* (0/0)=0 */

  PTR(a)[0] = 1;

  try_all (a, b, 0);   /* (0/0)=0 */

  mpz_set_ui (b, 2);

  PTR(a)[0] = 0;

  try_all (a, b, 0);   /* (0/2)=0 */

  PTR(a)[0] = 1;

  try_all (a, b, 0);   /* (0/2)=0 */

  mpz_clear (a);

  mpz_clear (b);

}

/* Assumes that b = prod p_k^e_k */

int

ref_jacobi (mpz_srcptr a, mpz_srcptr b, unsigned nprime,

	    mpz_t prime[], unsigned *exp)

{

  unsigned i;

  int res;

  for (i = 0, res = 1; i < nprime; i++)

    if (exp[i])

      {

	int legendre = refmpz_legendre (a, prime[i]);

	if (!legendre)

	  return 0;

	if (exp[i] & 1)

	  res *= legendre;

      }

  return res;

}

void

check_jacobi_factored (void)

{

#define PRIME_N 10

#define PRIME_MAX_SIZE 50

#define PRIME_MAX_EXP 4

#define PRIME_A_COUNT 10

#define PRIME_B_COUNT 5

#define PRIME_MAX_B_SIZE 2000

  gmp_randstate_ptr rands = RANDS;

  mpz_t prime[PRIME_N];

  unsigned exp[PRIME_N];

  mpz_t a, b, t, bs;

  unsigned i;

  mpz_init (a);

  mpz_init (b);

  mpz_init (t);

  mpz_init (bs);

  /* Generate primes */

  for (i = 0; i < PRIME_N; i++)

    {

      mp_size_t size;

      mpz_init (prime[i]);

      mpz_urandomb (bs, rands, 32);

      size = mpz_get_ui (bs) % PRIME_MAX_SIZE + 2;

      mpz_rrandomb (prime[i], rands, size);

      if (mpz_cmp_ui (prime[i], 3) <= 0)

	mpz_set_ui (prime[i], 3);

      else

	mpz_nextprime (prime[i], prime[i]);

    }

  for (i = 0; i < PRIME_B_COUNT; i++)

    {

      unsigned j, k;

      mp_bitcnt_t bsize;

      mpz_set_ui (b, 1);

      bsize = 1;

      for (j = 0; j < PRIME_N && bsize < PRIME_MAX_B_SIZE; j++)

	{

	  mpz_urandomb (bs, rands, 32);

	  exp[j] = mpz_get_ui (bs) % PRIME_MAX_EXP;

	  mpz_pow_ui (t, prime[j], exp[j]);

	  mpz_mul (b, b, t);

	  bsize = mpz_sizeinbase (b, 2);

	}

      for (k = 0; k < PRIME_A_COUNT; k++)

	{

	  int answer;

	  mpz_rrandomb (a, rands, bsize + 2);

	  answer = ref_jacobi (a, b, j, prime, exp);

	  try_all (a, b, answer);

	}

    }

  for (i = 0; i < PRIME_N; i++)

    mpz_clear (prime[i]);

  mpz_clear (a);

  mpz_clear (b);

  mpz_clear (t);

  mpz_clear (bs);

#undef PRIME_N

#undef PRIME_MAX_SIZE

#undef PRIME_MAX_EXP

#undef PRIME_A_COUNT

#undef PRIME_B_COUNT

#undef PRIME_MAX_B_SIZE

}

/* These tests compute (a|n), where the quotient sequence includes

   large quotients, and n has a known factorization. Such inputs are

   generated as follows. First, construct a large n, as a power of a

   prime p of moderate size.

   Next, compute a matrix from factors (q,1;1,0), with q chosen with

   uniformly distributed size. We must stop with matrix elements of

   roughly half the size of n. Denote elements of M as M = (m00, m01;

   m10, m11).

   We now look for solutions to

     n = m00 x + m01 y

     a = m10 x + m11 y

   with x,y > 0. Since n >= m00 * m01, there exists a positive

   solution to the first equation. Find those x, y, and substitute in

   the second equation to get a. Then the quotient sequence for (a|n)

   is precisely the quotients used when constructing M, followed by

   the quotient sequence for (x|y).

   Numbers should also be large enough that we exercise hgcd_jacobi,

   which means that they should be larger than

     max (GCD_DC_THRESHOLD, 3 * HGCD_THRESHOLD)

   With an n of roughly 40000 bits, this should hold on most machines.

*/

void

check_large_quotients (void)

{

#define COUNT 50

#define PBITS 200