计算几何级数之和 (mod m)

Question

我有一个系列，

S = i^(m) + i^(2m) + ...............  + i^(km)  (mod m)   

0 <= i < m, k may be very large (up to 100,000,000),  m <= 300000

我想求总和。我无法应用几何级数（GP）公式，因为结果将具有分母，然后我将必须找到可能不存在的模逆（如果分母和 m 不是互质）。

所以我做了一个替代算法，假设这些幂将使一个周期的长度远小于 k （因为它是一个模方程，所以我会得到类似 2,7,9,1,2,7,9， 1....) 并且该循环将在上述系列中重复。因此，我不会从 0 迭代到 k，而是只求一个循环中的数字之和，然后计算上述系列中的循环数并将它们相乘。所以我首先找到了i^m (mod m)，然后一次又一次地乘以这个数字，并在每一步取模，直到再次到达第一个元素。

但是当我实际编写算法时，对于 i 的某些值，我得到了非常大的循环。因此在终止之前花费了大量时间，因此我的假设是不正确的。

那么我们还能发现其他模式吗？ （基本上我不想迭代 k。） 所以请给我一个求和的有效算法的想法。

Answer 1

这是我遇到的类似问题的算法

您可能知道可以在对数时间内计算一个数字的幂。您也可以这样做来计算几何级数的总和。因为它认为

1 + a + a^2 + ... + a^(2*n+1) = (1 + a) * (1 + (a^2) + (a^2)^2 + ... + (a^2)^n),

可以递归计算右手上的几何级数来得到结果。

这样您就不需要除法，因此您可以将总和（以及中间结果）的余数对您想要的任何数字取模。

Answer 2

正如您所指出的，计算任意模 m 很困难，因为许多值可能没有乘法逆模 m。但是，如果您可以求解一组精心选择的替代模数，则可以将它们组合起来以获得模 m 的解。

将 m 分解为 p_1, p_2, p_3 ... p_n，使得每个 p_i 都是不同素数的幂

由于每个 p 都是不同素数的幂，因此它们是成对互质的。如果我们可以计算关于每个模数 p_i 的级数之和，我们可以使用中国剩余定理 将它们重新组装成 mod m 的解。

对于每个质数幂模，有两个简单的特殊情况：

如果 i^m 同余 0 mod p_i，则总和为 0。

如果 i^m 同余 1 mod p_i，那么总和等于 k ​​mod p_i。

对于其他值，可以应用几何数列之和的常用公式：

S = sum(j=0 to k, (i^m)^j) = ((i^m)^(k+1) - 1) / (i^m - 1)

TODO：证明 (i^m - 1) 与 p_i 互质，或者当它们具有非平凡的 GCD 时找到替代解决方案。希望 p_i 是素数幂并且也是 m 的约数这一事实会有一些用处......如果 p_i 是 i 的约数。条件成立。如果 p_i 是素数（与素数幂相对），则适用特殊情况 i^m = 1，或者 (i^m - 1) 具有乘法逆元。

如果几何和公式不适用于某些 p_i，您可以重新安排计算，这样您只需从 1 迭代到 p_i 而不是 1 到 k，从而利用事实上，这些术语的重复周期不超过 p_i。

（由于您的级数不包含 aj=0 项，因此您想要的值实际上是 S-1。）

这会产生一组同余 mod p_i，满足 CRT 的要求。
将它们组合成一个解 mod m 的过程在上面的链接中有描述，所以这里不再重复。

Answer 3

这可以通过重复平方的方法来完成，即O(log( k)) 时间，或 O(log(k)log(m)) 时间（如果您将 m 视为变量）。

一般来说，a[n]=1+b+b^2+... b^(n-1) mod m 可以通过以下方式计算：

1. a[j+k ]==b^{j}a[k]+a[j]
2. a[2n]==(b^n+1)a[n]

第二个是第一个推论。

在您的情况下，可以在 O(log m) 时间内计算出 b=i^m 。

下面的 Python 代码实现了这一点：

def geometric(n,b,m):
    T=1
    e=b%m
    total = 0
    while n>0:
        if n&1==1:
            total = (e*total + T)%m
        T = ((e+1)*T)%m
        e = (e*e)%m
        n = n/2
        //print '{} {} {}'.format(total,T,e)
    return total

这个魔法有一个数学原因 - 定义为结合的对的运算

(a,r)@(b,s)=(ab,as+r)

，规则 1 基本上意味着：

(b,1)@(b,1)@... n times ... @(b,1)=(b^n,1+b+b^2+...+b^(n-1))

当运算具有结合性时，重复平方总是有效。在本例中，@ 运算符的时间为 O(log(m))，因此重复平方需要 O(log(n)log(m))< /代码>。

看待这个问题的一种方法是矩阵求幂：

[[b,1],[0,1]]^n == [[b^n,1+b+...+b^(n-1))],[0,1]]

您可以使用类似的方法来计算 (a^nb^n)/(ab) 模 m 因为矩阵求幂给出：

[[b,1],[0,a]]^n == [[b^n,a^(n-1)+a^(n-2)b+...+ab^(n-2)+b^(n-1)],[0,a^n]]

Answer 4

基于 @braindoper 的方法，

1 + a + a^2 + ... +a^n mod m

在 Mathematica 中计算的完整算法如下所示：

geometricSeriesMod[a_, n_, m_] := 
   Module[ {q = a, exp = n, factor = 1, sum = 0, temp},

   While[And[exp > 0, q != 0],
     If[EvenQ[exp],
       temp = Mod[factor*PowerMod[q, exp, m], m];
       sum = Mod[sum + temp, m];
       exp--];
     factor = Mod[Mod[1 + q, m]*factor, m];
     q = Mod[q*q, m];
     exp = Floor[ exp /2];
   ];

   Return [Mod[sum + factor, m]]
]

参数：

• a 是级数的“比率”。它可以是任何整数（包括零和负值）。
• n 是该级数的最高指数。允许的整数 >= 0。
• m是整数模 != 0

注意：算法在每个算术运算后执行 Mod 运算。如果您将此算法转录为整数字长有限的语言，这一点至关重要。