Calculating sum of geometric series (mod m)

I have a series

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

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

I want to find the sum. I cannot apply the Geometric Progression (GP) formula because then result will have denominator and then I will have to find modular inverse which may not exist (if the denominator and m are not coprime).

So I made an alternate algorithm making an assumption that these powers will make a cycle of length much smaller than k (because it is a modular equation and so I would obtain something like 2,7,9,1,2,7,9,1....) and that cycle will repeat in the above series. So instead of iterating from 0 to k, I would just find the sum of numbers in a cycle and then ca开发者_JS百科lculate the number of cycles in the above series and multiply them. So I first found i^m (mod m) and then multiplied this number again and again taking modulo at each step until I reached the first element again.

But when I actually coded the algorithm, for some values of i, I got cycles which were of very large size. And hence took a large amount of time before terminating and hence my assumption is incorrect.

So is there any other pattern we can find out? (Basically I don't want to iterate over k.) So please give me an idea of an efficient algorithm to find the sum.

This is the algorithm for a similar problem I encountered

You probably know that one can calculate the power of a number in logarithmic time. You can also do so for calculating the sum of the geometric series. Since it holds that

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

you can recursively calculate the geometric series on the right hand to get the result.

This way you do not need division, so you can take the remainder of the sum (and of intermediate results) modulo any number you want.

As you've noted, doing the calculation for an arbitrary modulus m is difficult because many values might not have a multiplicative inverse mod m. However, if you can solve it for a carefully selected set of alternate moduli, you can combine them to obtain a solution mod m.

Factor m into p_1, p_2, p_3 ... p_n such that each p_i is a power of a distinct prime

Since each p is a distinct prime power, they are pairwise coprime. If we can calculate the sum of the series with respect to each modulus p_i, we can use the Chinese Remainder Theorem to reassemble them into a solution mod m.

For each prime power modulus, there are two trivial special cases:

If i^m is congruent to 0 mod p_i, the sum is trivially 0.

If i^m is congruent to 1 mod p_i, then the sum is congruent to k mod p_i.

For other values, one can apply the usual formula for the sum of a geometric sequence:

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

TODO: Prove that (i^m - 1) is coprime to p_i or find an alternate solution for when they have a nontrivial GCD. Hopefully the fact that p_i is a prime power and also a divisor of m will be of some use... If p_i is a divisor of i. the condition holds. If p_i is prime (as opposed to a prime power), then either the special case i^m = 1 applies, or (i^m - 1) has a multiplicative inverse.

If the geometric sum formula isn't usable for some p_i, you could rearrange the calculation so you only need to iterate from 1 to p_i instead of 1 to k, taking advantage of the fact that the terms repeat with a period no longer than p_i.

(Since your series doesn't contain a j=0 term, the value you want is actually S-1.)

This yields a set of congruences mod p_i, which satisfy the requirements of the CRT. The procedure for combining them into a solution mod m is described in the above link, so I won't repeat it here.

This can be done via the method of repeated squaring, which is O(log(k)) time, or O(log(k)log(m)) time, if you consider m a variable.

In general, a[n]=1+b+b^2+... b^(n-1) mod m can be computed by noting that:

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

The second just being the corollary for the first.

In your case, b=i^m can be computed in O(log m) time.

The following Python code implements this:

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

This bit of magic has a mathematical reason - the operation on pairs defined as

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

is associative, and the rule 1 basically means that:

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

Repeated squaring always works when operations are associative. In this case, the @ operator is O(log(m)) time, so repeated squaring takes O(log(n)log(m)).

One way to look at this is that the matrix exponentiation:

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

You can use a similar method to compute (a^n-b^n)/(a-b) modulo m because matrix exponentiation gives:

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

Based on the approach of @braindoper a complete algorithm which calculates

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

looks like this in 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]]
]

Parameters:

• a is the "ratio" of the series. It can be any integer (including zero and negative values).
• n is the highest exponent of the series. Allowed are integers >= 0.
• mis the integer modulus != 0

Note: The algorithm performs a Mod operation after every arithmetic operation. This is essential, if you transcribe this algorithm to a language with a limited word length for integers.