mirror of
https://gitlab.com/embeddable-common-lisp/ecl.git
synced 2026-01-16 06:12:25 -08:00
Removed Kronecker and Jacobi functions from GMP/mpz
This commit is contained in:
Juanjo Garcia-Ripoll
parent
08b5c0345b
commit
a73e7842ae
7 changed files with 19 additions and 786 deletions
|
|
@ -84,7 +84,7 @@ LIBMP_LT_REVISION = 10
|
|||
LIBMP_LT_AGE = 1
|
||||
|
||||
|
||||
SUBDIRS = mpn mpz # mpq cxx mpf printf scanf tests demos tune
|
||||
SUBDIRS = mpn mpz # mpq cxx mpf tests demos tune
|
||||
|
||||
EXTRA_DIST = macos configfsf.guess configfsf.sub .gdbinit INSTALL.autoconf
|
||||
|
||||
|
|
@ -160,8 +160,7 @@ MPZ_OBJECTS = mpz/abs$U.lo mpz/add$U.lo mpz/add_ui$U.lo \
|
|||
mpz/import$U.lo mpz/init$U.lo mpz/init2$U.lo mpz/inp_raw$U.lo \
|
||||
mpz/inp_str$U.lo mpz/invert$U.lo \
|
||||
mpz/ior$U.lo mpz/iset$U.lo mpz/iset_d$U.lo mpz/iset_si$U.lo \
|
||||
mpz/iset_str$U.lo mpz/iset_ui$U.lo mpz/jacobi$U.lo mpz/kronsz$U.lo \
|
||||
mpz/kronuz$U.lo mpz/kronzs$U.lo mpz/kronzu$U.lo \
|
||||
mpz/iset_str$U.lo mpz/iset_ui$U.lo \
|
||||
mpz/lcm$U.lo mpz/lcm_ui$U.lo mpz/lucnum_ui$U.lo mpz/lucnum2_ui$U.lo \
|
||||
mpz/millerrabin$U.lo mpz/mod$U.lo mpz/mul$U.lo mpz/mul_2exp$U.lo \
|
||||
mpz/mul_si$U.lo mpz/mul_ui$U.lo \
|
||||
|
|
@ -195,23 +194,6 @@ MPQ_OBJECTS = mpq/abs$U.lo mpq/aors$U.lo \
|
|||
|
||||
MPN_OBJECTS = mpn/fib_table$U.lo mpn/mp_bases$U.lo
|
||||
|
||||
PRINTF_OBJECTS = \
|
||||
printf/asprintf$U.lo printf/asprntffuns$U.lo \
|
||||
printf/doprnt$U.lo printf/doprntf$U.lo printf/doprnti$U.lo \
|
||||
printf/fprintf$U.lo \
|
||||
printf/obprintf$U.lo printf/obvprintf$U.lo printf/obprntffuns$U.lo \
|
||||
printf/printf$U.lo printf/printffuns$U.lo \
|
||||
printf/snprintf$U.lo printf/snprntffuns$U.lo \
|
||||
printf/sprintf$U.lo printf/sprintffuns$U.lo \
|
||||
printf/vasprintf$U.lo printf/vfprintf$U.lo printf/vprintf$U.lo \
|
||||
printf/vsnprintf$U.lo printf/vsprintf$U.lo \
|
||||
printf/repl-vsnprintf$U.lo
|
||||
|
||||
SCANF_OBJECTS = \
|
||||
scanf/doscan$U.lo scanf/fscanf$U.lo scanf/fscanffuns$U.lo \
|
||||
scanf/scanf$U.lo scanf/sscanf$U.lo scanf/sscanffuns$U.lo \
|
||||
scanf/vfscanf$U.lo scanf/vscanf$U.lo scanf/vsscanf$U.lo
|
||||
|
||||
# no $U for C++ files
|
||||
CXX_OBJECTS = \
|
||||
cxx/isfuns.lo cxx/ismpf.lo cxx/ismpq.lo cxx/ismpz.lo cxx/ismpznw.lo \
|
||||
|
|
@ -247,7 +229,7 @@ libgmp_la_SOURCES = gmp-impl.h longlong.h randmt.h \
|
|||
EXTRA_libgmp_la_SOURCES = tal-debug.c tal-notreent.c tal-reent.c
|
||||
libgmp_la_DEPENDENCIES = @TAL_OBJECT@ \
|
||||
$(MPZ_OBJECTS) $(MPN_OBJECTS) @mpn_objs_in_libgmp@ \
|
||||
# $(PRINTF_OBJECTS) $(SCANF_OBJECTS) $(MPF_OBJECTS) $(MPQ_OBJECTS)
|
||||
# $(MPF_OBJECTS) $(MPQ_OBJECTS)
|
||||
libgmp_la_LIBADD = $(libgmp_la_DEPENDENCIES)
|
||||
libgmp_la_LDFLAGS = $(GMP_LDFLAGS) $(LIBGMP_LDFLAGS) \
|
||||
-version-info $(LIBGMP_LT_CURRENT):$(LIBGMP_LT_REVISION):$(LIBGMP_LT_AGE)
|
||||
|
|
|
|||
|
|
@ -110,20 +110,19 @@ am__DEPENDENCIES_1 = mpz/abs$U.lo mpz/add$U.lo mpz/add_ui$U.lo \
|
|||
mpz/getlimbn$U.lo mpz/import$U.lo mpz/init$U.lo mpz/init2$U.lo \
|
||||
mpz/inp_raw$U.lo mpz/inp_str$U.lo mpz/invert$U.lo mpz/ior$U.lo \
|
||||
mpz/iset$U.lo mpz/iset_d$U.lo mpz/iset_si$U.lo \
|
||||
mpz/iset_str$U.lo mpz/iset_ui$U.lo mpz/jacobi$U.lo \
|
||||
mpz/kronsz$U.lo mpz/kronuz$U.lo mpz/kronzs$U.lo \
|
||||
mpz/kronzu$U.lo mpz/lcm$U.lo mpz/lcm_ui$U.lo \
|
||||
mpz/lucnum_ui$U.lo mpz/lucnum2_ui$U.lo mpz/millerrabin$U.lo \
|
||||
mpz/mod$U.lo mpz/mul$U.lo mpz/mul_2exp$U.lo mpz/mul_si$U.lo \
|
||||
mpz/mul_ui$U.lo mpz/n_pow_ui$U.lo mpz/neg$U.lo \
|
||||
mpz/nextprime$U.lo mpz/out_raw$U.lo mpz/out_str$U.lo \
|
||||
mpz/perfpow$U.lo mpz/perfsqr$U.lo mpz/popcount$U.lo \
|
||||
mpz/pow_ui$U.lo mpz/powm$U.lo mpz/powm_ui$U.lo \
|
||||
mpz/pprime_p$U.lo mpz/random$U.lo mpz/random2$U.lo \
|
||||
mpz/realloc$U.lo mpz/realloc2$U.lo mpz/remove$U.lo \
|
||||
mpz/root$U.lo mpz/rootrem$U.lo mpz/rrandomb$U.lo \
|
||||
mpz/scan0$U.lo mpz/scan1$U.lo mpz/set$U.lo mpz/set_d$U.lo \
|
||||
mpz/set_f$U.lo mpz/set_q$U.lo mpz/set_si$U.lo mpz/set_str$U.lo \
|
||||
mpz/iset_str$U.lo mpz/iset_ui$U.lo mpz/lcm$U.lo \
|
||||
mpz/lcm_ui$U.lo mpz/lucnum_ui$U.lo mpz/lucnum2_ui$U.lo \
|
||||
mpz/millerrabin$U.lo mpz/mod$U.lo mpz/mul$U.lo \
|
||||
mpz/mul_2exp$U.lo mpz/mul_si$U.lo mpz/mul_ui$U.lo \
|
||||
mpz/n_pow_ui$U.lo mpz/neg$U.lo mpz/nextprime$U.lo \
|
||||
mpz/out_raw$U.lo mpz/out_str$U.lo mpz/perfpow$U.lo \
|
||||
mpz/perfsqr$U.lo mpz/popcount$U.lo mpz/pow_ui$U.lo \
|
||||
mpz/powm$U.lo mpz/powm_ui$U.lo mpz/pprime_p$U.lo \
|
||||
mpz/random$U.lo mpz/random2$U.lo mpz/realloc$U.lo \
|
||||
mpz/realloc2$U.lo mpz/remove$U.lo mpz/root$U.lo \
|
||||
mpz/rootrem$U.lo mpz/rrandomb$U.lo mpz/scan0$U.lo \
|
||||
mpz/scan1$U.lo mpz/set$U.lo mpz/set_d$U.lo mpz/set_f$U.lo \
|
||||
mpz/set_q$U.lo mpz/set_si$U.lo mpz/set_str$U.lo \
|
||||
mpz/set_ui$U.lo mpz/setbit$U.lo mpz/size$U.lo \
|
||||
mpz/sizeinbase$U.lo mpz/sqrt$U.lo mpz/sqrtrem$U.lo \
|
||||
mpz/sub$U.lo mpz/sub_ui$U.lo mpz/swap$U.lo mpz/tdiv_ui$U.lo \
|
||||
|
|
@ -407,7 +406,7 @@ LIBGMPXX_LT_AGE = 1
|
|||
LIBMP_LT_CURRENT = 4
|
||||
LIBMP_LT_REVISION = 10
|
||||
LIBMP_LT_AGE = 1
|
||||
SUBDIRS = mpn mpz # mpq cxx mpf printf scanf tests demos tune
|
||||
SUBDIRS = mpn mpz # mpq cxx mpf tests demos tune
|
||||
|
||||
# The ansi2knr setups for the build programs are the same as the normal
|
||||
# automake ansi2knr rules, but using $(CC_FOR_BUILD) instead of $(CC).
|
||||
|
|
@ -498,8 +497,7 @@ MPZ_OBJECTS = mpz/abs$U.lo mpz/add$U.lo mpz/add_ui$U.lo \
|
|||
mpz/import$U.lo mpz/init$U.lo mpz/init2$U.lo mpz/inp_raw$U.lo \
|
||||
mpz/inp_str$U.lo mpz/invert$U.lo \
|
||||
mpz/ior$U.lo mpz/iset$U.lo mpz/iset_d$U.lo mpz/iset_si$U.lo \
|
||||
mpz/iset_str$U.lo mpz/iset_ui$U.lo mpz/jacobi$U.lo mpz/kronsz$U.lo \
|
||||
mpz/kronuz$U.lo mpz/kronzs$U.lo mpz/kronzu$U.lo \
|
||||
mpz/iset_str$U.lo mpz/iset_ui$U.lo \
|
||||
mpz/lcm$U.lo mpz/lcm_ui$U.lo mpz/lucnum_ui$U.lo mpz/lucnum2_ui$U.lo \
|
||||
mpz/millerrabin$U.lo mpz/mod$U.lo mpz/mul$U.lo mpz/mul_2exp$U.lo \
|
||||
mpz/mul_si$U.lo mpz/mul_ui$U.lo \
|
||||
|
|
@ -532,23 +530,6 @@ MPQ_OBJECTS = mpq/abs$U.lo mpq/aors$U.lo \
|
|||
mpq/set_f$U.lo mpq/swap$U.lo
|
||||
|
||||
MPN_OBJECTS = mpn/fib_table$U.lo mpn/mp_bases$U.lo
|
||||
PRINTF_OBJECTS = \
|
||||
printf/asprintf$U.lo printf/asprntffuns$U.lo \
|
||||
printf/doprnt$U.lo printf/doprntf$U.lo printf/doprnti$U.lo \
|
||||
printf/fprintf$U.lo \
|
||||
printf/obprintf$U.lo printf/obvprintf$U.lo printf/obprntffuns$U.lo \
|
||||
printf/printf$U.lo printf/printffuns$U.lo \
|
||||
printf/snprintf$U.lo printf/snprntffuns$U.lo \
|
||||
printf/sprintf$U.lo printf/sprintffuns$U.lo \
|
||||
printf/vasprintf$U.lo printf/vfprintf$U.lo printf/vprintf$U.lo \
|
||||
printf/vsnprintf$U.lo printf/vsprintf$U.lo \
|
||||
printf/repl-vsnprintf$U.lo
|
||||
|
||||
SCANF_OBJECTS = \
|
||||
scanf/doscan$U.lo scanf/fscanf$U.lo scanf/fscanffuns$U.lo \
|
||||
scanf/scanf$U.lo scanf/sscanf$U.lo scanf/sscanffuns$U.lo \
|
||||
scanf/vfscanf$U.lo scanf/vscanf$U.lo scanf/vsscanf$U.lo
|
||||
|
||||
|
||||
# no $U for C++ files
|
||||
CXX_OBJECTS = \
|
||||
|
|
@ -585,7 +566,7 @@ libgmp_la_SOURCES = gmp-impl.h longlong.h randmt.h \
|
|||
EXTRA_libgmp_la_SOURCES = tal-debug.c tal-notreent.c tal-reent.c
|
||||
libgmp_la_DEPENDENCIES = @TAL_OBJECT@ \
|
||||
$(MPZ_OBJECTS) $(MPN_OBJECTS) @mpn_objs_in_libgmp@ \
|
||||
# $(PRINTF_OBJECTS) $(SCANF_OBJECTS) $(MPF_OBJECTS) $(MPQ_OBJECTS)
|
||||
# $(MPF_OBJECTS) $(MPQ_OBJECTS)
|
||||
|
||||
libgmp_la_LIBADD = $(libgmp_la_DEPENDENCIES)
|
||||
libgmp_la_LDFLAGS = $(GMP_LDFLAGS) $(LIBGMP_LDFLAGS) \
|
||||
|
|
|
|||
|
|
@ -1,314 +0,0 @@
|
|||
/* mpz_jacobi, mpz_legendre, mpz_kronecker -- mpz/mpz Jacobi symbols.
|
||||
|
||||
Copyright 2000, 2001, 2002, 2005 Free Software Foundation, Inc.
|
||||
|
||||
This file is part of the GNU MP Library.
|
||||
|
||||
The GNU MP Library is free software; you can redistribute it and/or modify it
|
||||
under the terms of the GNU Lesser General Public License as published by the
|
||||
Free Software Foundation; either version 2.1 of the License, or (at your
|
||||
option) any later version.
|
||||
|
||||
The GNU MP Library 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 Lesser General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public License along
|
||||
with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free
|
||||
Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301,
|
||||
USA. */
|
||||
|
||||
#include <stdio.h>
|
||||
#include "gmp.h"
|
||||
#include "gmp-impl.h"
|
||||
#include "longlong.h"
|
||||
|
||||
|
||||
/* Change this to "#define TRACE(x) x" for some traces. */
|
||||
#define TRACE(x)
|
||||
|
||||
|
||||
#define MPN_RSHIFT_OR_COPY(dst,src,size,shift) \
|
||||
do { \
|
||||
if ((shift) != 0) \
|
||||
{ \
|
||||
ASSERT_NOCARRY (mpn_rshift (dst, src, size, shift)); \
|
||||
(size) -= ((dst)[(size)-1] == 0); \
|
||||
} \
|
||||
else \
|
||||
MPN_COPY (dst, src, size); \
|
||||
} while (0)
|
||||
|
||||
|
||||
/* This code does triple duty as mpz_jacobi, mpz_legendre and mpz_kronecker.
|
||||
|
||||
mpz_jacobi could assume b is odd, but the improvements from that seem
|
||||
small compared to other operations, and anything significant should be
|
||||
checked at run-time since we'd like odd b to go fast in mpz_kronecker
|
||||
too.
|
||||
|
||||
mpz_legendre could assume b is an odd prime, but knowing this doesn't
|
||||
present any obvious benefits. Result 0 wouldn't arise (unless "a" is a
|
||||
multiple of b), but the checking for that takes little time compared to
|
||||
other operations.
|
||||
|
||||
The main loop is just a simple binary GCD with the jacobi symbol result
|
||||
tracked during the reduction.
|
||||
|
||||
The special cases for a or b fitting in one limb let mod_1 or modexact_1
|
||||
get used, without any copying, and end up just as efficient as the mixed
|
||||
precision mpz_kronecker_ui etc.
|
||||
|
||||
When tdiv_qr is called it's not necessary to make "a" odd or make a
|
||||
working copy of it, but tdiv_qr is going to be pretty slow so it's not
|
||||
worth bothering trying to save anything for that case.
|
||||
|
||||
Enhancements:
|
||||
|
||||
mpn_bdivmod could be used instead of mpn_tdiv_qr, like in mpn_gcd.
|
||||
Currently tdiv_qr is preferred since it's sub-quadratic on big sizes,
|
||||
although bdivmod might be a touch quicker on small sizes. This can be
|
||||
revised when bdivmod becomes sub-quadratic too.
|
||||
|
||||
Some sort of multi-step algorithm should be used. The current subtract
|
||||
and shift for every bit is very inefficient. Lehmer (per current gcdext)
|
||||
would need some low bits included in its calculation to apply the sign
|
||||
change for reciprocity. Binary Lehmer keeps low bits to strip twos
|
||||
anyway, so might be better suited. Maybe the accelerated GCD style k-ary
|
||||
reduction would work, if sign changes due to the extra factors it
|
||||
introduces can be accounted for (or maybe they can be ignored). */
|
||||
|
||||
|
||||
int
|
||||
mpz_jacobi (mpz_srcptr a, mpz_srcptr b)
|
||||
{
|
||||
mp_srcptr asrcp, bsrcp;
|
||||
mp_size_t asize, bsize;
|
||||
mp_ptr ap, bp;
|
||||
mp_limb_t alow, blow, ahigh, bhigh, asecond, bsecond;
|
||||
unsigned atwos, btwos;
|
||||
int result_bit1;
|
||||
TMP_DECL;
|
||||
|
||||
TRACE (printf ("start asize=%d bsize=%d\n", SIZ(a), SIZ(b));
|
||||
mpz_trace (" a", a);
|
||||
mpz_trace (" b", b));
|
||||
|
||||
asize = SIZ(a);
|
||||
asrcp = PTR(a);
|
||||
alow = asrcp[0];
|
||||
|
||||
bsize = SIZ(b);
|
||||
if (bsize == 0)
|
||||
return JACOBI_LS0 (alow, asize); /* (a/0) */
|
||||
|
||||
bsrcp = PTR(b);
|
||||
blow = bsrcp[0];
|
||||
|
||||
if (asize == 0)
|
||||
return JACOBI_0LS (blow, bsize); /* (0/b) */
|
||||
|
||||
/* (even/even)=0 */
|
||||
if (((alow | blow) & 1) == 0)
|
||||
return 0;
|
||||
|
||||
/* account for effect of sign of b, then ignore it */
|
||||
result_bit1 = JACOBI_BSGN_SS_BIT1 (asize, bsize);
|
||||
bsize = ABS (bsize);
|
||||
|
||||
/* low zero limbs on b can be discarded */
|
||||
JACOBI_STRIP_LOW_ZEROS (result_bit1, alow, bsrcp, bsize, blow);
|
||||
|
||||
count_trailing_zeros (btwos, blow);
|
||||
TRACE (printf ("b twos %u\n", btwos));
|
||||
|
||||
/* establish shifted blow */
|
||||
blow >>= btwos;
|
||||
if (bsize > 1)
|
||||
{
|
||||
bsecond = bsrcp[1];
|
||||
if (btwos != 0)
|
||||
blow |= (bsecond << (GMP_NUMB_BITS - btwos)) & GMP_NUMB_MASK;
|
||||
}
|
||||
|
||||
/* account for effect of sign of a, then ignore it */
|
||||
result_bit1 ^= JACOBI_ASGN_SU_BIT1 (asize, blow);
|
||||
asize = ABS (asize);
|
||||
|
||||
if (bsize == 1 || (bsize == 2 && (bsecond >> btwos) == 0))
|
||||
{
|
||||
/* special case one limb b, use modexact and no copying */
|
||||
|
||||
/* (a/2)=(2/a) with a odd, and if b is even then a is odd here */
|
||||
result_bit1 ^= JACOBI_TWOS_U_BIT1 (btwos, alow);
|
||||
|
||||
if (blow == 1) /* (a/1)=1 always */
|
||||
return JACOBI_BIT1_TO_PN (result_bit1);
|
||||
|
||||
JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, alow, asrcp, asize, blow);
|
||||
TRACE (printf ("base (%lu/%lu) with %d\n",
|
||||
alow, blow, JACOBI_BIT1_TO_PN (result_bit1)));
|
||||
return mpn_jacobi_base (alow, blow, result_bit1);
|
||||
}
|
||||
|
||||
/* Discard low zero limbs of a. Usually there won't be anything to
|
||||
strip, hence not bothering with it for the bsize==1 case. */
|
||||
JACOBI_STRIP_LOW_ZEROS (result_bit1, blow, asrcp, asize, alow);
|
||||
|
||||
count_trailing_zeros (atwos, alow);
|
||||
TRACE (printf ("a twos %u\n", atwos));
|
||||
result_bit1 ^= JACOBI_TWOS_U_BIT1 (atwos, blow);
|
||||
|
||||
/* establish shifted alow */
|
||||
alow >>= atwos;
|
||||
if (asize > 1)
|
||||
{
|
||||
asecond = asrcp[1];
|
||||
if (atwos != 0)
|
||||
alow |= (asecond << (GMP_NUMB_BITS - atwos)) & GMP_NUMB_MASK;
|
||||
}
|
||||
|
||||
/* (a/2)=(2/a) with a odd */
|
||||
result_bit1 ^= JACOBI_TWOS_U_BIT1 (btwos, alow);
|
||||
|
||||
if (asize == 1 || (asize == 2 && (asecond >> atwos) == 0))
|
||||
{
|
||||
/* another special case with modexact and no copying */
|
||||
|
||||
if (alow == 1) /* (1/b)=1 always */
|
||||
return JACOBI_BIT1_TO_PN (result_bit1);
|
||||
|
||||
/* b still has its twos, so cancel out their effect */
|
||||
result_bit1 ^= JACOBI_TWOS_U_BIT1 (btwos, alow);
|
||||
|
||||
result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow); /* now (b/a) */
|
||||
JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, blow, bsrcp, bsize, alow);
|
||||
TRACE (printf ("base (%lu/%lu) with %d\n",
|
||||
blow, alow, JACOBI_BIT1_TO_PN (result_bit1)));
|
||||
return mpn_jacobi_base (blow, alow, result_bit1);
|
||||
}
|
||||
|
||||
|
||||
TMP_MARK;
|
||||
TMP_ALLOC_LIMBS_2 (ap, asize, bp, bsize);
|
||||
|
||||
MPN_RSHIFT_OR_COPY (ap, asrcp, asize, atwos);
|
||||
ASSERT (alow == ap[0]);
|
||||
TRACE (mpn_trace ("stripped a", ap, asize));
|
||||
|
||||
MPN_RSHIFT_OR_COPY (bp, bsrcp, bsize, btwos);
|
||||
ASSERT (blow == bp[0]);
|
||||
TRACE (mpn_trace ("stripped b", bp, bsize));
|
||||
|
||||
/* swap if necessary to make a longer than b */
|
||||
if (asize < bsize)
|
||||
{
|
||||
TRACE (printf ("swap\n"));
|
||||
MPN_PTR_SWAP (ap,asize, bp,bsize);
|
||||
MP_LIMB_T_SWAP (alow, blow);
|
||||
result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
|
||||
}
|
||||
|
||||
/* If a is bigger than b then reduce to a mod b.
|
||||
Division is much faster than chipping away at "a" bit-by-bit. */
|
||||
if (asize > bsize)
|
||||
{
|
||||
mp_ptr rp, qp;
|
||||
|
||||
TRACE (printf ("tdiv_qr asize=%ld bsize=%ld\n", asize, bsize));
|
||||
|
||||
TMP_ALLOC_LIMBS_2 (rp, bsize, qp, asize-bsize+1);
|
||||
mpn_tdiv_qr (qp, rp, (mp_size_t) 0, ap, asize, bp, bsize);
|
||||
ap = rp;
|
||||
asize = bsize;
|
||||
MPN_NORMALIZE (ap, asize);
|
||||
|
||||
TRACE (printf ("tdiv_qr asize=%ld bsize=%ld\n", asize, bsize);
|
||||
mpn_trace (" a", ap, asize);
|
||||
mpn_trace (" b", bp, bsize));
|
||||
|
||||
if (asize == 0) /* (0/b)=0 for b!=1 */
|
||||
goto zero;
|
||||
|
||||
alow = ap[0];
|
||||
goto strip_a;
|
||||
}
|
||||
|
||||
for (;;)
|
||||
{
|
||||
ASSERT (asize >= 1); /* a,b non-empty */
|
||||
ASSERT (bsize >= 1);
|
||||
ASSERT (ap[asize-1] != 0); /* a,b normalized (and hence non-zero) */
|
||||
ASSERT (bp[bsize-1] != 0);
|
||||
ASSERT (alow == ap[0]); /* low limb copies should be correct */
|
||||
ASSERT (blow == bp[0]);
|
||||
ASSERT (alow & 1); /* a,b odd */
|
||||
ASSERT (blow & 1);
|
||||
|
||||
TRACE (printf ("top asize=%ld bsize=%ld\n", asize, bsize);
|
||||
mpn_trace (" a", ap, asize);
|
||||
mpn_trace (" b", bp, bsize));
|
||||
|
||||
/* swap if necessary to make a>=b, applying reciprocity
|
||||
high limbs are almost always enough to tell which is bigger */
|
||||
if (asize < bsize
|
||||
|| (asize == bsize
|
||||
&& ((ahigh=ap[asize-1]) < (bhigh=bp[asize-1])
|
||||
|| (ahigh == bhigh
|
||||
&& mpn_cmp (ap, bp, asize-1) < 0))))
|
||||
{
|
||||
TRACE (printf ("swap\n"));
|
||||
MPN_PTR_SWAP (ap,asize, bp,bsize);
|
||||
MP_LIMB_T_SWAP (alow, blow);
|
||||
result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
|
||||
}
|
||||
|
||||
if (asize == 1)
|
||||
break;
|
||||
|
||||
/* a = a-b */
|
||||
ASSERT (asize >= bsize);
|
||||
ASSERT_NOCARRY (mpn_sub (ap, ap, asize, bp, bsize));
|
||||
MPN_NORMALIZE (ap, asize);
|
||||
alow = ap[0];
|
||||
|
||||
/* (0/b)=0 for b!=1. b!=1 when a==0 because otherwise would have had
|
||||
a==1 which is asize==1 and would have exited above. */
|
||||
if (asize == 0)
|
||||
goto zero;
|
||||
|
||||
strip_a:
|
||||
/* low zero limbs on a can be discarded */
|
||||
JACOBI_STRIP_LOW_ZEROS (result_bit1, blow, ap, asize, alow);
|
||||
|
||||
if ((alow & 1) == 0)
|
||||
{
|
||||
/* factors of 2 from a */
|
||||
unsigned twos;
|
||||
count_trailing_zeros (twos, alow);
|
||||
TRACE (printf ("twos %u\n", twos));
|
||||
result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, blow);
|
||||
ASSERT_NOCARRY (mpn_rshift (ap, ap, asize, twos));
|
||||
asize -= (ap[asize-1] == 0);
|
||||
alow = ap[0];
|
||||
}
|
||||
}
|
||||
|
||||
ASSERT (asize == 1 && bsize == 1); /* just alow and blow left */
|
||||
TMP_FREE;
|
||||
|
||||
/* (1/b)=1 always (in this case have b==1 because a>=b) */
|
||||
if (alow == 1)
|
||||
return JACOBI_BIT1_TO_PN (result_bit1);
|
||||
|
||||
/* swap with reciprocity and do (b/a) */
|
||||
result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
|
||||
TRACE (printf ("base (%lu/%lu) with %d\n",
|
||||
blow, alow, JACOBI_BIT1_TO_PN (result_bit1)));
|
||||
return mpn_jacobi_base (blow, alow, result_bit1);
|
||||
|
||||
zero:
|
||||
TMP_FREE;
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -1,129 +0,0 @@
|
|||
/* mpz_si_kronecker -- long+mpz Kronecker/Jacobi symbol.
|
||||
|
||||
Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
|
||||
|
||||
This file is part of the GNU MP Library.
|
||||
|
||||
The GNU MP Library is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU Lesser General Public License as published by
|
||||
the Free Software Foundation; either version 2.1 of the License, or (at your
|
||||
option) any later version.
|
||||
|
||||
The GNU MP Library 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 Lesser General Public
|
||||
License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public License
|
||||
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
|
||||
the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
||||
MA 02110-1301, USA. */
|
||||
|
||||
#include "gmp.h"
|
||||
#include "gmp-impl.h"
|
||||
#include "longlong.h"
|
||||
|
||||
|
||||
int
|
||||
mpz_si_kronecker (long a, mpz_srcptr b)
|
||||
{
|
||||
mp_srcptr b_ptr;
|
||||
mp_limb_t b_low;
|
||||
mp_size_t b_size;
|
||||
mp_size_t b_abs_size;
|
||||
mp_limb_t a_limb, b_rem;
|
||||
unsigned twos;
|
||||
int result_bit1;
|
||||
|
||||
#if GMP_NUMB_BITS < BITS_PER_ULONG
|
||||
if (a > GMP_NUMB_MAX || a < -GMP_NUMB_MAX)
|
||||
{
|
||||
mp_limb_t alimbs[2];
|
||||
mpz_t az;
|
||||
ALLOC(az) = numberof (alimbs);
|
||||
PTR(az) = alimbs;
|
||||
mpz_set_si (az, a);
|
||||
return mpz_kronecker (az, b);
|
||||
}
|
||||
#endif
|
||||
|
||||
b_size = SIZ (b);
|
||||
if (b_size == 0)
|
||||
return JACOBI_S0 (a); /* (a/0) */
|
||||
|
||||
/* account for the effect of the sign of b, then ignore it */
|
||||
result_bit1 = JACOBI_BSGN_SS_BIT1 (a, b_size);
|
||||
|
||||
b_ptr = PTR(b);
|
||||
b_low = b_ptr[0];
|
||||
b_abs_size = ABS (b_size);
|
||||
|
||||
if ((b_low & 1) != 0)
|
||||
{
|
||||
/* b odd */
|
||||
|
||||
result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_low);
|
||||
a_limb = (unsigned long) ABS(a);
|
||||
|
||||
if ((a_limb & 1) == 0)
|
||||
{
|
||||
/* (0/b)=1 for b=+/-1, 0 otherwise */
|
||||
if (a_limb == 0)
|
||||
return (b_abs_size == 1 && b_low == 1);
|
||||
|
||||
/* a even, b odd */
|
||||
count_trailing_zeros (twos, a_limb);
|
||||
a_limb >>= twos;
|
||||
/* (a*2^n/b) = (a/b) * twos(n,a) */
|
||||
result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, b_low);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
/* (even/even)=0, and (0/b)=0 for b!=+/-1 */
|
||||
if ((a & 1) == 0)
|
||||
return 0;
|
||||
|
||||
/* a odd, b even
|
||||
|
||||
Establish shifted b_low with valid bit1 for ASGN and RECIP below.
|
||||
Zero limbs stripped are acounted for, but zero bits on b_low are
|
||||
not because they remain in {b_ptr,b_abs_size} for the
|
||||
JACOBI_MOD_OR_MODEXACT_1_ODD. */
|
||||
|
||||
JACOBI_STRIP_LOW_ZEROS (result_bit1, a, b_ptr, b_abs_size, b_low);
|
||||
if ((b_low & 1) == 0)
|
||||
{
|
||||
if (UNLIKELY (b_low == GMP_NUMB_HIGHBIT))
|
||||
{
|
||||
/* need b_ptr[1] to get bit1 in b_low */
|
||||
if (b_abs_size == 1)
|
||||
{
|
||||
/* (a/0x80000000) = (a/2)^(BPML-1) */
|
||||
if ((GMP_NUMB_BITS % 2) == 0)
|
||||
result_bit1 ^= JACOBI_TWO_U_BIT1 (a);
|
||||
return JACOBI_BIT1_TO_PN (result_bit1);
|
||||
}
|
||||
|
||||
/* b_abs_size > 1 */
|
||||
b_low = b_ptr[1] << 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
count_trailing_zeros (twos, b_low);
|
||||
b_low >>= twos;
|
||||
}
|
||||
}
|
||||
|
||||
result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_low);
|
||||
a_limb = (unsigned long) ABS(a);
|
||||
}
|
||||
|
||||
if (a_limb == 1)
|
||||
return JACOBI_BIT1_TO_PN (result_bit1); /* (1/b)=1 */
|
||||
|
||||
/* (a/b*2^n) = (b*2^n mod a / a) * recip(a,b) */
|
||||
JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, b_rem, b_ptr, b_abs_size, a_limb);
|
||||
result_bit1 ^= JACOBI_RECIP_UU_BIT1 (a_limb, b_low);
|
||||
return mpn_jacobi_base (b_rem, a_limb, result_bit1);
|
||||
}
|
||||
|
|
@ -1,121 +0,0 @@
|
|||
/* mpz_ui_kronecker -- ulong+mpz Kronecker/Jacobi symbol.
|
||||
|
||||
Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
|
||||
|
||||
This file is part of the GNU MP Library.
|
||||
|
||||
The GNU MP Library is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU Lesser General Public License as published by
|
||||
the Free Software Foundation; either version 2.1 of the License, or (at your
|
||||
option) any later version.
|
||||
|
||||
The GNU MP Library 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 Lesser General Public
|
||||
License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public License
|
||||
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
|
||||
the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
||||
MA 02110-1301, USA. */
|
||||
|
||||
#include "gmp.h"
|
||||
#include "gmp-impl.h"
|
||||
#include "longlong.h"
|
||||
|
||||
|
||||
int
|
||||
mpz_ui_kronecker (unsigned long a, mpz_srcptr b)
|
||||
{
|
||||
mp_srcptr b_ptr;
|
||||
mp_limb_t b_low;
|
||||
int b_abs_size;
|
||||
mp_limb_t b_rem;
|
||||
int twos;
|
||||
int result_bit1;
|
||||
|
||||
/* (a/-1)=1 when a>=0, so the sign of b is ignored */
|
||||
b_abs_size = ABSIZ (b);
|
||||
|
||||
if (b_abs_size == 0)
|
||||
return JACOBI_U0 (a); /* (a/0) */
|
||||
|
||||
if (a > GMP_NUMB_MAX)
|
||||
{
|
||||
mp_limb_t alimbs[2];
|
||||
mpz_t az;
|
||||
ALLOC(az) = numberof (alimbs);
|
||||
PTR(az) = alimbs;
|
||||
mpz_set_ui (az, a);
|
||||
return mpz_kronecker (az, b);
|
||||
}
|
||||
|
||||
b_ptr = PTR(b);
|
||||
b_low = b_ptr[0];
|
||||
result_bit1 = 0;
|
||||
|
||||
if (! (b_low & 1))
|
||||
{
|
||||
/* (0/b)=0 for b!=+/-1; and (even/even)=0 */
|
||||
if (! (a & 1))
|
||||
return 0;
|
||||
|
||||
/* a odd, b even
|
||||
|
||||
Establish shifted b_low with valid bit1 for the RECIP below. Zero
|
||||
limbs stripped are accounted for, but zero bits on b_low are not
|
||||
because they remain in {b_ptr,b_abs_size} for
|
||||
JACOBI_MOD_OR_MODEXACT_1_ODD. */
|
||||
|
||||
JACOBI_STRIP_LOW_ZEROS (result_bit1, a, b_ptr, b_abs_size, b_low);
|
||||
if (! (b_low & 1))
|
||||
{
|
||||
if (UNLIKELY (b_low == GMP_NUMB_HIGHBIT))
|
||||
{
|
||||
/* need b_ptr[1] to get bit1 in b_low */
|
||||
if (b_abs_size == 1)
|
||||
{
|
||||
/* (a/0x80...00) == (a/2)^(NUMB-1) */
|
||||
if ((GMP_NUMB_BITS % 2) == 0)
|
||||
{
|
||||
/* JACOBI_STRIP_LOW_ZEROS does nothing to result_bit1
|
||||
when GMP_NUMB_BITS is even, so it's still 0. */
|
||||
ASSERT (result_bit1 == 0);
|
||||
result_bit1 = JACOBI_TWO_U_BIT1 (a);
|
||||
}
|
||||
return JACOBI_BIT1_TO_PN (result_bit1);
|
||||
}
|
||||
|
||||
/* b_abs_size > 1 */
|
||||
b_low = b_ptr[1] << 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
count_trailing_zeros (twos, b_low);
|
||||
b_low >>= twos;
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (a == 0) /* (0/b)=1 for b=+/-1, 0 otherwise */
|
||||
return (b_abs_size == 1 && b_low == 1);
|
||||
|
||||
if (! (a & 1))
|
||||
{
|
||||
/* a even, b odd */
|
||||
count_trailing_zeros (twos, a);
|
||||
a >>= twos;
|
||||
/* (a*2^n/b) = (a/b) * (2/a)^n */
|
||||
result_bit1 = JACOBI_TWOS_U_BIT1 (twos, b_low);
|
||||
}
|
||||
}
|
||||
|
||||
if (a == 1)
|
||||
return JACOBI_BIT1_TO_PN (result_bit1); /* (1/b)=1 */
|
||||
|
||||
/* (a/b*2^n) = (b*2^n mod a / a) * RECIP(a,b) */
|
||||
JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, b_rem, b_ptr, b_abs_size, a);
|
||||
result_bit1 ^= JACOBI_RECIP_UU_BIT1 (a, b_low);
|
||||
return mpn_jacobi_base (b_rem, (mp_limb_t) a, result_bit1);
|
||||
}
|
||||
|
|
@ -1,86 +0,0 @@
|
|||
/* mpz_kronecker_si -- mpz+long Kronecker/Jacobi symbol.
|
||||
|
||||
Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
|
||||
|
||||
This file is part of the GNU MP Library.
|
||||
|
||||
The GNU MP Library is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU Lesser General Public License as published by
|
||||
the Free Software Foundation; either version 2.1 of the License, or (at your
|
||||
option) any later version.
|
||||
|
||||
The GNU MP Library 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 Lesser General Public
|
||||
License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public License
|
||||
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
|
||||
the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
||||
MA 02110-1301, USA. */
|
||||
|
||||
#include "gmp.h"
|
||||
#include "gmp-impl.h"
|
||||
#include "longlong.h"
|
||||
|
||||
|
||||
/* After the absolute value of b is established it's treated as an unsigned
|
||||
long, because 0x80..00 doesn't fit in a signed long. */
|
||||
|
||||
int
|
||||
mpz_kronecker_si (mpz_srcptr a, long b)
|
||||
{
|
||||
mp_srcptr a_ptr;
|
||||
mp_size_t a_size;
|
||||
mp_limb_t a_rem, b_limb;
|
||||
int result_bit1;
|
||||
|
||||
a_size = SIZ(a);
|
||||
if (a_size == 0)
|
||||
return JACOBI_0S (b);
|
||||
|
||||
#if GMP_NUMB_BITS < BITS_PER_ULONG
|
||||
if (b > GMP_NUMB_MAX || b < -GMP_NUMB_MAX)
|
||||
{
|
||||
mp_limb_t blimbs[2];
|
||||
mpz_t bz;
|
||||
ALLOC(bz) = numberof (blimbs);
|
||||
PTR(bz) = blimbs;
|
||||
mpz_set_si (bz, b);
|
||||
return mpz_kronecker (a, bz);
|
||||
}
|
||||
#endif
|
||||
|
||||
result_bit1 = JACOBI_BSGN_SS_BIT1 (a_size, b);
|
||||
b_limb = (unsigned long) ABS (b);
|
||||
a_ptr = PTR(a);
|
||||
|
||||
if ((b_limb & 1) == 0)
|
||||
{
|
||||
mp_limb_t a_low = a_ptr[0];
|
||||
int twos;
|
||||
|
||||
if (b_limb == 0)
|
||||
return JACOBI_LS0 (a_low, a_size); /* (a/0) */
|
||||
|
||||
if (! (a_low & 1))
|
||||
return 0; /* (even/even)=0 */
|
||||
|
||||
/* (a/2)=(2/a) for a odd */
|
||||
count_trailing_zeros (twos, b_limb);
|
||||
b_limb >>= twos;
|
||||
result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, a_low);
|
||||
}
|
||||
|
||||
if (b_limb == 1)
|
||||
return JACOBI_BIT1_TO_PN (result_bit1); /* (a/1)=1 for any a */
|
||||
|
||||
result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a_size, b_limb);
|
||||
a_size = ABS(a_size);
|
||||
|
||||
/* (a/b) = (a mod b / b) */
|
||||
JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, a_rem, a_ptr, a_size, b_limb);
|
||||
return mpn_jacobi_base (a_rem, b_limb, result_bit1);
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -1,80 +0,0 @@
|
|||
/* mpz_kronecker_ui -- mpz+ulong Kronecker/Jacobi symbol.
|
||||
|
||||
Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
|
||||
|
||||
This file is part of the GNU MP Library.
|
||||
|
||||
The GNU MP Library is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU Lesser General Public License as published by
|
||||
the Free Software Foundation; either version 2.1 of the License, or (at your
|
||||
option) any later version.
|
||||
|
||||
The GNU MP Library 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 Lesser General Public
|
||||
License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public License
|
||||
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
|
||||
the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
||||
MA 02110-1301, USA. */
|
||||
|
||||
#include "gmp.h"
|
||||
#include "gmp-impl.h"
|
||||
#include "longlong.h"
|
||||
|
||||
|
||||
int
|
||||
mpz_kronecker_ui (mpz_srcptr a, unsigned long b)
|
||||
{
|
||||
mp_srcptr a_ptr;
|
||||
mp_size_t a_size;
|
||||
mp_limb_t a_rem;
|
||||
int result_bit1;
|
||||
|
||||
a_size = SIZ(a);
|
||||
if (a_size == 0)
|
||||
return JACOBI_0U (b);
|
||||
|
||||
if (b > GMP_NUMB_MAX)
|
||||
{
|
||||
mp_limb_t blimbs[2];
|
||||
mpz_t bz;
|
||||
ALLOC(bz) = numberof (blimbs);
|
||||
PTR(bz) = blimbs;
|
||||
mpz_set_ui (bz, b);
|
||||
return mpz_kronecker (a, bz);
|
||||
}
|
||||
|
||||
a_ptr = PTR(a);
|
||||
if ((b & 1) != 0)
|
||||
{
|
||||
result_bit1 = JACOBI_ASGN_SU_BIT1 (a_size, b);
|
||||
}
|
||||
else
|
||||
{
|
||||
mp_limb_t a_low = a_ptr[0];
|
||||
int twos;
|
||||
|
||||
if (b == 0)
|
||||
return JACOBI_LS0 (a_low, a_size); /* (a/0) */
|
||||
|
||||
if (! (a_low & 1))
|
||||
return 0; /* (even/even)=0 */
|
||||
|
||||
/* (a/2)=(2/a) for a odd */
|
||||
count_trailing_zeros (twos, b);
|
||||
b >>= twos;
|
||||
result_bit1 = (JACOBI_TWOS_U_BIT1 (twos, a_low)
|
||||
^ JACOBI_ASGN_SU_BIT1 (a_size, b));
|
||||
}
|
||||
|
||||
if (b == 1)
|
||||
return JACOBI_BIT1_TO_PN (result_bit1); /* (a/1)=1 for any a */
|
||||
|
||||
a_size = ABS(a_size);
|
||||
|
||||
/* (a/b) = (a mod b / b) */
|
||||
JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, a_rem, a_ptr, a_size, b);
|
||||
return mpn_jacobi_base (a_rem, (mp_limb_t) b, result_bit1);
|
||||
}
|
||||
Loading…
Add table
Add a link
Reference in a new issue