Finding modulus with very large exponentiation and divisor in C

Question

I need a function in C that calculates a^n mod q, where the divisor q is determined to be very big (15,383,399,235,709,406,497) and the exponent n may also be as large as that.

Based on the property of modulo multiplication, that (a * b) mod n = ((a mod n) * (b mod n)) mod n, my attempt is as follows:

#include <stdio.h>

unsigned long long int modExp(unsigned long long int base, unsigned long long int expo, unsigned long long int mod)
{
    unsigned long long int out = 1;
    unsigned long long int cnt;
    for (cnt=expo; cnt>0; cnt--)
    {
        out = (out * base) % mod;
    }
    return out;
}

int main()
{
    printf("%llu", modExp(3, (189 + 50 * 223), 15383399235709406497));
    return 0;
}

As seen in the main function above, I tested my function modExp with base 3, exponent (189 + 50 * 223), and divisor 15383399235709406497. It gave the output 3915400295876975163 with some warnings that

In function ‘findxA’:
findxA.c:17:32: warning: integer constant is so large that it is unsigned [enabled by default]
     unsigned long long int h = 12036625823877237123;
                                ^
findxA.c:17:5: warning: this decimal constant is unsigned only in ISO C90 [enabled by default]
     unsigned long long int h = 12036625823877237123;
     ^
findxA.c:18:32: warning: integer constant is so large that it is unsigned [enabled by default]
     unsigned long long int q = 15383399235709406497;
                                ^
findxA.c:18:5: warning: this decimal constant is unsigned only in ISO C90 [enabled by default]
     unsigned long long int q = 15383399235709406497;
     ^
findxA.c: In function ‘main’:
findxA.c:34:48: warning: integer constant is so large that it is unsigned [enabled by default]
     printf("%llu", modExp(3, (189 + 50 * 223), 15383399235709406497));
                                                ^
findxA.c:34:5: warning: this decimal constant is unsigned only in ISO C90 [enabled by default]
     printf("%llu", modExp(3, (189 + 50 * 223), 15383399235709406497));
     ^

To verify this answer, I compared the result (given by the C function) with the output given by evaluating the expression 3^(189 + 50 * 223) `mod` 15383399235709406497, written in Haskell. This Haskell expression evaluated to a different decimal, 12349118475990906085. I think it is my C function that is wrong, since I believe Haskell does a better job at handling such large decimals.

How can I improve my C function, modExp?

Edit: I have tried the option of the first answer of this question. However, as I am trying to deal with the decimal of unsigned long long int, I have switched every input and return types from int to unsigned long long int. This resulted in a Segmentation fault.

Edit2: I have used the function described on the above link in a wrong way. It does not give Segmentation fault; but it still does not output the correct decimal.

Answer 1

To reduce the chance of overflow, you can rely on (x*y)%z == ((x%z) * (y%z)) % z.

For example (untested):

unsigned long long int modExp(unsigned long long int base, unsigned long long int expo, unsigned long long int mod)
{
    unsigned long long int baseMod = base%mod;
    unsigned long long int out = 1;
    unsigned long long int cnt;
    for (cnt=expo; cnt>0; cnt--)
    {
        out = ( (out%mod) * baseMod) % mod;
    }
    return out;
}

Also note that exponentiation can be done more efficiently as "product of squares". For example, x ** 5 == 1*(x**4) * 0*(x**2) * 1*(x**1) because 5 == 1*4 + 0*2 + 1*1 (or 5 == 101b).

For example (untested):

unsigned long long int modExp(unsigned long long int base, unsigned long long int expo, unsigned long long int mod)
{
    unsigned long long int out = 1;
    unsigned long long int temp = base%mod;

    while(exp0 > 0) {
        if( (exp0 & 1) != 0) {
            out = (out * temp) % mod;
        }
        temp = (temp*temp) % mod;
        exp0 >>= 1;
    }
    return out;
}

For large exponents this should make a huge difference in performance (e.g. if the exponent is 123456 then the loop would have 17 iterations instead of 123456 iterations).

Also; if this is for some kind of cryptography (not unlikely); then you should state that clearly because you will probably want "constant time" (to reduce the chance of timing side-channels - e.g. inferring information about the exponent from how much time modExp() takes to execute).

Finally; even with the changes, numbers up to 15,383,399,235,709,406,497 are probably too large for unsigned long long (you'd need to ensure mod * mod <= ULLONG_MAX); which means that you'll probably need to use/implement "big integer" operations (e.g. maybe a typedef struct { uint32_t digit[4]; } MY_NUMBER_TYPE thing to handle 128-bit integers, with functions for multiplication and modulo).

Answer 2

This will never work this way!

This result of the multiplication here alone

out = (out * base) % mod;

needs potentially more than the 64 bit of the underlying data type. And if an integer overflow occurs (ie, the most significant bits get cut off), the mod operation turns out wrong!

Use larger data types, or use a different approach :-)

If you like, btw, please use this test code to see that there is actually an integer overflow occuring TWICE with your inputs:

#include <stdio.h>

unsigned long long int modExp(unsigned long long int base, unsigned long long int expo, unsigned long long int mod)
{
    unsigned long long int out = 1;
    while(expo>0)
    {
        if(expo % 2 == 1)
            out = (out * base) % mod;
        expo = expo >> 1;
        if (base * base < base) printf("WARNING! INTEGER OVERFLOW!!!!\n");
        base = (base * base) % mod;
    }
    return out;
}

int main()
{
    printf("%llu", modExp(3, (189 + 50 * 223), 15383399235709406497));
    return 0;
}