c++ - 使用 Boost 计算 Champernowne 常数 C10

Question

我正在尝试使用以下公式计算 Champernowne constant C10:

在上面的公式中，我用 b 代替 10 来计算 C10。我希望能够使用 Boost 的 cpp_dec_float 计算任意精度的常量。

这是我的代码:

#include <boost/multiprecision/cpp_dec_float.hpp>

const long long PRECISION = 100;

typedef boost::multiprecision::number<
    boost::multiprecision::cpp_dec_float<PRECISION> > arbFloat;

arbFloat champernowne()
{
    arbFloat c, sub, n, k;
    std::string precomp_c, postcomp_c;

    for(n = 1; n == 1 || precomp_c != postcomp_c; ++n) {
        for(k = 1; k <= n; ++k) {
            sub += floor(log10(k));
        }

        precomp_c = static_cast<std::string>(c);
        c += n / pow(10, n + sub);
        postcomp_c = static_cast<std::string>(c);
    }

    return c;
}

下面是代码的分解:

1. 我首先定义一个精度为 100 位的变量 arbFloat(这经常更改 — 所以我不想使用 cpp_dec_float_100)。

2. 公式有两个求和 block ，所以我使用两个 for 循环来实现它们。在最里面的 for 循环中，我计算以 k = 1 开始的总和，条件是 k <= n 的 floor(log10(k)) 。

我已验证在 floor() 上使用 log10() 和 cpp_dec_float 会返回具有正确精度的变量。

1. 因为最外层的求和一直到无穷大，我不得不在某个时候停止计算。为了检查是否超出精度，我在计算 c 之前将 c += n / pow(10, n + sub) 转换为字符串 - 然后在计算之后将其转换为字符串。如果字符串相同，我将结束计算，因为已超出精度(进一步的计算将是多余的)。

我还使用了这个设置(通过字符串转换和比较来检查超出的精度)来计算其他变量 - 它工作得很好。

1. 接下来我计算 c += n / pow(10, n + sub) 的最外层总和 - 以这种方式使用 pow() 确实保持了精度。最后，我返回 c 。

当我运行这个程序时，我得到以下变量:

0.1234567891001100120001300001400000150000001600000001700000000180000000001900000000002000000000000210

对比真正的 Champerowne 常数 C10:

0.1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253546

只有前 11 位数字是正确的，其余的都不是。我无法找到我要去哪里错了。我尝试了以下方法:

• 尝试用 c += n / pow(10, n + sub) 替换 c += n / pow(static_cast<arbFloat>(10), n + sub) 以检查 pow() 是否没有保持精度 - 但它没有改变任何东西。

• 尝试将 floor() 替换为将 log10(k) 转换为字符串并“舍入”字符串(仅保留 . 之前的字符)的方法 - 但它没有改变任何内容。

• 尝试将 k <= n 更改为 k < n k <= n + 1 - 以防万一我误解了总和 - 但这只会使它更加不准确。

如果我需要解释更多，请告诉我。任何帮助将不胜感激!

Answer 1

sub 的前一个值在每次迭代中被继承；在循环内声明它。

arbFloat champernowne() {
    arbFloat c;

    for (int n = 1;; ++n) {
        arbFloat sub;

        for (int k = 1; k <= n; ++k) {
            sub += floor(log10(k));
        }

        arbFloat const last = c;
        c += n / pow(10, n + sub);

        if (c == last) {
            break;
        }
    }

    return c;
}