How to get the first x leading binary digits of 5**x without big integer multiplication

Viewed 130

I want to efficiently and elegantly compute with perfect precision the first x leading binary digits of 5**x?

For example 5**20 is 10101101011110001110101111000101101011000110001. The first 8 leading binary digits is 10101101.

In my use case, x is only up to 1-60. I don't want to create a table. A solution using 64-bit integers would be fine. I just don't want to use big integers.

2 Answers

first x leading binary digits of 5**x without big integer multiplication

efficiently and elegantly compute with perfect precision the first x leading binary digits of 5x?

"compute with perfect precision" leaves out pow(). Too many implementations will return an imperfect result and FP math might not use 64 bit precision, even with long double.


Form an integer with a 64-bit whole number part .ms and a 64-bit fraction part .ls. Then loop 60 times, multiply by 5 and diving by 2 as needed, to keep the leading bits from growing too big.

Note there is some precision lost in the fraction, with N > 42, yet that is not significant enough to affect the whole number part OP is seeking.

#include <inttypes.h>
#include <stdint.h>
#include <stdio.h>

typedef struct {
  uint64_t ms, ls;
} uint128;

// Simplifications possible here, leave for OP
uint128 times5(uint128 x) {
  uint128 y = x;
  for (int i=1; i<5; i++) {
    // y += x
    y.ms += x.ms;  
    y.ls += x.ls;
    if (y.ls < x.ls) y.ms++;
  }
  return y;
}

uint128 div2(uint128 x) { 
  x.ls = (x.ls >> 1) | (x.ms << 63);
  x.ms >>= 1;
  return x;
}

int main(void) {
  uint128 y = {.ms = 1};
  uint64_t pow2 = 2;
  for (unsigned x = 1; x <= 60; x++) {
    y = times5(y);
    while (y.ms >= pow2) {
      y = div2(y);
    }
    printf("%2u %16" PRIX64 ".%016" PRIX64 "\n", x, y.ms, y.ls);
    pow2 <<= 1;
  }
}

Output

         whole part.fraction
 1                1.4000000000000000
 2                3.2000000000000000
 3                7.D000000000000000
 4                9.C400000000000000
...
57  14643E5AE44D12B.8F5FEE5AA432560D
58  32FA9BE33AC0AEC.E66FD3E29A7DD720
59  7F7285B812E1B50.401791B6823A99D0
60  9F4F2726179A224.501D762422C94044
    ^-------------^ This is the part OP is seeking.

The key to solving this task is: divide and conquer. Form an algorithm, (which is simply *5 and /2 as needed), and code a type and functions to do each small step.


Is a loop of 60 efficient? Perhaps not. Another approach would use Exponentiation by squaring. Certainly would be worth it for large N, yet for N == 60, a loop was simple enough for a quick turn.

5n = 2(-n) • 10n

Using this identity, we can easily compute the leading N base-2 digits of (the nearest integer to) any given power of 5.

This code example is in C, but it's the same idea in any other language.

Example output: https://wandbox.org/permlink/Fs205DDzQR0gaLSo

#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdint.h>

#define STATIC_ASSERT(CONDITION) ((void)sizeof(int[(CONDITION) ? 1 : -1]))

uint64_t pow5_leading_digits(double power, uint8_t ndigits)
{
    STATIC_ASSERT(DBL_MANT_DIG <= 64);

    double pow5 = exp2(-power) * pow(10, power);

    const double binary_digits = ceil(log2(pow5));

    assert(ndigits <= DBL_MANT_DIG);

    if (!ndigits || binary_digits < 0)
        return 0;

    // If pow5 can fit in the number of digits requested, return it
    if (binary_digits <= ndigits)
        return pow5;
    // If pow5 is too big to return, divide by 2 until it fits
    if (binary_digits > DBL_MANT_DIG)
        pow5 /= exp2(binary_digits - DBL_MANT_DIG + 1);

    return (uint64_t)pow5 >> (DBL_MANT_DIG - ndigits);
}

Edit: Now limits the returned value to those exactly representable with double's.

Related