Modular Exponentiation for High Numbers in C++

Modular Exponentiation for high numbers in C++

Exponentiation by squaring still "works" for modulo exponentiation. Your problem isn't that 2 ^ 168277 is an exceptionally large number, it's that one of your intermediate results is a fairly large number (bigger than 2^32), because 673109 is bigger than 2^16.

So I think the following will do. It's possible I've missed a detail, but the basic idea works, and this is how "real" crypto code might do large mod-exponentiation (although not with 32 and 64 bit numbers, rather with bignums that never have to get bigger than 2 * log (modulus)):

Start with exponentiation by squaring, as you have.
Perform the actual squaring in a 64-bit unsigned integer.
Reduce modulo 673109 at each step to get back within the 32-bit range, as you do.

Obviously that's a bit awkward if your C++ implementation doesn't have a 64 bit integer, although you can always fake one.

There's an example on slide 22 here: http://www.cs.princeton.edu/courses/archive/spr05/cos126/lectures/22.pdf, although it uses very small numbers (less than 2^16), so it may not illustrate anything you don't already know.

Your other example, 4000111222 ^ 3 % 1608 would work in your current code if you just reduce 4000111222 modulo 1608 before you start. 1608 is small enough that you can safely multiply any two mod-1608 numbers in a 32 bit int.

Power for big numbers modulo m in C

At GeeksForGeeks is this function:

// Returns (a * b) % mod
long long moduloMultiplication(long long a,
                               long long b,
                               long long mod)
{
    long long res = 0;  // Initialize result

    // Update a if it is more than
    // or equal to mod
    a %= mod;

    while (b)
    {
        // If b is odd, add a with result
        if (b & 1)
            res = (res + a) % mod;

        // Here we assume that doing 2*a
        // doesn't cause overflow
        a = (2 * a) % mod;

        b >>= 1;  // b = b / 2
    }

    return res;
}

Use it instead of

x = (x*x) % p;

I.e.,

x = moduloMultiplication(x, x, p);

and instead of

res = (res*x) % p

I.e.,

res = moduloMultiplication(res, x, p);

Modular Exponentiation with big numbers

The problem was in integer overflowing in multiplication.

Arbitrary-precision arithmetic Explanation

