Q:

Finding modulus with very large exponentiation and divisor in C

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.

A:

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).

