Big Integer Modular Exponentiation

Question

How to calculate (x^y) mod z with 1 <= x, y <= 10¹⁰⁰⁰ and z any positive integer 1 <= z < 2³¹ ?

What I have done so far is: scan x and y as a string, get the modulo, then calculate (x^y) mod z.

I know this is wrong because (x^y) mod z is not equal to ((x mod z)^{(y mod z)}) mod z. Then how do I solve this?

Edit: sorry I made the bottom constraint of x and y so high when creating the question. I just want to make other focus on the big integer problem, not the modular exponentiation :).

#define MOD z

long long power (long long k, long long n) {
    if (n == 1) return k;
    else {
        long long p = power (k, n/2);
        if (n % 2 == 0) return (p * p) % MOD;
        else return (((p * p) % MOD) * k) % MOD;
    }
}

long long convert (char *n) {
    long long number = 0;
    int ln = strlen (n);
    
    for (int x = 0; x < ln; x++) {
        number = number * 10;
        number = (number + (n[x] - '0')) % MOD;
    }
    
    return number % MOD;
}

int main () {
    char s_x[1111], s_y[1111];
    scanf ("%s %s", s_x, s_y);
    
    long long x, y, r;
    x = convert (s_x);
    y = convert (s_y);
    r = power (x, y);
        
    printf ("%lld\n", r);
}

Answer 1

Since modular exponentiaion is used so much, there are libraries for it. The following is an example that reads a, b, and c and outputs a^b mod c using GMP.

#include <stdio.h>
#include <gmp.h>

int main(void)
{
  mpz_t a, b, c, d;
  mpz_inits (a, b, c, d, NULL);
  printf ("a: ");
  mpz_inp_str (a, stdin, 10);
  printf ("b: ");
  mpz_inp_str (b, stdin, 10);
  printf ("c: ");
  mpz_inp_str (c, stdin, 10);
  mpz_powm (d, a, b, c); // compute d = a ^ b mod c
  gmp_printf ("a ^ b mod c = %Zd\n", d);
  return 0;
}

Compile it with -lgmp.

By the way, a^b ≡ a^{b mod Φ(c)} (mod c), where Φ is Euler's totient function

Answer 2

First, I assume that z is rather small (as in, fits into long). Also observe that

(x ^ y) % z = ((x % z) ^ y) % z

So it is fine to convert x the way you do, the only problem is y. Conveniently, you only do two things with y -- you divide it by two, and you check the remainder after division by two. Both of those things are trivial if you represent y as an array. First, for simplicity reverse the y, so that the least significant digit goes first, and also store digits, not digit characters in the array (as in, store 5, not '5'). You might also consider storing more than just one digit in each element, but this only improves it by a constant.

Now to check the remainder just check if the first element of the array is divisible by two (the number is even if its least significant digit is even). To divide by two, do something along the lines of:

for (int i = 0; i < y_len; ++ i) {
    if (i && y[i] % 2) y[i - 1] += 5;
    y[i] /= 2;
}
if (y_len && y[y_len - 1] == 0) -- y_len;

Plug this into your power routine, and it will work just fine. Note that your power method is logarithmic in y, so the fact that y can be up to 10^1000 doesn't make it unmanageably slow.

Answer 3

I guess you are trying to build a Diffie-Hellman Key Exchange algorithm. Try importing the OpenSSL library and then use it's BN_mod_exp() function.

BN_mod_exp() computes a to the p-th power modulo m (r=a^p % m). This function uses less time and space than BN_exp().

Source: https://www.openssl.org/docs/manmaster/crypto/BN_add.html

Answer 4

Thanks to @v7d8dpo4 explanation for Euler's Totient Function. I edited my code as of following:

#define MOD z

long long power (long long k, long long n) {
    if (n == 1) return k;
    else {
        long long p = power (k, n/2);
        if (n % 2 == 0) return (p * p) % MOD;
        else return (((p * p) % MOD) * k) % MOD;
    }
}

long long convert (char *n, int mod) {
    long long number = 0;
    int ln = strlen (n);

    for (int x = 0; x < ln; x++) {
        number = number * 10;
        number = (number + (n[x] - '0')) % mod;
    }

    return number % mod;
}

int main () {
    char s_x[1111], s_y[1111];
    scanf ("%s %s", s_x, s_y);

    long long x, y, r;
    x = convert (s_x, MOD);
    y = convert (s_y, totient (MOD)); // totient (x) is Euler's Totient Function of x
    r = power (x, y);

    printf ("%lld\n", r);
}