===================================================================
RCS file: /home/cvs/OpenXM_contrib/gmp/mpz/Attic/kronuz.c,v
retrieving revision 1.1.1.1
retrieving revision 1.1.1.2
diff -u -p -r1.1.1.1 -r1.1.1.2
--- OpenXM_contrib/gmp/mpz/Attic/kronuz.c	2000/09/09 14:12:54	1.1.1.1
+++ OpenXM_contrib/gmp/mpz/Attic/kronuz.c	2003/08/25 16:06:33	1.1.1.2
@@ -1,7 +1,6 @@
-/* mpz_ui_kronecker -- Kronecker/Jacobi symbol. */
+/* mpz_ui_kronecker -- ulong+mpz Kronecker/Jacobi symbol.
 
-/*
-Copyright (C) 1999, 2000 Free Software Foundation, Inc.
+Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
 
 This file is part of the GNU MP Library.
 
@@ -18,98 +17,93 @@ 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., 59 Temple Place - Suite 330, Boston,
-MA 02111-1307, USA.
-*/
+MA 02111-1307, USA. */
 
 #include "gmp.h"
 #include "gmp-impl.h"
 #include "longlong.h"
 
 
+/* This implementation depends on BITS_PER_MP_LIMB being even, so that
+   (a/2)^BITS_PER_MP_LIMB = 1 and so there's no need to pay attention to how
+   many low zero limbs are stripped.  */
+#if BITS_PER_MP_LIMB % 2 != 0
+Error, error, unsupported BITS_PER_MP_LIMB
+#endif
+
+
 int
-#if __STDC__
 mpz_ui_kronecker (unsigned long a, mpz_srcptr b)
-#else
-mpz_ui_kronecker (a, b)
-     unsigned long a;
-     mpz_srcptr    b;
-#endif
 {
-  int        b_abs_size;
-  mp_srcptr  b_ptr;
-  mp_limb_t  b_low;
+  mp_srcptr  b_ptr = PTR(b);
+  mp_limb_t  b_low = b_ptr[0];
+  mp_size_t  b_abs_size = ABSIZ (b);
+  mp_limb_t  b_rem;
   int        twos;
-  int        result_bit1;
+  int        result_bit1 = 0;
 
-  /* (a/0) */
-  b_abs_size = ABSIZ (b);
-  if (b_abs_size == 0)
-    return JACOBI_U0 (a);
-
   /* (a/-1)=1 when a>=0, so the sign of b is ignored */
-  b_ptr = PTR(b);
-  b_low = b_ptr[0];
 
-  /* (0/1)=1; (0/-1)=1; (0/b)=0 for b!=+/-1
-     (1/b)=1, for any b */
-  if (a <= 1)
-    return (a == 1) | ((b_abs_size == 1) & (b_low == 1));
+  if (b_abs_size == 0)
+    return JACOBI_U0 (a);  /* (a/0) */
 
-  if (b_low & 1)
+  if (a > GMP_NUMB_MAX)
     {
-      /* (a/1) = 1 for any a */
-      if (b_abs_size == 1 && b_low == 1)
-        return 1;
-
-      count_trailing_zeros (twos, a);
-      a >>= twos;
-      if (a == 1)
-        return JACOBI_TWOS_U (twos, b_low);  /* powers of (2/b) only */
-
-      /* powers of (2/b); reciprocity to (b/a); (b/a) == (b mod a / a) */
-      return mpn_jacobi_base (mpn_mod_1 (b_ptr, b_abs_size, a),
-                              a,
-                              JACOBI_TWOS_U_BIT1 (twos, b_low)
-                              ^ JACOBI_RECIP_UU_BIT1 (b_low, a));
+      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 is even; (a/2)=0 if a is even */
-  if ((a & 1) == 0)
-    return 0;
+  if (! (b_low & 1))
+    {
+      /* (0/b)=0 for b!=+/-1; (even/even)=0 */
+      if (! (a & 1))
+        return 0;
 
-  /* Require MP_BITS_PER_LIMB even, so (a/2)^MP_BITS_PER_LIMB = 1, and so we
-     don't have to pay attention to how many trailing zero limbs are
-     stripped.  */
-  ASSERT ((BITS_PER_MP_LIMB & 1) == 0);
+      /* a odd, b even */
 
-  MPN_STRIP_LOW_ZEROS_NOT_ZERO (b_ptr, b_abs_size);
-  b_low = b_ptr[0];
+      /* establish shifted b_low for use with reciprocity below */
+      MPN_STRIP_LOW_ZEROS_NOT_ZERO (b_ptr, b_abs_size, b_low);
+      if (! (b_low & 1))
+        {
+          if (b_low == GMP_NUMB_HIGHBIT)
+            {
+              if (b_abs_size == 1)   /* (a/0x80000000) == (a/2)^(BPML-1) */
+                return JACOBI_TWOS_U (GMP_NUMB_BITS-1, a);
 
-  if (b_low & 1)
-    /* reciprocity to (b/a); (b/a) == (b mod a / a) */
-    return mpn_jacobi_base (mpn_mod_1 (b_ptr, b_abs_size, a),
-                            a,
-                            JACOBI_RECIP_UU_BIT1 (b_low, a));
-
-  count_trailing_zeros (twos, b_low);
-
-  /* reciprocity to get (b/a) */
-  if (twos == BITS_PER_MP_LIMB-1)
+              /* b_abs_size > 1 */
+              b_low = b_ptr[1] << 1;
+            }
+          else
+            {
+              count_trailing_zeros (twos, b_low);
+              b_low >>= twos;
+            }
+        }
+    }
+  else
     {
-      if (b_abs_size == 1)
+      if (a == 0)        /* (0/b)=1 for b=+/-1, 0 otherwise */
+        return (b_abs_size == 1 && b_low == 1);
+
+      if (! (a & 1))
         {
-          /* b==0x800...00, one limb high bit only, so (a/2)^(BPML-1) */
-          return JACOBI_TWOS_U (BITS_PER_MP_LIMB-1, a);
+          /* 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);
         }
-
-      /* b_abs_size > 1 */
-      result_bit1 = JACOBI_RECIP_UU_BIT1 (a, b_ptr[1] << 1);
     }
-  else
-    result_bit1 = JACOBI_RECIP_UU_BIT1 (a, b_low >> twos);
 
-  /* powers of (a/2); reciprocity to (b/a); (b/a) == (b mod a / a) */
-  return mpn_jacobi_base (mpn_mod_1_rshift (b_ptr, b_abs_size, twos, a),
-                          a,
-                          JACOBI_TWOS_U_BIT1 (twos, a) ^ result_bit1);
+  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);
 }