18

我有一个系列

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。）所以请给我一个有效算法的想法来找到总和。

4

4 回答 4

15

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

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

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

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

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

于 2012-02-28T17:21:17.320 回答
13

正如您所指出的，对任意模 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)

待办事项：证明 (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。）

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

于 2009-10-05T23:21:27.190 回答
5

这可以通过重复平方的方法来完成，即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]

第二个只是第一个的推论。

在您的情况下，b=i^m可以及时计算O(log 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^n-b^n)/(a-b)模数m，因为矩阵求幂给出：

[[b,1],[0,a]]^n == [[b^n,a^(n-1)+a^(n-2)b+...+ab^(n-2)+b^(n-1)],[0,a^n]]
于 2017-02-03T21:37:54.940 回答
2

基于@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运算。这是必不可少的，如果您将此算法转录为整数字长有限的语言。

于 2015-12-25T11:41:32.180 回答