mpn_jacobi_base

[Niels Möller]

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.