third_party/gmp/tests/refmpz.c - RealtimeRoboticsGroup/test - Gitiles

 /* Reference mpz functions.

 Copyright 1997, 1999-2002 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/.  */

 /* always do assertion checking */
 #define WANT_ASSERT  1

 #include <stdio.h>
 #include <stdlib.h> /* for free */
 #include "gmp-impl.h"
 #include "longlong.h"
 #include "tests.h"


 /* Change this to "#define TRACE(x) x" for some traces. */
 #define TRACE(x)


 /* FIXME: Shouldn't use plain mpz functions in a reference routine. */
 void
 refmpz_combit (mpz_ptr r, unsigned long bit)
 {
   if (mpz_tstbit (r, bit))
     mpz_clrbit (r, bit);
   else
     mpz_setbit (r, bit);
 }


 unsigned long
 refmpz_hamdist (mpz_srcptr x, mpz_srcptr y)
 {
   mp_size_t      xsize, ysize, tsize;
   mp_ptr         xp, yp;
   unsigned long  ret;

   if ((SIZ(x) < 0 && SIZ(y) >= 0)
       || (SIZ(y) < 0 && SIZ(x) >= 0))
     return ULONG_MAX;

   xsize = ABSIZ(x);
   ysize = ABSIZ(y);
   tsize = MAX (xsize, ysize);

   xp = refmpn_malloc_limbs (tsize);
   refmpn_zero (xp, tsize);
   refmpn_copy (xp, PTR(x), xsize);

   yp = refmpn_malloc_limbs (tsize);
   refmpn_zero (yp, tsize);
   refmpn_copy (yp, PTR(y), ysize);

   if (SIZ(x) < 0)
     refmpn_neg (xp, xp, tsize);

   if (SIZ(x) < 0)
     refmpn_neg (yp, yp, tsize);

   ret = refmpn_hamdist (xp, yp, tsize);

   free (xp);
   free (yp);
   return ret;
 }

 void
 refmpz_gcd (mpz_ptr g, mpz_srcptr a_orig, mpz_srcptr b_orig)
 {
   mp_bitcnt_t a_twos, b_twos, common_twos;
   mpz_t a;
   mpz_t b;
   mpz_init (a);
   mpz_init (b);
   mpz_abs (a, a_orig);
   mpz_abs (b, b_orig);

   if (mpz_sgn (a) == 0)
     {
       mpz_set (g, b);
       return;
     }
   if (mpz_sgn (b) == 0)
     {
       mpz_set (g, a);
       return;
     }
   a_twos = mpz_scan1 (a, 0);
   mpz_tdiv_q_2exp (a, a, a_twos);

   b_twos = mpz_scan1 (b, 0);
   mpz_tdiv_q_2exp (b, b, b_twos);

   common_twos = MIN(a_twos, b_twos);
   for (;;)
     {
       int c;
       mp_bitcnt_t twos;
       c = mpz_cmp (a, b);
       if (c == 0)
 	break;
       if (c < 0)
 	mpz_swap (a, b);
       mpz_sub (a, a, b);
       twos = mpz_scan1 (a, 0);
       mpz_tdiv_q_2exp (a, a, twos);
     }
   mpz_mul_2exp (g, a, common_twos);

   mpz_clear (a);
   mpz_clear (b);
 }

 /* (0/b), with mpz b; is 1 if b=+/-1, 0 otherwise */
 #define JACOBI_0Z(b)  JACOBI_0LS (PTR(b)[0], SIZ(b))

 /* (a/b) effect due to sign of b: mpz/mpz */
 #define JACOBI_BSGN_ZZ_BIT1(a, b)   JACOBI_BSGN_SS_BIT1 (SIZ(a), SIZ(b))

 /* (a/b) effect due to sign of a: mpz/unsigned-mpz, b odd;
    is (-1/b) if a<0, or +1 if a>=0 */
 #define JACOBI_ASGN_ZZU_BIT1(a, b)  JACOBI_ASGN_SU_BIT1 (SIZ(a), PTR(b)[0])

 int
 refmpz_kronecker (mpz_srcptr a_orig, mpz_srcptr b_orig)
 {
   unsigned long  twos;
   mpz_t  a, b;
   int    result_bit1 = 0;

   if (mpz_sgn (b_orig) == 0)
     return JACOBI_Z0 (a_orig);  /* (a/0) */

   if (mpz_sgn (a_orig) == 0)
     return JACOBI_0Z (b_orig);  /* (0/b) */

   if (mpz_even_p (a_orig) && mpz_even_p (b_orig))
     return 0;

   if (mpz_cmp_ui (b_orig, 1) == 0)
     return 1;

   mpz_init_set (a, a_orig);
   mpz_init_set (b, b_orig);

   if (mpz_sgn (b) < 0)
     {
       result_bit1 ^= JACOBI_BSGN_ZZ_BIT1 (a, b);
       mpz_neg (b, b);
     }
   if (mpz_even_p (b))
     {
       twos = mpz_scan1 (b, 0L);
       mpz_tdiv_q_2exp (b, b, twos);
       result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]);
     }

   if (mpz_sgn (a) < 0)
     {
       result_bit1 ^= JACOBI_N1B_BIT1 (PTR(b)[0]);
       mpz_neg (a, a);
     }
   if (mpz_even_p (a))
     {
       twos = mpz_scan1 (a, 0L);
       mpz_tdiv_q_2exp (a, a, twos);
       result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]);
     }

   for (;;)
     {
       ASSERT (mpz_odd_p (a));
       ASSERT (mpz_odd_p (b));
       ASSERT (mpz_sgn (a) > 0);
       ASSERT (mpz_sgn (b) > 0);

       TRACE (printf ("top\n");
 	     mpz_trace (" a", a);
 	     mpz_trace (" b", b));

       if (mpz_cmp (a, b) < 0)
 	{
 	  TRACE (printf ("swap\n"));
 	  mpz_swap (a, b);
 	  result_bit1 ^= JACOBI_RECIP_UU_BIT1 (PTR(a)[0], PTR(b)[0]);
 	}

       if (mpz_cmp_ui (b, 1) == 0)
 	break;

       mpz_sub (a, a, b);
       TRACE (printf ("sub\n");
 	     mpz_trace (" a", a));
       if (mpz_sgn (a) == 0)
 	goto zero;

       twos = mpz_scan1 (a, 0L);
       mpz_fdiv_q_2exp (a, a, twos);
       TRACE (printf ("twos %lu\n", twos);
 	     mpz_trace (" a", a));
       result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]);
     }

   mpz_clear (a);
   mpz_clear (b);
   return JACOBI_BIT1_TO_PN (result_bit1);

  zero:
   mpz_clear (a);
   mpz_clear (b);
   return 0;
 }

 /* Same as mpz_kronecker, but ignoring factors of 2 on b */
 int
 refmpz_jacobi (mpz_srcptr a, mpz_srcptr b)
 {
   ASSERT_ALWAYS (mpz_sgn (b) > 0);
   ASSERT_ALWAYS (mpz_odd_p (b));

   return refmpz_kronecker (a, b);
 }

 /* Legendre symbol via powm. p must be an odd prime. */
 int
 refmpz_legendre (mpz_srcptr a, mpz_srcptr p)
 {
   int res;

   mpz_t r;
   mpz_t e;

   ASSERT_ALWAYS (mpz_sgn (p) > 0);
   ASSERT_ALWAYS (mpz_odd_p (p));

   mpz_init (r);
   mpz_init (e);

   mpz_fdiv_r (r, a, p);

   mpz_set (e, p);
   mpz_sub_ui (e, e, 1);
   mpz_fdiv_q_2exp (e, e, 1);
   mpz_powm (r, r, e, p);

   /* Normalize to a more or less symmetric range around zero */
   if (mpz_cmp (r, e) > 0)
     mpz_sub (r, r, p);

   ASSERT_ALWAYS (mpz_cmpabs_ui (r, 1) <= 0);

   res = mpz_sgn (r);

   mpz_clear (r);
   mpz_clear (e);

   return res;
 }


 int
 refmpz_kronecker_ui (mpz_srcptr a, unsigned long b)
 {
   mpz_t  bz;
   int    ret;
   mpz_init_set_ui (bz, b);
   ret = refmpz_kronecker (a, bz);
   mpz_clear (bz);
   return ret;
 }

 int
 refmpz_kronecker_si (mpz_srcptr a, long b)
 {
   mpz_t  bz;
   int    ret;
   mpz_init_set_si (bz, b);
   ret = refmpz_kronecker (a, bz);
   mpz_clear (bz);
   return ret;
 }

 int
 refmpz_ui_kronecker (unsigned long a, mpz_srcptr b)
 {
   mpz_t  az;
   int    ret;
   mpz_init_set_ui (az, a);
   ret = refmpz_kronecker (az, b);
   mpz_clear (az);
   return ret;
 }

 int
 refmpz_si_kronecker (long a, mpz_srcptr b)
 {
   mpz_t  az;
   int    ret;
   mpz_init_set_si (az, a);
   ret = refmpz_kronecker (az, b);
   mpz_clear (az);
   return ret;
 }


 void
 refmpz_pow_ui (mpz_ptr w, mpz_srcptr b, unsigned long e)
 {
   mpz_t          s, t;
   unsigned long  i;

   mpz_init_set_ui (t, 1L);
   mpz_init_set (s, b);

   if ((e & 1) != 0)
     mpz_mul (t, t, s);

   for (i = 2; i <= e; i <<= 1)
     {
       mpz_mul (s, s, s);
       if ((i & e) != 0)
 	mpz_mul (t, t, s);
     }

   mpz_set (w, t);

   mpz_clear (s);
   mpz_clear (t);
 }
	/* Reference mpz functions.

	Copyright 1997, 1999-2002 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/. */

	/* always do assertion checking */
	#define WANT_ASSERT 1

	#include <stdio.h>
	#include <stdlib.h> /* for free */
	#include "gmp-impl.h"
	#include "longlong.h"
	#include "tests.h"


	/* Change this to "#define TRACE(x) x" for some traces. */
	#define TRACE(x)


	/* FIXME: Shouldn't use plain mpz functions in a reference routine. */
	void
	refmpz_combit (mpz_ptr r, unsigned long bit)
	{
	if (mpz_tstbit (r, bit))
	mpz_clrbit (r, bit);
	else
	mpz_setbit (r, bit);
	}


	unsigned long
	refmpz_hamdist (mpz_srcptr x, mpz_srcptr y)
	{
	mp_size_t xsize, ysize, tsize;
	mp_ptr xp, yp;
	unsigned long ret;

	if ((SIZ(x) < 0 && SIZ(y) >= 0)
	\|\| (SIZ(y) < 0 && SIZ(x) >= 0))
	return ULONG_MAX;

	xsize = ABSIZ(x);
	ysize = ABSIZ(y);
	tsize = MAX (xsize, ysize);

	xp = refmpn_malloc_limbs (tsize);
	refmpn_zero (xp, tsize);
	refmpn_copy (xp, PTR(x), xsize);

	yp = refmpn_malloc_limbs (tsize);
	refmpn_zero (yp, tsize);
	refmpn_copy (yp, PTR(y), ysize);

	if (SIZ(x) < 0)
	refmpn_neg (xp, xp, tsize);

	if (SIZ(x) < 0)
	refmpn_neg (yp, yp, tsize);

	ret = refmpn_hamdist (xp, yp, tsize);

	free (xp);
	free (yp);
	return ret;
	}

	void
	refmpz_gcd (mpz_ptr g, mpz_srcptr a_orig, mpz_srcptr b_orig)
	{
	mp_bitcnt_t a_twos, b_twos, common_twos;
	mpz_t a;
	mpz_t b;
	mpz_init (a);
	mpz_init (b);
	mpz_abs (a, a_orig);
	mpz_abs (b, b_orig);

	if (mpz_sgn (a) == 0)
	{
	mpz_set (g, b);
	return;
	}
	if (mpz_sgn (b) == 0)
	{
	mpz_set (g, a);
	return;
	}
	a_twos = mpz_scan1 (a, 0);
	mpz_tdiv_q_2exp (a, a, a_twos);

	b_twos = mpz_scan1 (b, 0);
	mpz_tdiv_q_2exp (b, b, b_twos);

	common_twos = MIN(a_twos, b_twos);
	for (;;)
	{
	int c;
	mp_bitcnt_t twos;
	c = mpz_cmp (a, b);
	if (c == 0)
	break;
	if (c < 0)
	mpz_swap (a, b);
	mpz_sub (a, a, b);
	twos = mpz_scan1 (a, 0);
	mpz_tdiv_q_2exp (a, a, twos);
	}
	mpz_mul_2exp (g, a, common_twos);

	mpz_clear (a);
	mpz_clear (b);
	}

	/* (0/b), with mpz b; is 1 if b=+/-1, 0 otherwise */
	#define JACOBI_0Z(b) JACOBI_0LS (PTR(b)[0], SIZ(b))

	/* (a/b) effect due to sign of b: mpz/mpz */
	#define JACOBI_BSGN_ZZ_BIT1(a, b) JACOBI_BSGN_SS_BIT1 (SIZ(a), SIZ(b))

	/* (a/b) effect due to sign of a: mpz/unsigned-mpz, b odd;
	is (-1/b) if a<0, or +1 if a>=0 */
	#define JACOBI_ASGN_ZZU_BIT1(a, b) JACOBI_ASGN_SU_BIT1 (SIZ(a), PTR(b)[0])

	int
	refmpz_kronecker (mpz_srcptr a_orig, mpz_srcptr b_orig)
	{
	unsigned long twos;
	mpz_t a, b;
	int result_bit1 = 0;

	if (mpz_sgn (b_orig) == 0)
	return JACOBI_Z0 (a_orig); /* (a/0) */

	if (mpz_sgn (a_orig) == 0)
	return JACOBI_0Z (b_orig); /* (0/b) */

	if (mpz_even_p (a_orig) && mpz_even_p (b_orig))
	return 0;

	if (mpz_cmp_ui (b_orig, 1) == 0)
	return 1;

	mpz_init_set (a, a_orig);
	mpz_init_set (b, b_orig);

	if (mpz_sgn (b) < 0)
	{
	result_bit1 ^= JACOBI_BSGN_ZZ_BIT1 (a, b);
	mpz_neg (b, b);
	}
	if (mpz_even_p (b))
	{
	twos = mpz_scan1 (b, 0L);
	mpz_tdiv_q_2exp (b, b, twos);
	result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]);
	}

	if (mpz_sgn (a) < 0)
	{
	result_bit1 ^= JACOBI_N1B_BIT1 (PTR(b)[0]);
	mpz_neg (a, a);
	}
	if (mpz_even_p (a))
	{
	twos = mpz_scan1 (a, 0L);
	mpz_tdiv_q_2exp (a, a, twos);
	result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]);
	}

	for (;;)
	{
	ASSERT (mpz_odd_p (a));
	ASSERT (mpz_odd_p (b));
	ASSERT (mpz_sgn (a) > 0);
	ASSERT (mpz_sgn (b) > 0);

	TRACE (printf ("top\n");
	mpz_trace (" a", a);
	mpz_trace (" b", b));

	if (mpz_cmp (a, b) < 0)
	{
	TRACE (printf ("swap\n"));
	mpz_swap (a, b);
	result_bit1 ^= JACOBI_RECIP_UU_BIT1 (PTR(a)[0], PTR(b)[0]);
	}

	if (mpz_cmp_ui (b, 1) == 0)
	break;

	mpz_sub (a, a, b);
	TRACE (printf ("sub\n");
	mpz_trace (" a", a));
	if (mpz_sgn (a) == 0)
	goto zero;

	twos = mpz_scan1 (a, 0L);
	mpz_fdiv_q_2exp (a, a, twos);
	TRACE (printf ("twos %lu\n", twos);
	mpz_trace (" a", a));
	result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]);
	}

	mpz_clear (a);
	mpz_clear (b);
	return JACOBI_BIT1_TO_PN (result_bit1);

	zero:
	mpz_clear (a);
	mpz_clear (b);
	return 0;
	}

	/* Same as mpz_kronecker, but ignoring factors of 2 on b */
	int
	refmpz_jacobi (mpz_srcptr a, mpz_srcptr b)
	{
	ASSERT_ALWAYS (mpz_sgn (b) > 0);
	ASSERT_ALWAYS (mpz_odd_p (b));

	return refmpz_kronecker (a, b);
	}

	/* Legendre symbol via powm. p must be an odd prime. */
	int
	refmpz_legendre (mpz_srcptr a, mpz_srcptr p)
	{
	int res;

	mpz_t r;
	mpz_t e;

	ASSERT_ALWAYS (mpz_sgn (p) > 0);
	ASSERT_ALWAYS (mpz_odd_p (p));

	mpz_init (r);
	mpz_init (e);

	mpz_fdiv_r (r, a, p);

	mpz_set (e, p);
	mpz_sub_ui (e, e, 1);
	mpz_fdiv_q_2exp (e, e, 1);
	mpz_powm (r, r, e, p);

	/* Normalize to a more or less symmetric range around zero */
	if (mpz_cmp (r, e) > 0)
	mpz_sub (r, r, p);

	ASSERT_ALWAYS (mpz_cmpabs_ui (r, 1) <= 0);

	res = mpz_sgn (r);

	mpz_clear (r);
	mpz_clear (e);

	return res;
	}


	int
	refmpz_kronecker_ui (mpz_srcptr a, unsigned long b)
	{
	mpz_t bz;
	int ret;
	mpz_init_set_ui (bz, b);
	ret = refmpz_kronecker (a, bz);
	mpz_clear (bz);
	return ret;
	}

	int
	refmpz_kronecker_si (mpz_srcptr a, long b)
	{
	mpz_t bz;
	int ret;
	mpz_init_set_si (bz, b);
	ret = refmpz_kronecker (a, bz);
	mpz_clear (bz);
	return ret;
	}

	int
	refmpz_ui_kronecker (unsigned long a, mpz_srcptr b)
	{
	mpz_t az;
	int ret;
	mpz_init_set_ui (az, a);
	ret = refmpz_kronecker (az, b);
	mpz_clear (az);
	return ret;
	}

	int
	refmpz_si_kronecker (long a, mpz_srcptr b)
	{
	mpz_t az;
	int ret;
	mpz_init_set_si (az, a);
	ret = refmpz_kronecker (az, b);
	mpz_clear (az);
	return ret;
	}


	void
	refmpz_pow_ui (mpz_ptr w, mpz_srcptr b, unsigned long e)
	{
	mpz_t s, t;
	unsigned long i;

	mpz_init_set_ui (t, 1L);
	mpz_init_set (s, b);

	if ((e & 1) != 0)
	mpz_mul (t, t, s);

	for (i = 2; i <= e; i <<= 1)
	{
	mpz_mul (s, s, s);
	if ((i & e) != 0)
	mpz_mul (t, t, s);
	}

	mpz_set (w, t);

	mpz_clear (s);
	mpz_clear (t);
	}