CODEWITHSUNDEEP
pythonalgorithmmathlanguage-agnostic
dread_cat_pirate
dread_cat_pirate

Reputation: 171

Finding remainder mod involving exponent and division involving huge numbers

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)

Upvotes: 3

Views: 1061

Answers (3)

samgak
samgak

Reputation: 24417

(((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 computed by calculating this:

((a ** n) % ((a-1)*p))

and then adjusting for the -1 afterwards.

Raise a by to the nth power and mod by ((a-1)*p). This can be done using the Python pow() function. Then adjust for the -1 and divide by a-1.

Using the pow() function and passing a modulo value is faster than computing the full exponent and then taking the modulo, because the modulo can be applied to the partial products at each stage of the calculation, which stops the value from getting too large (106 to the power of 106 has 6 million decimal digits, with a modulo applied at each step the values never have to grow larger than the size of the modulo - 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.

Upvotes: 5

rici
rici

Reputation: 241671

(an−1) ⁄ (a−1) is the sum from i = 0 to n−1 of ai.

Computing the latter mod p is straightforward, based on the following:

let F(a, n) be Σ(i=0..n-1){ai} if n > 0, otherwise 0.

Now:

  1. F(a,n) = a×F(a,n−1) + 1

  2. F(a,2n) = (a+1)×F(a2,n)

The second identity is the divide-and-conquer recursion.

Since both of these only involve addition and multiplication, we can compute them mod p without needing an integer type larger than a×p by distributing the modulus operation. (See code below.)

With just the first recursion, we can code 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-conquer recursion as well, we arrive at something not 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 returns in less than a second with n == 106. The second one returns instantly, even with n as large as 109.

Upvotes: 3

Niklas B.
Niklas B.

Reputation: 95278

You can multiply by the modular inverse of p - 1. Due to Fermat's little theorem, you have xp-2 · x ≡ xp-1 ≡ 1 (mod p) for all 0 < x < p, so you don't even need extended Euclid to compute the inverse, just the pow function from standard Python:

(pow(a, n, p) - 1) * pow(a - 1, p - 2, p) % p

The algorithm has time complexity 𝒪(log p) because square-and-multiply is used.

Upvotes: 2

Related Questions