Big Data/Analytics Zone is brought to you in partnership with:

John Cook is an applied mathematician working in Houston, Texas. His career has been a blend of research, software development, consulting, and management. John is a DZone MVB and is not an employee of DZone and has posted 175 posts at DZone. You can read more from them at their website. View Full User Profile

Leading Digits and Quadmath

06.25.2013
| 2824 views |
  • submit to reddit

My previous post on Gelfand's Question looked at a problem that requires repeatedly finding the first digit of kn wherek is a single digit but n may be on the order of millions or billions.

The most direct approach would be to first compute kn 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 kn = 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 logk 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 arbitraryprecision. 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 gccwill give you misleading error messages.

The 128-bit floating point type __float128 has twice as many bits as a double. The quadmathfunctions 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 kn 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.

Published at DZone with permission of John Cook, author and DZone MVB. (source)

(Note: Opinions expressed in this article and its replies are the opinions of their respective authors and not those of DZone, Inc.)