Power for big numbers modulo m in C

Question

I am using the following function to compute powers of big numbers modulo m, where m is any integer,i.e. (a^b)%m

long long power(long long x, long long y, long long p)
{
    long long res = 1;      // Initialize result

    x = x % p;  // Update x if it is more than or
                // equal to p

    while (y > 0)
    {
        // If y is odd, multiply x with result
        if (y & 1)
            res = (res*x) % p;

        // y must be even now
        y = y>>1; // y = y/2
        x = (x*x) % p;
    }
    return res;
}

But, for some numbers even this function is not working. For example, if I call

 power(1000000000000,9897,52718071807);

I get a negative number as the output. It happens because of the following reason: There is a line in the power function:

  x = (x*x) % p;

When x is big, let's say x=46175307575, the value stored in x after executing x=(x * x)%p becomes negative. I don't understand why it happens. Even if the value of (x * x) crosses the upper range of long long int, I am not storing its value anywhere, I am just storing (x*x)%p whose value should lie between 0 to p. Also, since p doesn't crosses the long long range, how does x cross it? Please tell me why this problem occurs and how to solve this.

`x*x` is still evaluated as `long long`, you have an signed integer overflow which leads to undefined behavior. — Osiris, Aug 05 '18 at 14:45
Read about [*two's complement*](https://en.wikipedia.org/wiki/Two's_complement) which is the most common way to represent negative numbers. And perhaps use `unsigned` numbers to see if it helps? If not then you need a *bignum* library. — Some programmer dude, Aug 05 '18 at 14:45
You might want to use some [bignum](https://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic) library like [GMPlib](http://gmplib.org/) — Basile Starynkevitch, Aug 05 '18 at 14:46
You should check out this post by chux [Modular exponentiation without range restriction](https://codereview.stackexchange.com/q/187257/110790) — Joseph Wood, Aug 05 '18 at 16:16

score 3 · Accepted Answer · answered Aug 05 '18 at 14:55

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

`res + a` and `2 * a` both can still overflow. This approach reduces where OF occurs, yet not eliminated it. — chux - Reinstate Monica, Aug 06 '18 at 05:23

chux - Reinstate Monica · Answer 2 · 2018-08-06T06:02:19.243

Welcomed to signed integer overflow and undefined behavior (UB).

I am just storing (x*x)%p whose value should lie between 0 to p.

This is incorrect. x*x may overflow long long math and the result is UB. @Osiris. Sample UB includes a product that is negative with positive operands..

some_negative_value % some_positive_p results in a negative value. - see ref. This is outside the range [0...p).

The solution is to not overflow signed integer math.

A simple 1st step is to use unsigned integer math.

A full range solution without overflow problem is here Modular exponentiation without range restriction

Note OP's code also fails a corner case: power(some_x, 0, 1) as it returns 1 when 0 is expected.

// Fix
// long long res = 1;
long long res = 1%p;
// or 
long long res = p != 1;

GolamMazid Sajib · Answer 3 · 2018-08-06T06:03:58.973

0

If Mod == (2^63-1) then this solution will not work.

Solution: Mod <= 2^62

(p-1)*(p-1) > 2^63. so there will be overflow. you need to implement multiplication with modulo.

Try with this:

long long multiply(long long a,long long b,long long m){
if(b == 0){
    return 0;
}
if(b==1){
    return a%m;
}
if(b&1){
    return ((a%m)+multiply(a,b-1,m))%m;
}
long long x = multiply(a,b>>1,m);
return multiply(x,2,m);
}

long long bigmod(long long a,long long b, long long m){
if(b == 0){
    return 1;
}
if(b == 1){
    return a%m;
}
if(b & 1){
    return multiply(a%m,bigmod(a,b-1,m),m);
}
long long x = bigmod(a,b>>1,m);
return multiply(x,x,m);
}

edited Aug 06 '18 at 06:03

answered Aug 05 '18 at 15:01

GolamMazid Sajib

8,698
6
21
39

`(a%m)+multiply(a,b-1,m)` and `x<<1` can still overflow. This approach reduces where OF occurs, yet not eliminated it. – chux - Reinstate Monica Aug 06 '18 at 05:25
how (x<<1) be overflow? – GolamMazid Sajib Aug 06 '18 at 05:36
(a%m)+multiply(a,b-1,m) here used `+` not `*` – GolamMazid Sajib Aug 06 '18 at 05:36
multiply(a,b-1,m) return %m res – GolamMazid Sajib Aug 06 '18 at 05:37
how (x<<1) be overflow? --> `multiply(LLONG_MAX-1, 2, LLONG_MAX);` – chux - Reinstate Monica Aug 06 '18 at 05:46
`(a%m)+multiply(a,b-1,m)` can overflow too. Try `multiply(LLONG_MAX-1, 2, LLONG_MAX)`. – chux - Reinstate Monica Aug 06 '18 at 05:55
A solution, see [`addmodmax()`](https://codereview.stackexchange.com/questions/187257/modular-exponentiation-without-range-restriction) is mentioned in this [answer](https://stackoverflow.com/a/51701121/2410359). – chux - Reinstate Monica Aug 06 '18 at 06:03

score 0 · Answer 4 · answered Aug 05 '18 at 21:56

Apart from solution mentioned by @Doug Currie, you can also use 128 bits data type __int128.

long long pow(long long a, long long b, long long mod)
{
    __int128 res = 1;
    while(b > 0)
    {
        if(b&1)
        {
            res = (res*a);
            res = res%mod;
        }
        b = b>>1;
        a = ((__int128)a*a)%mod;
    }
    return res;
}

Power for big numbers modulo m in C

4 Answers4

Related