c++ - 仅需k位数字时的快速幂运算-继续

我需要帮助的地方...

我现在想做的是翻译此解决方案，该解决方案将数字的mantissa计算为c++:

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

从结果中查找前n位数字可以像这样完成:

"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)."

我的尝试(由数学和this code引发)失败了...:

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;
}

如果外面有人可以正确实现此解决方案，我需要您的帮助!

编辑

我的尝试输出示例:

n:2

m:0

n ^ m:1

计算的尾数:1.16334

n:2

m:1

n ^ m:2

计算的尾数:2.32667

n:2

m:2

n ^ m:4

计算的尾数:4.65335

n:2

m:98

n ^ m:3.16913e + 29

计算的尾数:8.0022

n:2

m:99

n ^ m:6.33825e + 29

计算的尾数:2.16596

最佳答案

我会避免使用pow。众所周知，正确实现很难。在很多SO问题中，标准库中的pow实现不当使人们为之困惑。

通过使用自然基数而不是基数10，您还可以避免很多麻烦。您将获得如下代码:

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 */

计算合适的m log(n)类似物。请注意，可能出现的最大m * logl(n)仅比2e10大一点。当您将其除以264并四舍五入到最接近的2的幂时，您会发现foo的ulp在最坏的情况下为2-29。特别是，这意味着即使使用完美的实现，使用long double也不能从该方法中获得超过8位的数字。
some_epsilon将是使long double始终超过数学上正确结果的最小expl(foo)；我还没有完全计算出来，但是应该在1e-9的数量级上。

鉴于这里的精度困难，我建议使用MPFR之类的库，而不是long double。您还可以使用double double技巧和四精度exp，log和fmod来使某些功能正常工作。