Faulhaber 公式的高效实施

Question

我想要有效实施 Faulhaber 公式

我想要的答案是

F(N,K) % P

其中 F(N,K) 是 faulhaber 公式的实现，P 是素数。

注意：N非常大，可达10^16，K可达3000

我在给定站点中尝试了双系列实现。但对于非常大的 n 和 k 来说，它太耗时了。任何人都可以帮助提高此实现的效率或描述实现该公式的其他方法吗？

Answer 1

如何使用 Schultz (1980) 的想法，下面概述了您提到的双级数实现 (mathworld.wolfram.com/PowerSum.html)？

来自 Wolfram MathWorld：

    Schultz (1980) 证明总和

S_p(n)

可以通过写

求出

    

    并求解 p+1 方程组

    

   针对 j=0, 1, ..., p 获得（Guo 和 Qi 1999），其中

delta (j,p)

是 Kronecker delta。

下面是 Haskell 中的一次尝试，似乎有效。在我的旧笔记本电脑上，它会在大约 36 秒内返回 n=10^16、p=1000 的结果。

{-# OPTIONS_GHC -O2 #-}

import Math.Combinatorics.Exact.Binomial
import Data.Ratio
import Data.List (foldl')

kroneckerDelta a b | a == b    = 1 % 1
                   | otherwise = 0 % 1

g a b = ((-1)^(a - b +1) * choose a b) % 1

coefficients :: Integral a => a -> a -> [Ratio a] -> [Ratio a]
coefficients p j cs
  | j < 0     = cs
  | otherwise = coefficients p (j - 1) (c:cs)
 where
   c = f / g (j + 1) j
   f = foldl h (kroneckerDelta j p) (zip [j + 2..p + 1] cs)
   h accum (i,cf) = accum - g i j * cf

trim r = let n = numerator r
             d = denominator r
             l = div n d
         in (mod l (10^9 + 7),(n - d * l) % d)

s n p = numerator (a % 1 + b) where
 (a,b) = foldl' (\(i',r') (i,r) -> (mod (i' + i) (10^9 + 7),r' + r)) (0,0) 
      (zipWith (\c i ->  trim (c * n^i)) (coefficients p p []) [1..p + 1])

main = print (s (10^16) 1000)

Answer 2

我发现了自己的算法来计算从 Faulhaber 公式获得的多项式的系数；它、它的证明和几个实现可以在 github.com/fcard/PolySum 找到。这个问题启发我包含一个 C++ 实现（使用 GMP 库来实现任意精度数字），截至撰写本文时，减去几个可用性功能，它是：

#include <gmpxx.h>
#include <vector>

namespace polysum {
  typedef std::vector<mpq_class> mpq_row;
  typedef std::vector<mpq_class> mpq_column;
  typedef std::vector<mpq_row>   mpq_matrix;

  mpq_matrix make_matrix(size_t n) {
    mpq_matrix A(n+1, mpq_row(n+2, 0));
    A[0] = mpq_row(n+2, 1);

    for (size_t i = 1; i < n+1; i++) {
      for (size_t j = i; j < n+1; j++) {
        A[i][j] += A[i-1][j];
        A[i][j] *= (j - i + 2);
      }
      A[i][n+1] = A[i][n-1];
    }
    A[n][n+1] = A[n-1][n+1];
    return A;
  }

  void reduced_row_echelon(mpq_matrix& A) {
    size_t n = A.size() - 1;
    for (size_t i = n; i+1 > 0; i--) {
      A[i][n+1] /= A[i][i];
      A[i][i] = 1;
      for (size_t j = i-1; j+1 > 0; j--) {
        auto p = A[j][i];
        A[j][i] = 0;
        A[j][n+1] -= A[i][n+1] * p;
      }
    }
  }

  mpq_column sum_coefficients(size_t n) {
    auto A = make_matrix(n);
    reduced_row_echelon(A);

    mpq_column result;
    for (auto row: A) {
      result.push_back(row[n+1]);
    }
    return result;
  }
}

我们可以像这样使用上面的内容：

#include <cmath>
#include <gmpxx.h>
#include <polysum.h>

mpq_class power_sum(size_t K, unsigned int N) {
  auto coeffs = polysum::sum_coefficients(K)

  mpq_class result(0);
  for (size_t i = 0; i <= K; i++) {
    result += A[i][n+1] * pow(N, i+1);
  }
  return result;
}

完整的实现提供了一个可打印和可调用的

Polynomial

类，以及一个

polysum

函数来将一个多项式构造为另一个多项式的和。

#include "polysum.h"

void power_sum_print(size_t K, unsigned int N) {
  auto F = polysum::polysum(K);
  std::cout << "polynomial: " << F;
  std::cout << "result: " << F(N);
}

至于效率，上面在我的计算机上计算出

K=1000

和

N=1e16

的结果大约需要

1.75

秒，相比之下，更成熟和优化的SymPy实现在同一台机器上需要大约

90

秒，和 mathematica 需要

30

秒。对于

K=3000

，上面大约需要

4

分钟，mathematica 几乎花了

20

分钟（但使用的内存少得多），我让 sympy 运行了一整夜，但它没有完成，可能是因为它耗尽了内存.

这里可以进行的优化包括使矩阵稀疏并利用仅需要计算一半的行和列的事实。链接存储库中的 Rust 版本实现了稀疏和行优化，大约需要

0.7

秒来计算

K=1000

，大约需要

45

来计算

K=3000

（分别使用

105mb

和

2.9gb

内存） 。 Haskell 版本实现了所有三种优化，并且

1

大约需要

K=1000

秒，

34

大约需要

K=3000

秒。 （分别使用

60mb

和

880mb

内存）并且完全未优化的 Python 实现大约需要

12

秒（对于

K=1000

）但会耗尽内存

K=3000

。

无论使用哪种语言，这种方法看起来都是最快的，但研究仍在进行中。由于舒尔茨的方法也归结为求解

n+1

方程组，并且应该能够以相同的方式进行优化，因此这将取决于他的矩阵是否计算得更快。此外，内存使用根本无法很好地扩展，Mathematica 仍然是这里明显的赢家，仅使用

80mb

来表示

K=3000

。我们拭目以待。

Answer 3

我从 Faulhaber 错过的公式构建了一个简单的 Julia 脚本（就像他错过的学科中其他 20,000 个公式一样，后来发现了这些公式，对他来说很丢脸）。

function faulhaber(d)
    polyVector = [big(0)//big(1) for n = 1:d+2]
    #initialisation at d = 1, indices shifted by 1!
    polyVector[2] = 1/2
    polyVector[3] = 1/2
    for j=2:d
        oldUnit = polyIntEval(polyVector, 1)
        for k=j+2:-1:3
            polyVector[k] = j * polyVector[k-1]/(k-1)
        end
        polyVector[2] = 1 - j*oldUnit
    end
    return polyVector
end

function polyIntEval(v, x)
    d = length(v) - 1
    result = 0
    for j=0:d
        result += v[j+1]*x^(j+1) / (j+1)
    end
    return result
end

在我的笔记本电脑上，大约需要 5 秒才能获得 d = 1000 的完整多项式。