My previous post looked at a problem that requires repeatedly finding the first digit of *k*^{n} where *k* is a single digit but *n* may be on the order of millions or billions.

The most direct approach would be to first compute *k*^{n} as a very large integer, then find it’s first digit. That approach is slow, and gets slower as *n* increases. A faster way is to look at the fractional part of log *k*^{n} = *n* log *k* and see which digit it corresponds to.

If *n* is not terribly big, this can be done in ordinary precision. But when *n* is large, multiplying log *k* by *n* and taking the fractional part brings less significant digits into significance. So for very large *n*, you need extra precision. I first did this in Python using SymPy, then switched to C++ for more speed. There I used the `quadmath`

library for `gcc`

. (If *n* is big enough, even quadruple precision isn’t enough. An advantage to SymPy over `quadmath`

is that the former has *arbitrary* precision. You could, for example, set the precision to be 10 more than the number of decimal places in *n* in order to retain 10 significant figures in the fractional part of *n* log *k*.)

The `quadmath.h`

header file needs to be wrapped in an `extern C `

declaration. Otherwise `gcc`

will give you misleading error messages.

The 128-bit floating point type `__float128`

has twice as many bits as a `double`

. The `quadmath`

functions have the same name as their standard `math.h`

counterparts, but with a `q`

added on the end, such as `log10q`

and `fmodq`

below.

Here’s code for computing the leading digit of *k*^{n} that illustrates using `quadmath`

.

#include <cmath > extern "C" { #include <quadmath.h> } __float128 logs[11]; for (int i = 2; i <= 10; i++) logs[i] = log10q(i + 0.0); int first_digit(int base, long long exponent) { __float128 t = fmodq(exponent*logs[base], 1.0); for (int i = 2; i <= 10; i++) if (t < logs[i]) return i-1; }

The code always returns because `t`

is less than 1.

Caching the values of `log10q`

saves repeated calls to a relatively expensive function. So does using the search at the bottom rather than computing `powq(10, t)`

.

The linear search at the end is more efficient than it may seem. First, it’s only search a list of length 9. Second, because of Benford’s law, the leading digits are searched in order of decreasing frequency, i.e. most inputs will cause `first_digit`

to return early in the search.

When you compile code using `quadmath`

, be sure to add `-lquadmath`

to the compile command.

Check out Boost.Multiprecision: http://boost.org/libs/multiprecision/

In addition to floating-point-fixed-precision types like, for instance, “cpp_dec_float_100” (floating point with 100 decimal digits precision; of course the “100” part is customizable), it also provides variable (not fixed at compile-time) precision types like “mpfr_float”.

Calculating at 1000 digits of precision can be fun! ;-)