It's all a matter of adequate storage and algorithms to treat numbers as smaller parts. Let's assume you have a compiler in which an int can only be 0 through 99 and you want to handle numbers up to 999999 (we'll only worry about positive numbers here to keep it simple).

You do that by giving each number three ints and using the same rules you (should have) learned back in primary school for addition, subtraction and the other basic operations.

In an arbitrary precision library, there's no fixed limit on the number of base types used to represent our numbers, just whatever memory can hold.

Addition for example: 123456 + 78:

Working from the least significant end:

initial carry = 0.
56 + 78 + 0 carry = 134 = 34 with 1 carry
34 + 00 + 1 carry = 35 = 35 with 0 carry
12 + 00 + 0 carry = 12 = 12 with 0 carry

This is, in fact, how addition generally works at the bit level inside your CPU.

Subtraction is similar (using subtraction of the base type and borrow instead of carry), multiplication can be done with repeated additions (very slow) or cross-products (faster) and division is trickier but can be done by shifting and subtraction of the numbers involved (the long division you would have learned as a kid).

I've actually written libraries to do this sort of stuff using the maximum powers of ten that can be fit into an integer when squared (to prevent overflow when multiplying two ints together, such as a 16-bit int being limited to 0 through 99 to generate 9,801 (<32,768) when squared, or 32-bit int using 0 through 9,999 to generate 99,980,001 (<2,147,483,648)) which greatly eased the algorithms.

Some tricks to watch out for.

1/ When adding or multiplying numbers, pre-allocate the maximum space needed then reduce later if you find it's too much. For example, adding two 100-"digit" (where digit is an int) numbers will never give you more than 101 digits. Multiply a 12-digit number by a 3 digit number will never generate more than 15 digits (add the digit counts).

2/ For added speed, normalise (reduce the storage required for) the numbers only if absolutely necessary - my library had this as a separate call so the user can decide between speed and storage concerns.

3/ Addition of a positive and negative number is subtraction, and subtracting a negative number is the same as adding the equivalent positive. You can save quite a bit of code by having the add and subtract methods call each other after adjusting signs.

4/ Avoid subtracting big numbers from small ones since you invariably end up with numbers like:

         10
         11-
-- -- -- --
99 99 99 99 (and you still have a borrow).

Instead, subtract 10 from 11, then negate it:

11
10-
--
 1 (then negate to get -1).

Here are the comments (turned into text) from one of the libraries I had to do this for. The code itself is, unfortunately, copyrighted, but you may be able to pick out enough information to handle the four basic operations. Assume in the following that -a and -b represent negative numbers and a and b are zero or positive numbers.

For addition, if signs are different, use subtraction of the negation:

-a +  b becomes b - a
 a + -b becomes a - b

For subtraction, if signs are different, use addition of the negation:

 a - -b becomes   a + b
-a -  b becomes -(a + b)

Also special handling to ensure we're subtracting small numbers from large:

small - big becomes -(big - small)

Multiplication uses entry-level math as follows:

475(a) x 32(b) = 475 x (30 + 2)
               = 475 x 30 + 475 x 2
               = 4750 x 3 + 475 x 2
               = 4750 + 4750 + 4750 + 475 + 475

The way in which this is achieved involves extracting each of the digits of 32 one at a time (backwards) then using add to calculate a value to be added to the result (initially zero).

ShiftLeft and ShiftRight operations are used to quickly multiply or divide a LongInt by the wrap value (10 for "real" math). In the example above, we add 475 to zero 2 times (the last digit of 32) to get 950 (result = 0 + 950 = 950).

Then we left shift 475 to get 4750 and right shift 32 to get 3. Add 4750 to zero 3 times to get 14250 then add to result of 950 to get 15200.

Left shift 4750 to get 47500, right shift 3 to get 0. Since the right shifted 32 is now zero, we're finished and, in fact 475 x 32 does equal 15200.

Division is also tricky but based on early arithmetic (the "gazinta" method for "goes into"). Consider the following long division for 12345 / 27:

       457
   +-------
27 | 12345    27 is larger than 1 or 12 so we first use 123.
     108      27 goes into 123 4 times, 4 x 27 = 108, 123 - 108 = 15.
     ---
      154     Bring down 4.
      135     27 goes into 154 5 times, 5 x 27 = 135, 154 - 135 = 19.
      ---
       195    Bring down 5.
       189    27 goes into 195 7 times, 7 x 27 = 189, 195 - 189 = 6.
       ---
         6    Nothing more to bring down, so stop.

Therefore 12345 / 27 is 457 with remainder 6. Verify:

  457 x 27 + 6
= 12339    + 6
= 12345

This is implemented by using a draw-down variable (initially zero) to bring down the segments of 12345 one at a time until it's greater or equal to 27.

Then we simply subtract 27 from that until we get below 27 - the number of subtractions is the segment added to the top line.

When there are no more segments to bring down, we have our result.

Keep in mind these are pretty basic algorithms. There are far better ways to do complex arithmetic if your numbers are going to be particularly large. You can look into something like GNU Multiple Precision Arithmetic Library - it's substantially better and faster than my own libraries.

It does have the rather unfortunate misfeature in that it will simply exit if it runs out of memory (a rather fatal flaw for a general purpose library in my opinion) but, if you can look past that, it's pretty good at what it does.

If you cannot use it for licensing reasons (or because you don't want your application just exiting for no apparent reason), you could at least get the algorithms from there for integrating into your own code.

I've also found that the bods over at MPIR (a fork of GMP) are more amenable to discussions on potential changes - they seem a more developer-friendly bunch.

Big Integer Modular Exponentiation

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);
}

Issue with Modular Exponentiation C++

First of all, some general advice:

It's better not to use strings when working with integers, as operations with strings are much slower and might become a bottleneck for performance. It's also less clear what is actually being done when strings are involved.
You shouldn't use std::pow with integers, because it operates on floating-point numbers and loses precision.

For the main question, as a workaround, you can use this O(log^2(n)) solution, which should work for arguments up to 63 bits (since it only ever uses addition and multiplication by 2). Note how all that string magic is unnecessary if you just iterate over the bits in small-to-large order:

#include <cstdint>

uint64_t modular_mul(uint64_t a, uint64_t b, uint64_t mod) {
    uint64_t result = 0;
    for (uint64_t current_term = a; b; b >>= 1) {
        if (b & 1) {
            result = (result + current_term) % mod;
        }
        current_term = 2 * current_term % mod;
    }
    return result;
}

uint64_t modular_pow(uint64_t base, uint64_t exp, uint64_t mod) {
    uint64_t result = 1;
    for (uint64_t current_factor = base; exp; exp >>= 1) {
        if (exp & 1) {
            result = modular_mul(result, current_factor, mod);
        }
        current_factor = modular_mul(current_factor, current_factor, mod);
    }
    return result;
}

Also, in gcc a (non-standard) __uint128_t is available for some targets. (which can be used to replace modular_mul with normal multiplication)

Modular exponentiation overflows when using multiple of two uint_32 numbers

The problem is that, due to the modulus being > 2⁶³, the a + sum step in your modmult routine might overflow, losing a bit. The same might happen with 2*a.

One way to fix the problem would be to add a modadd routine:

uint64_t modadd(uint_64_t a, uint64_t b, uint64_t mod) {
    if (a >= mod) a %= mod;    // precondition -- might not be needed if the caller can guarentee it.
    if (b >= mod) b %= mod;    // precondition -- might not be needed if the caller can guarentee it

    a += b;
    if (a >= mod || a < b) a -= mod;
    return a;
}

Then, use modadd(sum, a, mod) and modadd(a, a, mod) in your modmult routine.

Modular Exponentiation for High Numbers in C++