mpn_jacobi_base
Niels Möller
nisse at lysator.liu.se
Wed Mar 10 17:21:37 CET 2010
For fun, I've adapted the branch-free binary gcd loop from mpn_gcd_1 to
mpn_jacobi_base. Seems to work, disabled by default, and pushed into the
repository. Included below.
Depends on a reasonably fast count_trailing_zeros. One could try
replacing that with a loop or with a table lookup. Although the
corresponding table-based code in gcd_1.c didn't seem so promising, but
maybe it uses a too small table. Considerations should be the same for
gcd and jacobi.
On shell:
$ ./speed -c -s10-60 -t10 mpn_jacobi_base_1 mpn_jacobi_base_2 mpn_jacobi_base_3 mpn_jacobi_base_4
overhead 6.06 cycles, precision 10000 units of 3.74e-10 secs, CPU freq 2672.99 MHz
mpn_jacobi_base_1 mpn_jacobi_base_2 mpn_jacobi_base_3 mpn_jacobi_base_4
10 #104.73 143.50 146.80 112.92
20 240.51 383.43 370.81 #217.87
30 364.91 610.03 591.08 #315.81
40 492.92 844.08 810.62 #406.76
50 624.25 1081.43 1037.55 #500.09
60 756.88 1316.02 1259.01 #593.46
Regards,
/Niels
/* Computes (a/b) for odd b and any a. The initial bit is taken as a
* parameter. We have no need for the convention that the sign is in
* bit 1, internally we use bit 0. */
/* FIXME: Could try table-based count_trailing_zeros. */int
mpn_jacobi_base (mp_limb_t a, mp_limb_t b, int bit)
{
int c;
ASSERT (b & 1);
if (a == 0)
/* Ok, here we use that the sign is bit 1, after all. */
return b == 1 ? (1-(bit & 2)) : 0;
bit >>= 1;
/* Below, we represent a and b shifted right so that the least
significant one bit is implicit. */
b >>= 1;
count_trailing_zeros (c, a);
bit ^= c & (b ^ (b >> 1));
/* We may have c==GMP_LIMB_BITS-1, so we can't use a>>c+1. */
a >>= c;
a >>= 1;
while (a != b)
{
mp_limb_t t = a - b;
mp_limb_t bgta = LIMB_HIGHBIT_TO_MASK (t);
/* If b > a, invoke reciprocity */
bit ^= (bgta & a & b);
/* b <-- min (a, b) */
b += (bgta & t);
/* a <-- |a - b| */
a = (t ^ bgta) - bgta;
/* Number of trailing zeros is the same no matter if we look at
* t or a, but using t gives more parallelism. */
count_trailing_zeros (c, t);
c ++;
/* (2/b) = -1 if b = 3 or 5 mod 8 */
bit ^= c & (b ^ (b >> 1));
a >>= c;
}
return a == 0 ? 1-2*(bit & 1) : 0;
}
--
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.
More information about the gmp-devel
mailing list