logo标签列表

计算 1^X + 2^X + ... + N^X 模 1000000007

algorithmalgebranumber-theorymodulo
19

有没有算法可以计算 (1^x + 2^x + 3^x + ... + n^x) mod 1000000007？
注意：这里的a^b表示a的b次方。

限制条件为1 <= n <= 10^16, 1 <= x <= 1000，因此N的值非常大。

如果m = 1000000007，我只能解决O(m log m)的问题。但是由于时间限制为2秒，速度非常慢。

你有任何有效的算法吗？

有评论说这可能是这个问题的重复，但它肯定是不同的。

1
n^x mod m 等同于 ((n^(x-1) mod m) * (n mod m)) mod m; (a + b) mod m 等同于 ((a mod m) + (b mod m)) mod m。 - Vatine
@Vatine 但是n非常大！！！请阅读限制条件。如果我天真地计算1^x，2^x，3^x，...，n^x，那么计算时间将会很长！ - square1001
1
@CiaPan 请注意，10^9+7是一个质数。 - square1001
3
请查看https://en.wikipedia.org/wiki/Faulhaber's_formula。如果您的问题来自像HackerRank这样的挑战网站，我打赌内积总是可以简化为整数，并且有一种方法也可以摆脱全局因子。 - Fabian Pijcke
1
对于小的 n，您可以在合理的时间内直接计算总和。对于大的 n，我认为您可以找到这样的数字对，它们相加等于 1000000007（例如 1 和 1000000006，7 和 1000000000...）。这些数字分别等于 1 和（-1）mod 1000000007，因此它们在求和中被抵消 - 同样，它们各自的幂也会被抵消。然后您可以跳过它们，这样您只需直接计算不超过 1000000007 个总和项即可... - CiaPan
显示剩余13条评论
4个回答
19

你可以总结这个系列。

1**X + 2**X + ... + N**X

利用Faulhaber公式，可以得到一个具有X + 1次幂的多项式来计算任意的N。

如果您不想计算Bernoulli数，可以通过解决X + 2个线性方程（对于N = 1、N = 2、N = 3、...、N = X + 2）来找到这个多项式，虽然这种方法比较慢，但更容易实现。

让我们举一个X = 2的例子。在这种情况下，我们有一个X + 1 = 3阶多项式：

    A*N**3 + B*N**2 + C*N + D

线性方程为：

      A +    B +   C + D = 1              =  1 
    A*8 +  B*4 + C*2 + D = 1 + 4          =  5
   A*27 +  B*9 + C*3 + D = 1 + 4 + 9      = 14
   A*64 + B*16 + C*4 + D = 1 + 4 + 9 + 16 = 30

解决方程后，我们将得到

  A = 1/3
  B = 1/2
  C = 1/6
  D =   0

最终公式为

  1**2 + 2**2 + ... + N**2 == N**3 / 3 + N**2 / 2 + N / 6

现在，您需要做的就是将一个任意大的N放入公式中。到目前为止，该算法的复杂度为O(X**2)（因为它不依赖于N）。

但不幸的是，x可能是1000... 解决线性方程可以是O(X^3)，但我认为它很慢。我现在正在检查Faulhaber公式。 - square1001
如果你需要非常快的速度，我建议预计算伯努利数（最多到1000）；这不是一个简单的任务（如线性方程求解），但在这种情况下，你将得到几乎可以直接使用的多项式。 - Dmitry Bychenko
我查了一下伯努利数，最后明白了B[1]、B[2]、...、B[x]的数量级是O(x^2)。所以这个问题可以用O(N^2)的复杂度解决！ - square1001
@square1001：在“O(X ** 2)”中，要成立，它不取决于N（要求和的项目数），而是取决于它们应该被提高的幂次X。 - Dmitry Bychenko
2
已实现。我被接受了。 - square1001
显示剩余2条评论
3

有几种加速模幂运算的方法。从现在开始，我将使用**代表"求幂"，%代表"取余"。

首先，我们有一些观察结果。对于任意的a、b和m，恒有 (a * b) % m等于((a % m) * (b % m)) % m。同时，恒有a ** n等于(a ** floor(n / 2)) * (a ** (n - floor(n/2)))。这意味着，对于一个小于等于1000的指数，我们可以通过最多20次乘法（和21次取模）完成幂运算。

我们也可以跳过许多计算，因为(a ** b) % m等于((a % m) ** b) % m，如果m明显小于n，则我们只需要多次重复求和，并在末尾截取部分重复即可。

但是你的算法必须计算1^x，2^x，...，1000000006^x。它有超过10^9个术语，因此您必须迭代10^10次。所以，我认为它无法通过，因为时间限制为2秒。 - square1001
1
@square1001 刚在我的笔记本电脑上测试了一个 NOP 循环，循环10 ** 10次需要大约0.1秒，这可能主要是由 “fork” 和 “exec” 所占用的时间。 - Vatine
1
如果 x 是奇数，这就变得*容易多了，因为例如 1000000006 同余于 -1，所以 1^x 和 1000000006^x 相互抵消，其他很多项也是如此。但是如果 x 是偶数，我没有看到任何捷径。 - Tom Zych
如果x事先已知，您可以为指数的（非模）总和构建一个闭合公式。但是，是的，您可能会将其简化为“计算500000003 ** x％m;乘以floor（n / m）;循环从1到n％m，进行指数运算和求和”。 - Vatine
1
@Vatine 500000003被500000004取消，因为x是奇数。因此，对于奇数x和n=m，总和等于零。 - Sjoerd
1

