准确计算 1 到 10 之间值的谐波数

Question

我一直在用 C++ 开发一个项目，作为项目的一部分，我需要计算实数值的谐波数。对于超过 40 的值，我有一个非常准确的公式，并且使用 Kahan 和我可以准确地获得整数的谐波数，但我不知道如何准确计算像 HarmonicNumber(1.5) 这样的值。我该怎么做？

注意：如果有人有用于计算 Digamma 函数的快速 C++ 代码，我也可以使用它，因为 Digamma 函数很容易转换为 HarmonicNumber 函数。

编辑：我已经编写了以下代码，尽管我希望有更快、更简单的代码。

编辑 2：对于双精度值，我需要一个小于 10^-15 的相对误差，对于像“long double”这样的 80 位扩展精度值，我需要一个小于 10^-18 的相对误差

long double HarmonicNumber(long double n)
{
//Absolute error is smaller than or equal to 2^-61, or 4.33681e-19, for n<5000
//For n>5000, all but the last two bits are correct. 
constexpr long double m1 = 1.0L / 24;
constexpr long double m2 = -7.0L / 960;
constexpr long double m3 = 31.0L / 8064;
constexpr long double m4 = -127.0L / 30720;
constexpr long double m5 = 511.0L / 67584;
constexpr long double m6 = -1414477.0L / 67092480;
constexpr long double m7 = 8191.0L / 98304;
constexpr long double EulerGamma = 0.5772156649015328606065120900824024310421L;
long double v = n + 0.5L;
long double v2 = 1.0L / (v * v);
//Uses asymptotic expansion with progressively more terms.
//Fewer terms are needed for larger inputs.
if(n >= 10000.L) return m1*v2 + log(v) + EulerGamma;
if(n >= 450.00L) return (m2*v2 + m1)*v2 + log(v) + EulerGamma;
if(n >= 110.00L) return ((m3*v2 + m2)*v2 + m1)*v2 + log(v) + EulerGamma;
if(n >= 42.000L) return (((m4*v2 + m3)*v2 + m2)*v2 + m1)*v2 + log(v) + EulerGamma;
if(n >= 24.000L) return ((((m5*v2 + m4)*v2 + m3)*v2 + m2)*v2 + m1)*v2 + log(v) + EulerGamma;
if(n >= 17.000L) return (((((m6*v2 + m5)*v2 + m4)*v2 + m3)*v2 + m2)*v2 + m1)*v2 + log(v) + EulerGamma;
if(n >= 13.000L) return ((((((m7*v2 + m6)*v2 + m5)*v2 + m4)*v2 + m3)*v2 + m2)*v2 + m1)*v2 + log(v) + EulerGamma;
if(n >= 6.0L)
{
    //Calculates HarmonicNumber(n+7) and then subtracts fraction to find HarmonicNumber(n)
    v = n + 7.5L;
    v2 = 1.0L / (v * v);
    auto base = ((((((m7*v2 + m6)*v2 + m5)*v2 + m4)*v2 + m3)*v2 + m2)*v2 + m1)*v2 + log(v) + EulerGamma;
    long double n2 = n + 4.0L;
    auto n2sq = n2*n2;
    auto n2sh = n2sq - 7.0L;
    auto shft = 1.0L / n2 + (2*n2*n2sh*(3.0L*n2sq-7.0L))/(n2sq*n2sh*n2sh-36.0L);
    return base - shft;
}
else
{
    //Calculates HarmonicNumber(n+14) and then subtracts fraction to find HarmonicNumber(n)
    v = n + 14.5L;
    v2 = 1.0L / (v * v);
    auto base = ((((((m7*v2 + m6)*v2 + m5)*v2 + m4)*v2 + m3)*v2 + m2)*v2 + m1)*v2 + log(v) + EulerGamma;
    long double n2 = n + 4.0L;
    auto n2sq = n2*n2;
    auto n2sh = n2sq - 7.0L;
    auto shft = 1.0L / n2 + (2*n2*n2sh*(3.0L*n2sq-7.0L))/(n2sq*n2sh*n2sh-36.0L);
    n2 = n + 11.0L;
    n2sq = n2*n2;
    n2sh = n2sq - 7.0L;
    shft += 1.0L / n2 + (2*n2*n2sh*(3.0L*n2sq-7.0L))/(n2sq*n2sh*n2sh-36.0L);
    return base - shft;
}

}`

Answer 1

通常harmonic number 函数被理解为只为整数定义。正如您所指出的，您可以使用关系 H_{n} = \psi(n+1) + \gamma 转换为为非整数定义的 digamma 函数。您可以将其称为将调和数函数泛化为非整数，但如果您不希望与您交谈的数学人士挠头并交叉地看着您，您最好只是说您想要计算 digamma 函数。

因此，您希望计算小值和大值的 digamma 函数。幸运的是，Wikipedia tells you how。我在这里简单总结一下。对于较大的 x，您希望使用伯努利数的渐近展开式。

对于所有 x >~ 16，您可以使用前十个伯努利数的存储表获得全双精度。对于较小的 x，只需使用

将 x 移到足够大的值以应用第一种技术。例如，对于 x = 12.5，只需求和 s = 1/12.5 + 1/13.5 + 1/14.5 + 1/15.5。然后通过前面的技术计算 \psi(16.5) 并减去 s 得到 \psi(12.5)。

在我自己的library 中，我对小 x 使用了一种更复杂的技术，这种技术速度更快，消除错误的影响更小，但差异非常小。这种简单的技术对于大多数用途来说已经足够了。

最后，对于否定参数，以及非常接近零的参数，您应该使用reflection formula

转化为积极的论点。

将所有这些放在一起，您会得到以下代码：

static const int bernoulli_length = 8;
static const double bernoulli[] = {
    1.0 / 6.0, -1.0 / 30.0, 1.0 / 42.0, -1.0 / 30.0,
    5.0 / 66.0, -691.0 / 2730.0, 7.0 / 6.0, -3617.0 / 510.0
};

double Psi(double x) {

    // Reflect to positive x
    if (x < 0.25) {
        // For x very close to a negative integer, this will loose accuracy
        // due to finite PI. To fix this, we need sinpi and cospi functions. 
        return (Psi(1.0 - x) - M_PI * cos(M_PI * x) / sin(M_PI * x));
    }

    // Shift out to large enough x
    double s = 0.0;
    while (x < 16.0) {
        s += 1.0 / x;
        x += 1.0;
    }

    // Use the asymptotic expansion
    double psi = log(x) - 1.0 / (2.0 * x);
    double x2 = x * x;
    double x2k = 1.0;
    for (int k = 0; k < bernoulli_length; k++) {
        double psi_old = psi;
        x2k *= x2;
        psi -= bernoulli[k] / (2 * (k + 1) * x2k);
        if (psi == psi_old) {
            return(psi - s);
        }
    }
    throw std::range_error("Convergence failure.");
}

我会说这更清楚一点。