Fast exponentiation when only k digits are required - continued

Question

Where do I need help ...

Now I want to translate this solution, which calculates mantissa from number to C ++:

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

The search for the first n digits from the result can be performed as follows:

 "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 (caused by math and this code ) failed ...:

 ull 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); ull temp = floor(pow(10.0, exp(q * modf( (pow * log(2)/q), &dummy) + length - 1)); return temp; } 

If someone can correctly implement this solution, I need your help !!

---

EDIT

Example output from my attempts:

---

n: 2

m: 0

n ^ m: 1

Estimated Mantissa: 1.16334

---

n: 2

m: 1

n ^ m: 2

Estimated Mantissa: 2.32667

---

n: 2

m: 2

n ^ m: 4

Estimated Mantissa: 4.65335

---

n: 2

m: 98

n ^ m: 3.16913e + 29

Estimated Mantissa: 8.0022

---

n: 2

m: 99

n ^ m: 6.33825e + 29

Estimated Mantissa: 2.16596

Answer 1

I would avoid pow for this. This, as you know, is difficult to implement correctly. There are many SO questions in which people were burned due to poor implementation of pow in their standard library.

You can also save a lot of pain by working in a natural base instead of base 10. You will get a 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 corresponding analogue of m log(n) . Note that the largest m * logl(n) that can occur is slightly larger than 2e10 . When you divide this by 2 ⁶⁴ and round to the nearest power of two, you will see that ulp foo in the worst case is 2 ^-29 . This means, in particular, that you cannot get more than 8 digits from this method using long double s, even with an ideal implementation.

some_epsilon will be the smallest long double that makes expl(foo) always superior to the mathematically correct result; I did not calculate it exactly, but it should be of the order of 1e-9 .

In light of the difficulties with precision, I can suggest using a library like MPFR instead of long double s. You can also get something to work with the double double trick and the four-point precision exp , log and fmod .