I need a fast algorithm to evaluate the following
((a^n-1)/(a-1)) % p
Both a
and n
are nearly equal but less to 10^6 and p
is a fixed prime number (let's say p=1000003
). I need to compute it under 1 second. I am using python. Wolfram Mathematica computes it instantly. It takes 35.2170000076 seconds with following code
print (((10**6)**(10**6)-1)/((10**6)-1))%1000003
If that denominator a-1
were not present, I could group the powers into smaller order and use the relation a*b (mod c) = (a (mod c) * b (mod c)) (mod c)
but denominator is present.
How to evaluate this with a fast algorithm? No numpy/scipy are available.
UPDATE:: Here is the final code I came up with
def exp_div_mod(a, n, p):
r = pow(a, n, p*(a-1)) - 1
r = r - 1 if r == -1 else r
return r/(a-1)
1/(a-1) mod p
; then apply the last identity you wrote. – hiro protagonist Jul 4 '15 at 13:14gmpy.divm
"returns x such that b*x==a modulo m, or else raises a ZeroDivisionError exception if no such value x exists". Witha
,n
&p
as in the OP,(pow(a,n,p)-1)*gmpy.divm(1,a-1,p) % p
returns 444446. – PM 2Ring Jul 4 '15 at 13:26r = r - 1 if r == -1 else r
in your update is wrong. – PM 2Ring Jul 7 '15 at 3:15