Calculate A*A Mod N Without Overflow

How to do arithmetic modulo another number, without overflow?

Richard Rast pointed out that Wikipedia version works only with 63-bit integers. I extended the code provided by Boiethios to work with full range of 64-bit unsigned integers.

fn mul_mod64(mut x: u64, mut y: u64, m: u64) -> u64 {
    let msb = 0x8000_0000_0000_0000;
    let mut d = 0;
    let mp2 = m >> 1;
    x %= m;
    y %= m;

    if m & msb == 0 {
        for _ in 0..64 {
            d = if d > mp2 {
                (d << 1) - m
            } else {
                d << 1
            };
            if x & msb != 0 {
                d += y;
            }
            if d >= m {
                d -= m;
            }
            x <<= 1;
        }
        d
    } else {
        for _ in 0..64 {
            d = if d > mp2 {
                d.wrapping_shl(1).wrapping_sub(m)
            } else {
                // the case d == m && x == 0 is taken care of 
                // after the end of the loop
                d << 1
            };
            if x & msb != 0 {
                let (mut d1, overflow) = d.overflowing_add(y);
                if overflow {
                    d1 = d1.wrapping_sub(m);
                }
                d = if d1 >= m { d1 - m } else { d1 };
            }
            x <<= 1;
        }
        if d >= m { d - m } else { d }
    }
}

#[test]
fn test_mul_mod64() {
    let half = 1 << 16;
    let max = std::u64::MAX;

    assert_eq!(mul_mod64(0, 0, 2), 0);
    assert_eq!(mul_mod64(1, 0, 2), 0);
    assert_eq!(mul_mod64(0, 1, 2), 0);
    assert_eq!(mul_mod64(1, 1, 2), 1);
    assert_eq!(mul_mod64(42, 1, 2), 0);
    assert_eq!(mul_mod64(1, 42, 2), 0);
    assert_eq!(mul_mod64(42, 42, 2), 0);
    assert_eq!(mul_mod64(42, 42, 42), 0);
    assert_eq!(mul_mod64(42, 42, 41), 1);
    assert_eq!(mul_mod64(1239876, 2948635, 234897), 163320);

    assert_eq!(mul_mod64(1239876, 2948635, half), 18476);
    assert_eq!(mul_mod64(half, half, half), 0);
    assert_eq!(mul_mod64(half+1, half+1, half), 1);

    assert_eq!(mul_mod64(max, max, max), 0);
    assert_eq!(mul_mod64(1239876, 2948635, max), 3655941769260);
    assert_eq!(mul_mod64(1239876, max, max), 0);
    assert_eq!(mul_mod64(1239876, max-1, max), max-1239876);
    assert_eq!(mul_mod64(max, 2948635, max), 0);
    assert_eq!(mul_mod64(max-1, 2948635, max), max-2948635);
    assert_eq!(mul_mod64(max-1, max-1, max), 1);
    assert_eq!(mul_mod64(2, max/2, max-1), 0);
}

Modulus Function to Avoid Integer Overflow in C++

a %= p
b %= p
c = p-a

if(b==c)
  sum = 0
if (b<c)
 sum = a+b
if (b>c)
 sum = b-c

EDIT: The trick is to avoid any calculation that might cause overflow, without knowing where the limit is. All we know is that the given a, b and p are below the limit -- maybe just below the limit.

After the first two steps (a%=p;b%=p;) we know a<p and b<p. We still daren't add a+b, because that sum might exceed p and break the limit^*. But we can see how much room we have left with c = p-a, which is safe because we know that c<=p and c>0. (The stated types are unsigned, but we may as well avoid negative numbers, if only because their limits are sometimes off by one from the negatives of the positive limits, in ways I can never remember.)

If b=c, then b=p-a, so a+b=p, so the sum (mod p) is zero.

If b<c, then a+b<p, so we can safely compute a+b (and need not apply the modulo).

if b>c, then it is not safe to compute a+b, but we know that the number we're looking for is a+b-p, which we can rewrite as b-(p-a), and we already have b and p-a, so we can safely perform that subtraction.

(*) _{That's right, I said "daren't". It's a perfectly good word.}

How to efficiently verify whether pow(a, b) % b == a in C (without overflow)

You can use the Exponentiation by squaring method.

The idea is the following:

• Decompose b in binary form and decompose the product
• Notice that we always use %b which is below 32768, so the result will always fit in a 32 bit number.

So the C code is:

/*
 * this function computes (num ** pow) % mod
 */
int pow_mod(int num, int pow, int mod)
{
    int res = 1

    while (pow>0)
    {
        if (pow & 1)
        {
            res = (res*num) % mod;
        }
        pow /= 2;
        num = (num*num)%mod;
    }

    return res;
}

How to compute this modulus when there is an integer overflow

You can use the identity you quoted, you just need another similar identity: (a+b) mod m = (a mod m) + (b mod m).

The goal is to multiply x*y mod m without any intermediate values exceeding the overflow limit (in this case 2^64), where x is starting less than m (if it isn't, reduce it mod m), y is possibly larger than m, and x*y can potentially overflow. We can do this if m is less than half of the overflow limit.

A solution is simple: Just perform basic multiplication bit-by-bit for x*y and do every step modulo m.

Start with x and y less than m (if either isn't, reduce it first). Write y in the form a_0 * 2^0 + a_1 * 2^1 + a_2 * 2^2 + ... , where a_n is either 0 or 1 (indicating the term is present or not). (Aka, write y in binary form.) Now we have:

x * (a_0 * 2^0 + a_1 * 2^1 + a_2 * 2^2 + ...) mod m

Distribute x over each of the terms of y:

(x * a_0 * 2^0) + (x * a_1 * 2^1) + (x * a_2 * 2^2) + ... mod m

Then use the original multiplication identity: For each term above, multiply x by 2 mod m until you reach the desired power of 2 for that term. (Since x < m and 2 * m < 2^64, then 2 * x < 2^64, so we can multiply by 2 without overflowing.) When you are done, add the result for each term mod m (you can keep a running sum as you go).

None of those operations will exceed 2^64 and thus will not overflow. This will work for any value of m less than 2^64 / 2 = 2^63 and any integers x and y less than m.

This is not necessarily the fastest way to do it, feel free to find something more efficient. For starters, the smaller m is compared to the overflow limit, the bigger the radix for the terms we can rewrite y as.

How can we compute N choose K modulus a prime number without overflow?

To compute (n choose k) % M, you can separately compute the nominator (n!) modulus M and the denominator (k!*(n - k)!) modulus M and then multiply the nominator by the denominator's modular multiplicative inverse (in M). Since M is prime, you can use Fermat's Little Theorem to calculate the multiplicative inverse.

There is a nice explanation, with sample code, on the following link (problem SuperSum):

http://www.topcoder.com/wiki/display/tc/SRM+467

Avoiding sum, multiplication overflow with modulo

To calculate (100003 - 200003*x + 300007*x*x*x) % 1000000 I would do:

1. y = x % 1000000
2. res = (100003 - (200003 * y) % 1000000 + (((300007 * y) % 1000000) * y) % 1000000) * y) % 1000000
3. if (x > 0 && res < 0) res += 1000000
4. if (x < 0 && res > 0) res -= 1000000

The way the formula is written we know that the result has to have the same sign as x. This is probably where you are going wrong, you omit step 3 and 4. Although as mentioned using 32 bit it can still overflow. You can omit step 4 if you need a positive modulo and just add 1000000 if res is negative.

---

Related Topics