Faulhaber 公式的有效实现答案

【问题标题】：Efficient implementation fo Faulhaber's FormulaFaulhaber 公式的有效实现
【发布时间】：2014-03-29 03:44:28
【问题描述】：

我想要Faulhaber's Formula的高效实现

我想回答为

F(N,K) % P

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

注意：N 很大，最大可达 10^16，K 最大可达 3000

我在给定站点中尝试了double series implementation。但是对于非常大的 n 和 k 来说太耗时了。任何人都可以帮助提高此实现效率或描述一些其他方式来实现该公式。

【问题讨论】：

这是来自 Project Euler 吗？我怀疑您需要更深入地利用计算是 mod P 的事实。
是的，我想过，但伯努利数是分数，所以很难利用。但我不知道...
我想之前有人在这里评论过这个问题可能来自这里：hackerrank.com/contests/infinitum-mar14/challenges/…。在那里，K 最大为 10^3 (1000) 且 P = 10^9 + 7。
我想我成功了。请查看我的更新答案

标签： c++ performance algorithm combinations

【解决方案1】：

如何使用 Schultz (1980) 的想法，在您提到的双系列实现 (mathworld.wolfram.com/PowerSum.html) 下方概述？

来自 Wolfram MathWorld：

Schultz (1980) 表明，S_p(n) 的总和可以通过写作找到

并求解 p+1 方程组

为 j=0, 1, ..., p (Guo and 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)

【讨论】：

【解决方案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 版本实现了所有三个优化，对于 K=1000 大约需要 1 秒，对于 K=3000 大约需要 34 秒。（分别使用内存的60mb 和880mb）和完全未优化的python 实现对于K=1000 大约需要12 秒，但对于K=3000 内存不足。

无论使用何种语言，这种方法似乎都是最快的，但研究仍在进行中。由于舒尔茨的方法也归结为求解n+1 方程组，并且应该能够以相同的方式进行优化，这将取决于他的矩阵计算速度是否更快。此外，内存使用量根本无法很好地扩展，Mathematica 仍然是这里的明显赢家，仅使用 80mb 代替 K=3000。我们拭目以待。

【讨论】：