9

I need to calculate a*a mod n but a is fairly large, resulting in overflow when I square it. Doing ((a % n)*(a % n)) % n doesn't work because (n-1)2 can overflow. This is in C++ and I'm using int64_t

Edit:

Example value: a = 821037907258 and n = 800000000000, which overflows if you square it.

I am using DevCPP and I've already tried getting big-integer libraries working to no avail.

Edit 2:

No, there's no pattern to these numbers.

Community
  • 1
  • 1
WhatsInAName
  • 969
  • 2
  • 9
  • 12
  • Are you able to use a big-integer library? – Oliver Charlesworth Apr 09 '12 at 16:06
  • I tried installing BigInteger, GMP, etc, but I use Windows and those libraries seem to require a billion steps that all fail for me at various points as I crawl back the list of needed dependencies. – WhatsInAName Apr 09 '12 at 16:07
  • Any number over 32bits (~4billion unsigned, ~2billion signed) squared will overflow. You can detect this condition beforehand and handle it yourself using some relatively simple 'big integer' routines. – std''OrgnlDave Apr 09 '12 at 16:08
  • Yes, but I think the question is whether there is a mathematical simplification that can be applied. – Matt Apr 09 '12 at 16:08
  • 1
    @Mat Simplification? Hmph, no. But the multiplication can be handled with shifts and adds, which isn't too hard to implement. – std''OrgnlDave Apr 09 '12 at 16:09
  • Can you use `double` and cast that after calculation to `int64`? – juergen d Apr 09 '12 at 16:10
  • You'll probably need to look at an arbitrary precision math library as found in some of the various crypto packages (e.g. OpenSSL, libgcrypt). Big integers are used a lot in crypto. – jbruni Apr 09 '12 at 16:11
  • 3
    If you are use `g++`, you could use `__uint128_t` and `__int128_t` extensions. – Sergey Kalinichenko Apr 09 '12 at 16:13
  • 2
    @juergend Double is 64-bits, which is what, 13 bits for the mantissa? That leaves 2 bits for signs (one for integer part, one for mantissa) and you get like, 51 bits of integer precision. That won't help the overflow much, in fact, it will overflow sooner - assuming the OP wants accuracy – std''OrgnlDave Apr 09 '12 at 16:13
  • @OrgnlDave: Actually the integer precision of `double` is 53 bit. – kennytm Apr 09 '12 at 16:14
  • @KennyTM N is fixed, and it fits within an int64 fine. No problems there. – WhatsInAName Apr 09 '12 at 16:20
  • @Mat not sure how you're getting that the remainder is a^2-n. In a smaller case, 4^2 mod 3 = 1. But 4^2-3 = 13. – WhatsInAName Apr 09 '12 at 16:22
  • Any bigint library should work fine with this, or you had used it the wrong way. But what's the reason for avoiding `__int128_t` when it's available? – phuclv Mar 31 '15 at 03:39

5 Answers5

15

If you can't use a big-integer library, and you don't have a native uint128_t (or similar), you'll need to do this manually.

One option is to express a as the sum of two 32-bit quantities, i.e. a = 232b + c, where b contains the 32 msbs, and c contains the 32 lsbs. Squaring is then a set of four cross-multiplications; each result is guaranteed to fit into a 64-bit type. You then do the modulo operation as you recombine the individual terms (carefully taking into account the shifts needed to realign everything).

