Modular Exponentiation doesn't work for large numbers in C

Question

I don't know why it doesn't work well in large numbers since I've maximized the use of unsigned long long. Please help me. Or if this code is a total mess please do suggest a better and more efficient algorithm. This is a homework and we are only allowed to use
#include <stdio.h> and #include <stdlib.h>. I hope you could help me.

unsigned long long modularexp(unsigned long long a,unsigned long long b,unsigned long long mod)
{
    //a^b % mod
    unsigned long long product,pseq;
    product=1;
    pseq=a%mod;
    while(b>0)
    {
        if(b&1)
            product=(product*pseq)%mod;
        pseq=(pseq*pseq)%mod;
        b>>=1;
    }
    return product;
}

Answer 1

Here is a test program that shows where the overflow occurs:

#include <stdio.h>
#include <stdlib.h>

static
unsigned long long modularexp(unsigned long long a,
                              unsigned long long b,
                              unsigned long long mod)
{
    unsigned long long product, pseq;
    product = 1;
    pseq = a % mod;
    while (b > 0) {
        if (b & 1) {
            if (pseq != 0 && product * pseq / pseq != product) {
                printf("multiplication overflow at b=%llu:"
                       " product=%llu, pseq=%llu\n",
                       b, product, pseq);
            }
            product = (product * pseq) % mod;
        }
        if (pseq != 0 && pseq * pseq / pseq != pseq) {
            printf("multiplication overflow at b=%llu:"
                   " pseq=%llu\n",
                   b, pseq);
        }
        pseq = (pseq * pseq) % mod;
        b >>= 1;
    }
    return product;
}

int main(int argc, char *argv[]) {
    unsigned long long a, b, mod, res;
    if (argc > 3) {
        a = strtoull(argv[1], NULL, 0);
        b = strtoull(argv[2], NULL, 0);
        mod = strtoull(argv[3], NULL, 0);
        res = modularexp(a, b, mod);
        printf("%llu ** %llu %% %llu = %llu\n", a, b, mod, res);
    }
}

Running it with your values produces this:

~/dev/stackoverflow > ./t47 1727 62483 7491569251
multiplication overflow at b=3905: product=5150827583, pseq=5846688665
multiplication overflow at b=3905: pseq=5846688665
multiplication overflow at b=488: pseq=6744552923
multiplication overflow at b=244: pseq=4362729489
multiplication overflow at b=61: pseq=6056813204
multiplication overflow at b=30: pseq=6747710262
multiplication overflow at b=15: pseq=5762416609
multiplication overflow at b=7: pseq=5992947929
multiplication overflow at b=3: pseq=6028562833
multiplication overflow at b=1: product=7015761607, pseq=7047166445
multiplication overflow at b=1: pseq=7047166445
1727 ** 62483 % 7491569251 = 4674220559

You have to protect the computations with a function like that provided by M Oehm.

Answer 2

Your code is good pre se, but it works reliably only if the modulus is less than 2^32. In your case, the modulus excees 2^32 and the following:

a * b < ulong.max 

is not necessarily true, even if a < mod and b < mod.

You can fix this by adding an overflow check to all multiplications. If the multiplication a * b overflows, subdivide a = a1 + a2 and calculate a1 * b + a2 * b.

The code below does this by introducing an additional function for modular multiplication.

#include <stdio.h>

typedef unsigned long long ulong;

#define ULONG_MAX ((ulong) -1)

ulong modmul(ulong a, ulong b, ulong mod)
{
    if (a == 0) return 0;

    if (ULONG_MAX / a > b) {
        return (a * b) % mod;
    } else {
        ulong a1 = a / 2;
        ulong a2 = a - a1;

        return (modmul(a1, b, mod) + modmul(a2, b, mod)) % mod;
    }  
}

ulong modexp(ulong a, ulong b, ulong mod)
{
    ulong product = 1;

    a = a % mod;

    while (b > 0) {
        if (b & 1) {
            product = modmul(product, a, mod);
        }
        a = modmul(a, a, mod);
        b >>= 1;
    }

    return product;
}

int main()
{
    ulong a = modexp(1727, 62483, 7491569251LL);
    ulong expect = 5500747491LL;

    printf("%llu\n%llu\n", a, expect);

    return 0;
}