Modular exponentiation function returns the wrong int value

Question

I got a modular exponentiation function in C which looks like this.

int modexp(int m, int e, int n)
{
    printf("%i\n", m);
    printf("%i\n", e);
    printf("%i\n", n);
    printf("\n");

    if(e == 0)
    {
        return 1;
    }
    if(e%2 == 1)
    {
        return modexp(m, e-1, n) * m%n;
    }
    else
    {
        int modTemp = modexp(m, (int)(e/2), n);
        modTemp = modTemp * modTemp;
        return modTemp % n;
    }
}

which I am calling in my main() function like this

int p = 3511;
printf("%i\n", modexp(2, p-1, p*p));

When printing the values for m, e and n I get the correct recursion values up to e = 0 in the end. This is when the function should return 1. It definitely returns at that position in the code, however instead of the expected integer 1, I get -6593454 and I have no idea why.

The full code can be found here: https://gist.github.com/anonymous/7024ac77b2432a381968

any input is highly appreciated...

Answer 1

Multipying an n-bit value with an m-bit value produces a (n+m)-bit result. That's why modexp(m, e-1, n) * m and modTemp * modTemp overflow. You need a widening multiplication to get the full 64-bit result from 32-bit values

return ((int64_t)modexp(m, e-1, n) * m) % n;
...
int64_t modTemp = modexp(m, (int)(e/2), n);
modTemp = modTemp * modTemp;
return modTemp % n;

Answer 2

Your numbers are overflowing. They are most likely 32 or 64-bit signed integers here, which means the maximum size can be 2³¹=2147483648 or 2⁶³ = 9223372036854775807. Just look at the size of the modexp(...) multiplied by m%n. That's quite a large number :)

The same goes for modTemp * modTemp

Answer 3

1 is indeed returned when e == 0, but this return value is a term in the calculations made on previous levels of recursion. If you modify your code like this:

#include <stdio.h>
int modexp(int m, int e, int n)
{
    printf("%i\n", m);
    printf("%i\n", e);
    printf("%i\n", n);
    printf("\n");

    if(e == 0)
    {
        printf("returning 1\n");
        return 1;
    }
    if(e%2 == 1)
    {
        int tmp=modexp(m, e-1, n) * m%n;
        printf("returning %d\n",tmp);
        return tmp; 
    }
    else
    {
        int modTemp = modexp(m, (int)(e/2), n);
        modTemp = modTemp * modTemp;
        modTemp= modTemp % n;
        printf("returning %d\n",modTemp);
        return modTemp;;
    }
}


int main(){
    int p = 3511;
    printf("%d\n", modexp(2, p-1, p*p));
    return 0;
}

You will get this output at the end:

returning 1
returning 2
returning 4
returning 8
returning 64
returning 4096
returning 8192
returning 5473259
returning 10946518
returning 2217782
returning 5855618
returning 11711236
returning 8916378
returning 5505635
returning -1026672
returning 975519
returning 1951038
returning 195330
returning 390660
returning -6593454
-6593454

This shows that the 1 returned for e == 0 is used for further calculations.

At some point your code overflows int, that's why you get a negative number at the end. You should probably use a longer type.

Answer 4

I agree it's an overflow of the modTemp variable. However, for clarity, I'd suggest to only change the type of modTemp from Int to long, and change m/e/n with b/e/m like the following:

#include <stdio.h>

int modexp(int b, int e, int m); // b: base; e: exp; m: modulus


int main()
{
    int p = 3511; // given number p
    printf("%i\n", modexp(2, p-1, p*p)); // last line printed is the remainder of mod exp
}


int modexp(int b, int e, int m)
{
    printf("%i\n", b);
    printf("%i\n", e);
    printf("%i\n", m);
    printf("\n");

    if(e == 0)
    {
        return 1;
    }
    if(e%2 == 1)
    {
        return modexp(b, e-1, m) * b%m;
    }
    else
    {
        long modTemp = modexp(b, (int)(e/2), m); // !!here modTemp is declared as long!!
        modTemp = modTemp * modTemp;
        return modTemp % m;
    }
}