Oliver Charlesworth
  • 267,707
  • 33
  • 569
  • 680
  • @OrgnlDave: Nope, each term of the cross-multiplication fits into 64 bits. – Oliver Charlesworth Apr 09 '12 at 16:11
  • Sorry I am a beginner and have no idea what this means – WhatsInAName Apr 09 '12 at 16:15
  • @WhatsInAName: No idea about which part? Are you comfortable with the fact that a 64-bit integer can be expressed as the combination of two 32-bit integers? And that (a+b) squared can be expressed as the sum of four multiplications? – Oliver Charlesworth Apr 09 '12 at 16:16
  • @OliCharlesworth Maybe show a sample implementation with 4- or 8-bit integers so that it's short and concise but clear? – std''OrgnlDave Apr 09 '12 at 16:17
  • @OrgnlDave: This is the same thing, really. Producing that 128-bit integer requires 4 individual multiplications. – Oliver Charlesworth Apr 09 '12 at 16:17
  • @OilCharlesworth Splitting some value into a_hi and a_lo (how do I split?) and how to do the necessary multiplications and recombinations/shifts... all that stuff is beyond my scope. – WhatsInAName Apr 09 '12 at 16:19
  • 3
    @WhatsInAName Basically, Oli is telling you to implement long multiplication and long division the same way you would do it grade-school. – Mysticial Apr 09 '12 at 16:21
  • 2
    @Mysticial I'm pretty sure they don't teach bitshifting and modular arithmetic in gradeschool. Please keep the snarkiness to a minimum. That is clearly not what I am asking about. – WhatsInAName Apr 09 '12 at 16:24
  • @WhatsInAName: Oli means rewriting *a* into the form *a* = 2³² *b* + *c* where 0 ≤ *b*, *c* < 2³². *b* is `a_hi` and *c* is `a_lo`. – kennytm Apr 09 '12 at 16:26
  • @KennyTM: Yep. I've modified my answer, because that's a lot clearer. – Oliver Charlesworth Apr 09 '12 at 16:27
  • 1
    @WhatsInAName I'm not trying to be snarky - I apologize if you took it that way. In grade-school arithmetic, you "shift" the digits. In computers, you "shift" the bits. It's all the same except in a different base. – Mysticial Apr 09 '12 at 16:27
  • I think I'm going to have to just figure out a different way to handle this entire issue in the first place. Sorry for wasting everyone's time. I'll mark this as the answer anyway because I'm sure it's right, I just have no clue how to do this. – WhatsInAName Apr 09 '12 at 16:28
  • 1
    Is there any guarantee that `mod(2^64, n)*b*b` (or each subterm as we calculate the end result) won't overflow? Or should we recursively apply the multiplication function described above? – vhallac Apr 09 '12 at 16:37
  • @vhallac: Yes, it's possible that `mod(term,n)*mod(1< – Oliver Charlesworth Apr 09 '12 at 16:42
  • What is msbs and lsbs? Can you use his number in an example? Your answer doesn't seem to apply to the question... Please expand on it or provide a sample... – Jason Goemaat Feb 08 '15 at 08:04
  • @JasonGoemaat: https://www.google.co.uk/search?q=msb+lsb. It does apply to the question; the crux is "decompose your input into chunks that fit within the types you have available". – Oliver Charlesworth Feb 08 '15 at 10:51
  • Ah, never seen that before: https://www.google.com/search?q=msbs - Can you expand on `You then do the modulo operation as you recombine the individual terms (carefully taking into account the shifts needed to realign everything).` and give a sample? – Jason Goemaat Feb 08 '15 at 13:41
  • @JasonGoemaat: I haven't got the patience to work a full example ;) However, the result of the cross-multiplication gives us something of the form `2^64 * a + 2^32 * (b + c) + d`. Each of these remaining operations can be done in modulo arithmetic without overflowing. – Oliver Charlesworth Feb 08 '15 at 17:06
4

I know you no longer need this, and there is an alternative solution, but I want to add an alternative method to implement it. It provides two different techniques: the double and add algorithm, and the method to handle mod(a + b, n) with overflow detection.

Double and add algorithm is usually used in fields where multiplication is not possible or too costly to calculate directly (such as elliptic curves), but we could adopt it to handle it in our situation to handle overflows instead.

The following code is probably slower than the accepted solution (even when you optimize it), but if speed is not critical, you may prefer it for clarity.

unsigned addmod(unsigned x, unsigned y, unsigned n)
{
    // Precondition: x<n, y<n
    // If it will overflow, use alternative calculation
    if (x + y <= x) x = x - (n - y);
    else x = (x + y) % n;
    return x;
}

unsigned sqrmod(unsigned a, unsigned n)
{
    unsigned b;
    unsigned sum = 0;

    // Make sure original number is less than n
    a = a % n;

    // Use double and add algorithm to calculate a*a mod n
    for (b = a; b != 0; b >>= 1) {
        if (b & 1) {
            sum = addmod(sum, a, n);
        }
        a = addmod(a, a, n);
    }
    return sum;
}
vhallac
  • 13,301
  • 3
  • 25
  • 36
  • `x = x - (n - y)` - are you sure that does anything useful? – Oliver Charlesworth Apr 09 '12 at 19:05
  • 1
    @OliCharlesworth It calculates `(x+y)%n` without overflowing. The non-overflow is simple to prove. both x and y are less than `n` (and hence maximum word size), and we are only subtracting numbers, so they only go downwards, which means they can never overflow. To prove `x - (n - y)` is correct is a little trickier. First, since it overflows, `x+y` is clearly greater than `n`. Also, since both `x` and `y` are less than `n`, `x+y` is less than `2*n`, so the modulus is equal to `(x+y)-n`. If you write `x+y` as `x-(-y)` and rearrange terms, you get the formula I used. – vhallac Apr 09 '12 at 20:17
  • @OliCharlesworth I think you can use the alternative formula `(x+y)+2^W-n` instead of `x-(n-y)`, but I didn't try it. This would have the advantage of letting us calculate `x+y` only once in `addmul`. – vhallac Apr 09 '12 at 20:23
  • `addmod()` needs a comment saying a precondition is that `x < n` and `y < n`. – Craig McQueen Oct 29 '13 at 05:52
  • @CraigMcQueen Yes. It is not obvious. I am editing the answer to add it. – vhallac Oct 30 '13 at 17:51
