Likely GMP bug

[Marco Bodrato]

Ciao,

Il Lun, 28 Maggio 2018 10:00 pm, Niels Möller ha scritto:
> I'd suggest the below (complete file, I think that's more readable

The code is clean.

You removed all gotos...

> The last part of the function requires vlimb odd, but tolerates
> arbitrary u, including 0.

... the effect is that in many cases (if U don't need reduction modulo V)
the trailing zeros of U are removed twice.

> This would be a candidate gcd_11 or gcd_11_odd.

Maybe, if we keep both parts in the same function, some code replication
can be good?

Is something like the following too redundant?

mp_limb_t
mpn_gcd_1 (mp_srcptr up, mp_size_t size, mp_limb_t vlimb)
{
  mp_limb_t      ulimb;
  unsigned long  zero_bits, u_low_zero_bits;
  int c;

  ASSERT (size >= 1);
  ASSERT (vlimb != 0);
  ASSERT_MPN_NONZERO_P (up, size);

  ulimb = up[0];

  /* Need vlimb odd for modexact, want it odd to get common zeros. */
  count_trailing_zeros (zero_bits, vlimb);
  vlimb >>= zero_bits;

  if (size > 1)
    {
      /* Must get common zeros before the mod reduction.  If ulimb==0 then
         vlimb already gives the common zeros.  */
      if (ulimb != 0)
        {
          count_trailing_zeros (u_low_zero_bits, ulimb);
          zero_bits = MIN (zero_bits, u_low_zero_bits);
        }

      ulimb = MPN_MOD_OR_MODEXACT_1_ODD (up, size, vlimb);
      if (ulimb == 0)
        return vlimb << zero_bits;

     /* Note that if ulimb == GMP_LIMB_HIGHBIT, c+1 is an invalid shift
count. */
      count_trailing_zeros (c, ulimb);
      ulimb = (ulimb >> 1) >> c;
    }
  else
    {
      /* size==1, so up[0]!=0 */
      count_trailing_zeros (u_low_zero_bits, ulimb);
      ulimb >>= u_low_zero_bits;
      zero_bits = MIN (zero_bits, u_low_zero_bits);

      /* make u bigger */
      if (vlimb > ulimb)
        MP_LIMB_T_SWAP (ulimb, vlimb);

      /* if u is much bigger than v, reduce using a division rather than
         chipping away at it bit-by-bit */
      if ((ulimb >> 16) > vlimb)
        {
          ulimb %= vlimb;
          if (ulimb == 0)
            return vlimb << zero_bits;

          count_trailing_zeros (c, ulimb);
          ulimb >>= (c + 1);
        {
      else
        ulimb >>= 1;
    }

  ASSERT (vlimb & 1);

  /* Represent the odd numbers ulimb and vlimb without the redundant
     least significant one bit. This reduction in size by one bit
     ensures that the high bit of t, below, is set if and only if
     vlimb > ulimb. And in addition, t != GMP_LIMB_HIGHBIT. */
  vlimb >>= 1;

  while (ulimb != vlimb)
    {
...



Ĝis,
m

-- 
http://bodrato.it/