Finding residual modes involving exponential and huge fission
I need a fast algorithm to evaluate the following
((a^n-1)/(a-1)) % p
Both a
and are n
almost equal, but less than 10 ^ 6, and p
is a fixed prime (say p=1000003
). I need to compute it within 1 second. I am using python. Wolfram Mathematica calculates it instantly . It takes 35.2170000076 seconds with the following code
print (((10**6)**(10**6)-1)/((10**6)-1))%1000003
If this denominator was a-1
not present, I could group cardinalities in a lower order and use a ratio a*b (mod c) = (a (mod c) * b (mod c)) (mod c)
, but the denominator is present.
How do you evaluate this using a fast algorithm? No numpy / scipy.
UPDATE :: Here is the last 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)
source to share
(((a ** n) - 1) / (a-1))% p
can be rewritten as
(((a ** n) - 1)% ((a-1) * p)) / (a-1)
This part:
(((a ** n) - 1)% ((a-1) * p))
can be calculated by calculating this:
((a ** n)% ((a-1) * p))
and then adjust the value -1 after.
Raise a to the nth power and mod to ((a-1) * p). This can be done using the Python pow () function. Then adjust to -1 and divide by a-1.
Using pow () and passing a value modulo is faster than calculating the full exponent and then taking modulo, since modulus can be applied to partial products at each step of the calculation, which stops the value from getting too large (10 6 to the power of 10 6 has 6 million decimal digits, with the modulo applied at each step, the values โโshould never increase more than the modulo size - about 13 digits in this example).
code:
def exp_div_mod(a, n, p):
m = p * (a - 1)
return ((pow(a, n, m) - 1) % m) // (a - 1);
print exp_div_mod((10**6), (10**6), 1000003)
output:
444446
Note: this approach only works if a, n and p are integers.
source to share
(a n & minus; 1) & frasl; (a & minus; 1) is the sum of i = 0 to n & minus 1 of i .
Calculating the latter mod p is straightforward based on the following:
let F (a, n) - & Sigma; (i = 0..n-1) {a i } if n> 0, otherwise 0.
Now:
-
F (a, n) = a x; F (a, n & minus; 1) + 1
-
F (a, 2 n ) = (a + 1) x; F (a 2 n)
The second identity is divide-and-win recursion.
Since both of them only deal with addition and multiplication, we can compute them mod p without requiring an integer type greater than a x; p, distributing the unit operation. (See code below.)
With only the first recursion, we can program an iterative solution:
def sum_of_powers(a, n, p):
sum = 0
for i in range(n): sum = (a * sum + 1) % p
return sum
Using the divide and win recursion, we end up with something much more complicated:
def sum_of_powers(a, n, p):
if n % 2 == 1:
return (a * sum_of_powers(a, n-1, p) + 1) % p
elif n > 0:
return ((a + 1) * sum_of_powers(a * a % p, n // 2, p)) % p
else:
return 0
The first solution comes back in less than a second with n == 10 6 . The second returns instantly, even if n is 10 9 .
source to share
You can multiply by modular inverse p - 1
. Because of Fermat's little theorem, you have x p-2 x โก x p-1 โก 1 (mod p) for all 0 <x <p, so you don't even need the extended Euclid to compute the inverse function, just pow
from standard Python :
(pow(a, n, p) - 1) * pow(a - 1, p - 2, p) % p
The algorithm has a time complexity of ๐ช (log p) since it uses square multiplication .
source to share