4

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

Jonathan
  • 1,126
  • 10
  • 19
  • 2
    Define "failed". Can you give some example input, the output you got, and the output you expected? – Oliver Charlesworth Jul 06 '13 at 20:46
  • Maybe this is more suitable for mathematics part of the stack exchange. – huseyin tugrul buyukisik Jul 06 '13 at 20:46
  • The math isn't the problem, its converting it into code. If that still applies to stack exchange I'll be happy to transfer my question. – Jonathan Jul 06 '13 at 20:48
  • I run a for loop that passes the index of the loop to the function as the exponent (pow), and a variable input by the user as the prefix length desired (length). What I get back are prefixes to numbers that aren't even powers of 2 (which should have been hard coded into the function according to the math). I can't see where I'm going wrong, so I can't explain how the output correlates to the input. – Jonathan Jul 06 '13 at 21:11
  • 2
    @Jonathan: You should find the first index for which you get the incorrect result, and then break that out as a standalone [test case](http://sscce.org) that we might be able to help with. – Oliver Charlesworth Jul 06 '13 at 22:15
  • I've verified that passing in i=0 through i=n as the power calculates all of the powers of two using `long double x = exp((m*log(2)/log(10))*log(10));`. But I still can't manage to arrive at the correct decimal values using the extra multiplications in the equation. – Jonathan Jul 06 '13 at 23:12
  • Thank you all for attempting to help me here. For now I'm putting this one to rest. If I find a solution I'll be sure to post it here! **Edit** _The problem is that all of the solutions suggested to me, and that I've found are limited to 32/64 bit numbers. I'll be working with arbitrarily large numbers (like powers of 2^n where n > 10000). I'm interested in solutions that are compatible with numbers that large._ – Jonathan Jul 07 '13 at 01:11

1 Answers1

3

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 264 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.

tmyklebu
  • 13,915
  • 3
  • 28
  • 57
  • Thank you for that explanation. I won't lie... half of that didn't make much sense to me yet, but I'll look into this! And I've noticed the precision issue. In the end, I'll probably have to find another way to solve my problem... But seriously, thanks again! Ps. I'm sorry I can't up your answer... – Jonathan Jul 07 '13 at 00:43
  • Is there any robust way to deal with the edge case where result starts with `9999999999999999999` or `10000000000000000`? Arbitrarily small relative error can change result. – zch Jul 08 '13 at 11:17
  • @zch: No. The absurdly small relative errors like that *probably* don't happen in this case, but I don't have a rigorous proof of why. (The recent proof of Catalan's conjecture shows that `999...9` never comes up, though --- this leads me to suspect that anything better than a heuristic argument for why this hack works given enough precision is going to be really, really hard.) – tmyklebu Jul 08 '13 at 15:36