-2

I want to calculate the modulus of a square: n^2 % m, where both n and m are large numbers (but less than the maximum of a 64 bits integer). The problem arrises when n^2 gets larger than the 64 bits maximum.

Is there an algorithm available to do this calculation? I know n^2 % m = (n % m)^2 %m, but this doesn't help me when m > n.

Spifff
  • 748
  • 6
  • 14
  • 1
    Yes, but depending on the language you're using to implement it there may be simpler tricks (for example in x64 assembly (and by extensions anything that has direct access to it) this can be implemented in 2 instructions) – harold Mar 19 '17 at 11:21
  • 1
    Does 2*m fit into a standard integer type? – kraskevich Mar 19 '17 at 11:48
  • As the previous comments state/imply, we need more information. Do you want a general algorithm that will work in (practically) any programming language on (practically) any hardware that have 64-bit integers, that will work for any positive values of `m` and `n` that fit in 64 bits? That is the most general case--more specific cases will allow faster or easier algorithms. – Rory Daulton Mar 19 '17 at 11:52
  • Why this limit of 64 ?, Do you want really an algorithm, or just a way to calculate. For the way, you can use Java BigInteger (and other libraries). For the algorithm, you can try to do everything in binary form. – guillaume girod-vitouchkina Mar 19 '17 at 11:57
  • One more question: are these signed or unsigned 64-bit integers? That matters. – Rory Daulton Mar 19 '17 at 12:06
  • I'm using Delphi, which only has signed 64 bits integers. I don't want to use biginteger libraries because speed is important. – Spifff Mar 19 '17 at 13:18

3 Answers3

2

Let n = k * 232 + j where j, k < 232. Then n ^ 2 % m = (264k2 + 2 * k * j * 232 + j2) % m

#include <iostream>

int main()
{
    uint64_t n = 17179874627;
    uint64_t m = 27778894627;

    uint64_t k = n >> 32;
    uint64_t j = n & 4294967295;

    uint64_t a = (k * k) % m;   // k^2
    a = (65536 * a) % m;        // 2^16 * k^2
    a = (65536 * a) % m;        // 2^32 * k^2
    a = (65536 * a) % m;        // 2^48 * k^2
    a = (65536 * a) % m;        // 2^64 * k^2

    uint64_t b = (j * 65536) % m;
    b = (b * 65536) % m;        // j * 2^32
    b = (b * k) % m;            // k * j * 2^32
    b = (2 * b) % m;            // 2 * k * j * 2^32

    uint64_t c = (j * j) % m;   // j^ 2

    std::cout << "Result " << (a + b + c) % m;
}
invisal
  • 11,075
  • 4
  • 33
  • 54
0

Based on invisal's answer, I expanded it into Delphi code for the situation (n1 * n2) % m:

k1 := n1 shr 30;
j1 := n1 and 1073741823; // = 2^30 - 1
k2 := n2 shr 30;
j2 := n2 and 1073741823;

a1 := (k1 * k2) mod m;
for i := 1 to 4 do
  a1 := (32768 * a1) mod m;

a2 := (k1 * j2) mod m;
for i := 1 to 2 do
  a2 := (32768 * a2) mod m;

a3 := (k2 * j1) mod m;
for i := 1 to 2 do
  a3 := (32768 * a3) mod m;

a4 := (j1 * j2) mod m;

Result := (a1 + a2 + a3 + a4) mod m;

Note that Delphi does not have an unsigned 64 bit integer, so I used 2^30 instead of 2^32.

Spifff
  • 748
  • 6
  • 14
0

Here is an algorithm that should work. I'll show it by example.

Let's say we need 29^2 mod 41, but cannot go beyond 256 (8 bits)

First, write 29 in binary: 11101

Now proceed like this: for each bit:

  1. double the result
  2. add (bitvalue *29) to result
  3. Take result mod 41

So starting with result = 0:

1 : 0 ; 29 ; 29
1 : 58 ; 87 ; 5
1 : 10 ; 39 ; 39
0 : 78 ; 78 ; 37
1 : 74 ; 103 ; 21

And of course 841 mod 41 does indeed equal 21, qed

tinnderbox
  • 1