CODEWITHSUNDEEP
calgorithmfunctionlong-long
Klasen
Klasen

Reputation: 159

Power for big numbers modulo m in C

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.

Upvotes: 1

Views: 594

Answers (5)

Michel
Michel

Reputation: 97

You can perform a power modulo in C99 :

typedef unsigned long long int ulong;

I'm using this utility function :

ulong mul_mod(ulong a, ulong b, const ulong mod) {
    ulong res = 0, c; // return (a * b) % mod, avoiding overflow errors while doing modular multiplication.
    for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (c = b) >= mod - b ? c -= mod : 0, b += c);
    return res % mod;
}

My pow_mod function is designed to compute (n ^ exp) % mod efficiently using the method of exponentiation by squaring, to handle large exponents in modular arithmetic :

ulong pow_mod(ulong n, ulong exp, const ulong mod) {
    ulong res = 1; // return (n ^ exp) % mod
    for (n %= mod; exp; exp & 1 ? res = mul_mod(res, n, mod) : 0, n = mul_mod(n, n, mod), exp >>= 1);
    return res;
}

The function reduce the time complexity to O(log exp) (it's faster than the naive approach of multiplying n by itself exp times), and should perform well in scenarios requiring modular exponentiation, like cryptographic algorithms or large number computations under modulus.

Upvotes: 0

Golam Mazid Sajib
Golam Mazid Sajib

Reputation: 9447

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

Upvotes: 0

chux
chux

Reputation: 154075

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;

Upvotes: 1

Shashwat Kumar
Shashwat Kumar

Reputation: 5297

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

Upvotes: 0

Doug Currie
Doug Currie

Reputation: 41200

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

Upvotes: 3

Related Questions