["Followup-To:" header set to comp.theory.]
On 2017-02-07, Julio P. Di Egidio <
ju...@diegidio.name> wrote:
> [Cross-posted to comp.theory, sci.math, comp.lang.prolog. Please set
> follow-ups as needed.]
>
> Hello all,
>
> I have put together an algorithm (plus reference implementation in Prolog)
> for the arbitrary precision integer binary logarithm with remainder: I
> needed this myself but could not find anything on the web, so I thought I'd
> share.
Pragmatically speaking, you probably won't run into an implementation of
arbitrary precision ("bignum") integers in which you cannot obtain an
estimate of this value from the implementation.
E.g. if a bignum consists of seventeen limbs of sixty-four bits,
then it consists of at least 16 x 64 bits, and you can apply a fast
algorithm on just the most significant 64 bit word to determine the
precise value.
Below is some C code for that which I hadded to a bignum library
called MPI (written by M. Fromberger).
We implement a binary search with bitmasks, rather than shifting,
which is practical due to the fixed size. The binary search tree
is expanded into open coded cases, so there is no loop, and all
the masks are literal constants:
static int s_highest_bit(mp_digit n)
{
#if MP_DIGIT_SIZE == 8
if (n & 0xFFFFFFFF00000000) {
if (n & 0xFFFF000000000000) {
if (n & 0xFF00000000000000) {
if (n & 0xF000000000000000) {
if (n & 0xC000000000000000)
return (n & 0x8000000000000000) ? 64 : 63;
else
return (n & 0x2000000000000000) ? 62 : 61;
} else {
if (n & 0x0C00000000000000)
return (n & 0x0800000000000000) ? 60 : 59;
else
return (n & 0x0200000000000000) ? 58 : 57;
}
} else {
if (n & 0x00F0000000000000) {
if (n & 0x00C0000000000000)
return (n & 0x0080000000000000) ? 56 : 55;
else
return (n & 0x0020000000000000) ? 54 : 53;
} else {
if (n & 0x000C000000000000)
return (n & 0x0008000000000000) ? 52 : 51;
else
return (n & 0x0002000000000000) ? 50 : 49;
}
}
} else {
if (n & 0x0000FF0000000000) {
if (n & 0x0000F00000000000) {
if (n & 0x0000C00000000000)
return (n & 0x0000800000000000) ? 48 : 47;
else
return (n & 0x0000200000000000) ? 46 : 45;
} else {
if (n & 0x00000C0000000000)
return (n & 0x0000080000000000) ? 44 : 43;
else
return (n & 0x0000020000000000) ? 42 : 41;
}
} else {
if (n & 0x000000F000000000) {
if (n & 0x000000C000000000)
return (n & 0x0000008000000000) ? 40 : 39;
else
return (n & 0x0000002000000000) ? 38 : 37;
} else {
if (n & 0x0000000C00000000)
return (n & 0x0000000800000000) ? 36 : 35;
else
return (n & 0x0000000200000000) ? 34 : 33;
}
}
}
} else {
if (n & 0x00000000FFFF0000) {
if (n & 0x00000000FF000000) {
if (n & 0x00000000F0000000) {
if (n & 0x00000000C0000000)
return (n & 0x0000000080000000) ? 32 : 31;
else
return (n & 0x0000000020000000) ? 30 : 29;
} else {
if (n & 0x000000000C000000)
return (n & 0x0000000008000000) ? 28 : 27;
else
return (n & 0x0000000002000000) ? 26 : 25;
}
} else {
if (n & 0x0000000000F00000) {
if (n & 0x0000000000C00000)
return (n & 0x0000000000800000) ? 24 : 23;
else
return (n & 0x0000000000200000) ? 22 : 21;
} else {
if (n & 0x00000000000C0000)
return (n & 0x0000000000080000) ? 20 : 19;
else
return (n & 0x0000000000020000) ? 18 : 17;
}
}
} else {
if (n & 0x000000000000FF00) {
if (n & 0x000000000000F000) {
if (n & 0x000000000000C000)
return (n & 0x0000000000008000) ? 16 : 15;
else
return (n & 0x0000000000002000) ? 14 : 13;
} else {
if (n & 0x0000000000000C00)
return (n & 0x0000000000000800) ? 12 : 11;
else
return (n & 0x0000000000000200) ? 10 : 9;
}
} else {
if (n & 0x00000000000000F0) {
if (n & 0x00000000000000C0)
return (n & 0x0000000000000080) ? 8 : 7;
else
return (n & 0x0000000000000020) ? 6 : 5;
} else {
if (n & 0x000000000000000C)
return (n & 0x0000000000000008) ? 4 : 3;
else
return (n & 0x0000000000000002) ? 2 : (n ? 1 : 0);
}
}
}
}
#elif MP_DIGIT_SIZE == 4
if (n & 0xFFFF0000) {
if (n & 0xFF000000) {
if (n & 0xF0000000) {
if (n & 0xC0000000)
return (n & 0x80000000) ? 32 : 31;
else
return (n & 0x20000000) ? 30 : 29;
} else {
if (n & 0x0C000000)
return (n & 0x08000000) ? 28 : 27;
else
return (n & 0x02000000) ? 26 : 25;
}
} else {
if (n & 0x00F00000) {
if (n & 0x00C00000)
return (n & 0x00800000) ? 24 : 23;
else
return (n & 0x00200000) ? 22 : 21;
} else {
if (n & 0x000C0000)
return (n & 0x00080000) ? 20 : 19;
else
return (n & 0x00020000) ? 18 : 17;
}
}
} else {
if (n & 0x0000FF00) {
if (n & 0x0000F000) {
if (n & 0x0000C000)
return (n & 0x00008000) ? 16 : 15;
else
return (n & 0x00002000) ? 14 : 13;
} else {
if (n & 0x00000C00)
return (n & 0x00000800) ? 12 : 11;
else
return (n & 0x00000200) ? 10 : 9;
}
} else {
if (n & 0x000000F0) {
if (n & 0x000000C0)
return (n & 0x00000080) ? 8 : 7;
else
return (n & 0x00000020) ? 6 : 5;
} else {
if (n & 0x0000000C)
return (n & 0x00000008) ? 4 : 3;
else
return (n & 0x00000002) ? 2 : (n ? 1 : 0);
}
}
}
#elif MP_DIGIT_SIZE == 2
if (n & 0xFF00) {
if (n & 0xF000) {
if (n & 0xC000)
return (n & 0x8000) ? 16 : 15;
else
return (n & 0x2000) ? 14 : 13;
} else {
if (n & 0x0C00)
return (n & 0x0800) ? 12 : 11;
else
return (n & 0x0200) ? 10 : 9;
}
} else {
if (n & 0x00F0) {
if (n & 0x00C0)
return (n & 0x0080) ? 8 : 7;
else
return (n & 0x0020) ? 6 : 5;
} else {
if (n & 0x000C)
return (n & 0x0008) ? 4 : 3;
else
return (n & 0x0002) ? 2 : (n ? 1 : 0);
}
}
#elif MP_DIGIT_SIZE == 1
if (n & 0xF0) {
if (n & 0xC0)
return (n & 0x80) ? 8 : 7;
else
return (n & 0x20) ? 6 : 5;
} else {
if (n & 0x0C)
return (n & 0x08) ? 4 : 3;
else
return (n & 0x02) ? 2 : (n ? 1 : 0);
}
#else
#error fixme: unsupported MP_DIGIT_SIZE
#endif
/* notreached */
abort();
}
With this, the function for the bits in the whole bignum is
just:
int s_highest_bit_mp(mp_int *a)
{
int nd1 = USED(a) - 1;
return s_highest_bit(DIGIT(a, nd1)) + nd1 * MP_DIGIT_BIT;
}