我认为Vatine的答案可能是正确的方法，但我已经打出了这个答案，对于这个问题或其他类似的问题可能会有用。

今天早上我没有时间提供详细的答案，但请考虑以下内容。计算1^2 + 2^2 + 3^2 + ... + n^2需要O(n)步。然而，它等同于(n(n+1)(2n+1)/6)，可以在O(1)时间内计算出来。对于任何更高的整数幂x，都存在类似的等价关系。

可能存在一般性的解决方案；我不知道有没有，Wolfram Alpha似乎也不知道。然而，通常等价表达式的次数为x+1，可以通过计算一些示例值并解决一组线性方程来计算系数。

但是，对于像你的问题中的1000这样大的x，这将是困难的，并且可能无法在2秒钟内完成。

也许比我更懂数学的人可以将其转化为可行的解决方案？

编辑：糟糕，我看到Fabian Pijcke在我发帖之前已经发布了一个关于Faulhaber's formula的有用链接。

0
如果你想要一些易于实现且快速的东西，请尝试这个：
Function Sum(x: Number, n: Integer) -> Number
  P := PolySum(:x, n)
  return P(x)
End

Function PolySum(x: Variable, n: Integer) -> Polynomial
  C := Sum-Coefficients(n)
  P := 0
  For i from 1 to n + 1
    P += C[i] * x^i
  End
  return P
End

Function Sum-Coefficients(n: Integer) -> Vector of Rationals
  A := Create-Matrix(n)
  R := Reduced-Row-Echelon-Form(A)
  return last column of R
End

Function Create-Matrix(n: Integer) -> Matrix of Integers
  A := New (n + 1) x (n + 2) Matrix of Integers
  Fill A with 0s
  Fill first row of A with 1s

  For i from 2 to n + 1
    For j from i to n + 1
      A[i, j] := A[i-1, j] * (j - i + 2)
    End
    A[i, n+2] := A[i, n]
  End
  A[n+1, n+2] := A[n, n+2]
  return A
End

解释

我们的目标是获得一个多项式 Q，使得 Q(x) = sum i^n for i from 1 to x。已知 Q(x) = Q(x - 1) + x^n => Q(x) - Q(x - 1) = x^n，因此我们可以像下面这样建立一个方程组：

d^0/dx^0( Q(x) - Q(x - 1) ) = d^0/dx^0( x^n )
d^1/dx^1( Q(x) - Q(x - 1) ) = d^1/dx^1( x^n )
d^2/dx^2( Q(x) - Q(x - 1) ) = d^2/dx^2( x^n )
                           ...                            .
d^n/dx^n( Q(x) - Q(x - 1) ) = d^n/dx^n( x^n )

假设 Q(x) = a_1*x + a_2*x^2 + ... + a_(n+1)*x^(n+1)，我们将会有 n+1 个线性方程，未知数为 a1, ..., a_(n+1)，而且系数 cj 乘以方程 i 中的未知数 aj 遵循以下模式（其中 (k)_p = (k!)/(k - p)!）：
  • 如果 j < i，则 cj = 0
  • 否则，cj = (j)_(i - 1)
第 i 个方程的独立值为 (n)_(i - 1)。解释起来有点麻烦，但你可以在 这里 查看证明。

上述算法等同于解决这个线性方程组。在https://github.com/fcard/PolySum中可以找到大量的实现和进一步的解释。这个算法的主要缺点是它消耗了很多内存，即使我最节省内存的版本在n=3000时也使用了几乎1gb的内存。但它比SymPy和Mathematica都要快，所以我认为这没问题。与Schultz's method相比，后者使用了一组备选方程。

示例

对于小的n，手动应用这种方法很容易。当n=1时，矩阵如下：

| (1)_0   (2)_0   (1)_0 |     |   1       1       1   |
|   0     (2)_1   (1)_1 |  =  |   0       2       1   |

应用高斯-约旦消元法，我们得到

|   1       0      1/2  |
|   0       1      1/2  |

结果 = {a1 = 1/2, a2 = 1/2} => Q(x) = x/2 + (x^2)/2

请注意，矩阵总是已经处于行阶梯形式，我们只需要将其简化。

对于 n=2：

| (1)_0   (2)_0   (3)_0   (2)_0 |     |   1       1       1       1   |
|   0     (2)_1   (3)_1   (2)_1 |  =  |   0       2       3       2   |
|   0       0     (3)_2   (2)_2 |     |   0       0       6       2   |

应用高斯-约旦消元法，我们得到

|   1       1       0       2/3 |      |   1       0       0       1/6 |
|   0       2       0         1 |  =>  |   0       1       0       1/2 |
|   0       0       1       1/3 |      |   0       0       1       1/3 |

结果 = {a1 = 1/6, a2 = 1/2, a3 = 1/3} => Q(x) = x/6 + (x^2)/2 + (x^3)/3

该算法速度的关键在于它不会为矩阵的每个元素计算阶乘，而是知道 (k)_p = (k)_(p-1) * (k - (p - 1))，因此 A[i,j] = (j)_(i-1) = (j)_(i-2) * (j - (i - 2)) = A[i-1, j] * (j - (i - 2))，所以它使用上一行来计算当前行。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接