Fast Exponentiation when only k digits are required - continued

Question

Where I need help...

What I want to do now is translate this solution, which calculates the mantissaof a number to c++:

n^m = exp10(m log10(n)) = exp(q (m log(n)/q)) where q = log(10)

Finding the first n digits from the result can be done like this:

"the first K digits of exp10(x) = the first K digits of exp10(frac(x))
 where frac(x) = the fractional part of x = x - floor(x)."

My attempts (sparked by the math and this code) failed...:

u l l function getPrefix(long double pow /*exponent*/, long double length /*length of prefix*/)
{
   long double dummy; //unused but necessary for modf
   long double q = log(10);

   u l l temp = floor(pow(10.0, exp(q * modf( (pow * log(2)/q), &dummy) + length - 1));
   return temp;
}

If anyone out there can correctly implement this solution, I need your help!!

---

EDIT

Example output from my attempts:

---

n: 2

m: 0

n^m: 1

Calculated mantissa: 1.16334

---

n: 2

m: 1

n^m: 2

Calculated mantissa: 2.32667

---

n: 2

m: 2

n^m: 4

Calculated mantissa: 4.65335

---

n: 2

m: 98

n^m: 3.16913e+29

Calculated mantissa: 8.0022

---

n: 2

m: 99

n^m: 6.33825e+29

Calculated mantissa: 2.16596

Answer 1

I'd avoid pow for this. It's notoriously hard to implement correctly. There are lots of SO questions where people got burned by a bad pow implementation in their standard library.

You can also save yourself a good deal of pain by working in the natural base instead of base 10. You'll get code that looks like this:

long double foo = m * logl(n);
foo = fmodl(foo, logl(10.0)) + some_epsilon;
sprintf(some_string, "%.9Lf", expl(foo));
/* boring string parsing code here */

to compute the appropriate analogue of m log(n). Notice that the largest m * logl(n) that can arise is just a little bigger than 2e10. When you divide that by 2⁶⁴ and round up to the nearest power of two, you see that an ulp of foo is 2^-29 at worst. This means, in particular, that you cannot get more than 8 digits out of this method using long doubles, even with a perfect implementation.

some_epsilon will be the smallest long double that makes expl(foo) always exceed the mathematically correct result; I haven't computed it exactly, but it should be on the order of 1e-9.

In light of the precision difficulties here, I might suggest using a library like MPFR instead of long doubles. You may also be able to get something to work using a double double trick and quad-precision exp, log, and fmod.