1

Here is a double-and-add implementation of a multiplication a * b % m, without overflows, whatever the size of a, b and m.

(Note that the res -= m and temp_b -= m lines rely on 64-bit unsigned integer overflow in order to give the expected results. This should be fine since unsigned integer overflow is well-defined in C and C++. For this reason it's important to use unsigned integer types.)

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
    uint64_t res = 0;
    uint64_t temp_b;

    /* Only needed if b may be >= m */
    if (b >= m) {
        if (m > UINT64_MAX / 2u)
            b -= m;
        else
            b %= m;
    }

    while (a != 0) {
        if (a & 1) {
            /* Add b to res, modulo m, without overflow */
            if (b >= m - res) /* Equiv to if (res + b >= m), without overflow */
                res -= m;
            res += b;
        }
        a >>= 1;

        /* Double b, modulo m */
        temp_b = b;
        if (b >= m - b)       /* Equiv to if (2 * b >= m), without overflow */
            temp_b -= m;
        b += temp_b;
    }
    return res;
}

This is my modification of another answer to another similar question.

Community
  • 1
  • 1
Craig McQueen
  • 41,871
  • 30
  • 130
  • 181
-3

You can implement the multiplication yourself, adding n each run, and mod'ing the result right away.

Alberto
  • 499
  • 4
  • 23
  • 3
    You realize he's talking about numbers over 4 billion, right? That's what it takes to overflow a 64-bit int... – std''OrgnlDave Apr 09 '12 at 16:07
  • It would take way too long for me to iterate through the a-value to do this – WhatsInAName Apr 09 '12 at 16:08
  • 1
    @OrgnlDave: that would be enough to overflow a 32-bit int. For a 64-bit int, you need something about 4 billion times that size. (Or do you mean the two inputs are over 4 billion, so multiplying them overflows a 64-bit int?) – Jerry Coffin Apr 09 '12 at 16:09
-5

I really think that ((a % n)*(a % n)) % n should work for positive integers. Why do you think it doesn't work? Do you have a counter-example? If n could be negative, then the % operator is undefined.

phuclv
  • 37,963
  • 15
  • 156
  • 475
danimator
  • 16
  • 5
  • 2
    a is smaller than n. a^2 is larger than n. However, a^2 overflows. – WhatsInAName Apr 09 '12 at 16:27
  • @WhatsInAName: If `(n-1)` squared can overflow, you should have said that in the question, because the question doesn't say that. Your question merely says that `a` squared overflows. This answer depends on `n` squared not overflowing. – Mooing Duck Apr 09 '12 at 16:37
  • I didn't know it mattered because I'm not squaring n. But yes, n-1 squared would overflow, too. – WhatsInAName Apr 09 '12 at 16:40
  • 1
    If `n` is very large, then `a%n` can be very large. – Oliver Charlesworth Apr 09 '12 at 16:43
  • @WhatsInAName: if `(n-1)` squared didn't overflow, then this answer would work just fine. Please put that information in the question, so that future answerers know. – Mooing Duck Apr 09 '12 at 16:46
  • 2
    @MooingDuck: I don't think that's necessary. The fact that the OP is asking the question at all implies that `n-1` squared can overflow. This answer is based on an unwarranted assumption. – Oliver Charlesworth Apr 09 '12 at 16:49
  • Well I mean of course the answer would work if that didn't overflow (because a^2 would be smaller than that) and I've already stated that a%n does nothing. – WhatsInAName Apr 09 '12 at 16:50
  • @OliCharlesworth: He said that the process "didn't work", (which I admit, is a good sign that n-1 squared can overflow), but I think the question would be improved if he used better wording that "didn't work". There's many questions here that say "didn't work" when they mean "I copy-pasted this and got a compiler error and gave up". – Mooing Duck Apr 09 '12 at 17:10