Anomaly in mpn_sqrtrem and mpn_rootrem

[bodrato at mail.dm.unipi.it]

Ciao,

Il Mar, 7 Luglio 2015 11:12 am, Torbjörn Granlund ha scritto:
> bodrato at mail.dm.unipi.it writes:

>   Maybe we should improve mpn_rootrem for
>   small sizes in general...

> Tabling start values is hard, but we should consider doing it for k <
> some limit, and perhaps provide just 4 bits.

We can use something based on the logarithm.
The code above should give a range for the first 5 bits (the first is 1
anyway), and it should work for any k.

void
root (mp_limb_t *r, mp_limb_t n, unsigned k)
{
  /* vector(32,i,floor((log(31+i)/log(2)-5)*256)) */
  unsigned char v[] = {0, 11, 22, 33, 43, 53, 63, 73,
		 	82, 91, 100, 109, 117, 125, 134, 141,
			149, 157, 164, 172, 179, 186, 193, 200,
			206, 213, 219, 225, 232, 238, 244, 250, 255};
  unsigned c;
  mp_limb_t low, hi;

  count_leading_zeros (c, n);
  c -= GMP_NAIL_BITS;
  n <<= c;
  n >>= GMP_NUMB_BITS - 6;
  n -= 32;
  /* Now 0 <= n < 32 contains the higest bits of the original n, with
     the first 1 removed. */
  c = GMP_NUMB_BITS - c-1;
  /* c is the number of bits of the original number. */
  low = v[n] + (c%k) * 256;
  hi = v[n+1] + (c%k) * 256;
  low = low / k;
  hi = hi / k;
  for (c=32; low < v[--c];) ;
  r[0] = c+32;
  for (;hi > v[c]; ++c) ;
  r[1] = c+32;
}

Regards,
m

-- 
http://bodrato.it/papers/