Anomaly in mpn_sqrtrem and mpn_rootrem
bodrato at mail.dm.unipi.it
bodrato at mail.dm.unipi.it
Wed Jul 8 10:11:10 UTC 2015
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/
More information about the gmp-devel
mailing list