Restore a number from several its remainders (chinese remainder theorem)

Question

I have a long integer number, but it is stored not in decimal form, but as set of remainders.

So, I have not the N number, but set of such remainders:

r_1 = N % 2147483743
r_2 = N % 2147483713
r_3 = N % 2147483693
r_4 = N % 2147483659
r_5 = N % 2147483647
r_6 = N % 2147483629

I know, that N is less than multiplication of these primes, so chinese remainder theorem does work here ( http://en.wikipedia.org/wiki/Chinese_remainder_theorem ).

How can I restore N in decimal, if I have this 6 remainders? The wonderful will be any program to do this (C/C+GMP/C++/perl/java/bc).

For example, what minimal N can have this set of remainders:

r_1 = 1246736738 (% 2147483743)
r_2 = 748761 (% 2147483713)
r_3 = 1829651881 (% 2147483693)
r_4 = 2008266397 (% 2147483659)
r_5 = 748030137 (% 2147483647)
r_6 = 1460049539 (% 2147483629)

Waistcoat

What? no dc? Oh, well ... +1 for bc :) – Highsmith

Answer 1

Here the code (C+GMP), based on this LGPL code by Ben Lynn blynn@github; stanford of Garner algorithm (found with RIP Google Code Search by query garner mpz_t): https://github.com/blynn/pbc/blob/master/guru/indexcalculus.c (Part of his The PBC (Pairing-Based Crypto) library)

Compile with gcc -std=c99 -lgmp. Also change size for your case.

#include <gmp.h>
#include <stdlib.h>
#include <stdio.h>
#include <malloc.h>

// Garner's Algorithm.
// See Algorithm 14.71, Handbook of Cryptography.

//    x - result    v residuals    m - primes   t-size of vectors
static void CRT(mpz_t x, mpz_ptr *v, mpz_ptr *m, int t) {
  mpz_t u;
  mpz_t C[t];
  int i, j;

  mpz_init(u);
  for (i=1; i<t; i++) {
    mpz_init(C[i]);
    mpz_set_ui(C[i], 1);
    for (j=0; j<i; j++) {
      mpz_invert(u, m[j], m[i]);
      mpz_mul(C[i], C[i], u);
      mpz_mod(C[i], C[i], m[i]);
    }
  }
  mpz_set(u, v[0]);
  mpz_set(x, u);
  for (i=1; i<t; i++) {
    mpz_sub(u, v[i], x);
    mpz_mul(u, u, C[i]);
    mpz_mod(u, u, m[i]);
    for (j=0; j<i; j++) {
      mpz_mul(u, u, m[j]);
    }
    mpz_add(x, x, u);
  }

  for (i=1; i<t; i++) mpz_clear(C[i]);
  mpz_clear(u);
}

const int size=6; // Change this please

int main()
{
    mpz_t res;
    mpz_ptr t[size], p[size];
    for(int i=0;i<size;i++) { 
        t[i]=(mpz_ptr)malloc(sizeof(mpz_t));
        p[i]=(mpz_ptr)malloc(sizeof(mpz_t));
        mpz_init(p[i]);
        mpz_init(t[i]);
    }
    mpz_init(res);

    for(int i=0;i<size;i++){
        unsigned long rr,pp;
        scanf("%*c%*c%*c = %lu (%% %lu)\n",&rr,&pp);
        printf("Got %lu res on mod %% %lu \n",rr,pp);
        mpz_set_ui(p[i],pp);
        mpz_set_ui(t[i],rr);
    }

    CRT(res,t,p,size);

    gmp_printf("N = %Zd\n", res);
}

Example is solved:

$ ./a.out
r_1 = 1246736738 (% 2147483743)
r_2 = 748761 (% 2147483713)
r_3 = 1829651881 (% 2147483693)
r_4 = 2008266397 (% 2147483659)
r_5 = 748030137 (% 2147483647)
r_6 = 1460049539 (% 2147483629)

Got 1246736738 res on mod % 2147483743 
Got 748761 res on mod % 2147483713 
Got 1829651881 res on mod % 2147483693 
Got 2008266397 res on mod % 2147483659 
Got 748030137 res on mod % 2147483647 
Got 1460049539 res on mod % 2147483629 
N = 703066055325632897509116263399480311

N is 703066055325632897509116263399480311

Waistcoat

Answer 2

The article you link already provides a constructive algorithm to find the solution.

Basically, for each i you solve integer equation ri*ni + si*(N/ni) = 1 where N = n1*n2*n3*.... The ri and si are unknowns here. This can be solved by extended euclidean algorithm. It's very popular and you'll have no problem finding sample implementations in any language.

Then, assuming ei = si*(N/ni), the answer is sum(ei*ai) for every i.
All this is described in that article, with proof and example.

Berezniki

But I can't program this algorithm. Can you help? – Waistcoat

@Waistcoat What exactly is the difficulty? Most people here are busy enough and won't write the complete solution for you, but they might help with specific problem areas. – Berezniki

Answer 3

Here's a Python 3 implementation, based on this Rosetta Code task: https://rosettacode.org/wiki/Chinese_remainder_theorem

from functools import reduce
from operator import mul    

def chinese_remainder(n, a):
    """
    Chinese Remainder Theorem.

    :param n: list of pairwise relatively prime integers
    :param a: remainders when x is divided by n
    """
    s = 0
    prod = reduce(mul, n)
    for n_i, a_i in zip(n, a):
        p = prod // n_i
        s += a_i * inverse(p, n_i) * p
    return s % prod    

def inverse(a, b):
    """
    Modular multiplicative inverse.
    """
    b0 = b
    x0, x1 = 0, 1
    if b == 1:
        return 1
    while a > 1:
        q = a // b
        a, b = b, a % b
        x0, x1 = x1 - q * x0, x0
    if x1 < 0:
        x1 += b0
    return x1    

n = [2147483743, 2147483713, 2147483693, 2147483659, 2147483647, 2147483629]
a = [1246736738, 748761, 1829651881, 2008266397, 748030137, 1460049539]

print(chinese_remainder(n, a))  # 703066055325632897509116263399480311

A nice feature of Python is that it supports arbitrarily large integers naturally.

Ira

Answer 4

Just adding a little bit of theory from here and its verbatim implementation with python:

from math import prod
    
def CRT_residue(a, m):
    # assert len(a) == len(m)
    n, P = len(a), prod(m)
    M = [P // m[i] for i in range(n)]  # M_i s
    # modular inverse of M_i to be computed with Extended Euclid algorithm
    # extract the Bézout coefficient corresponding to M_i
    N = [gcd_ext(M[i], m[i])[1] for i in range(n)] # N_i s
    return sum([a[i]*M[i]*N[i] for i in range(n)]) % P

def gcd_ext(a, b): # Extended Euclid algorithm
    x0, y0 = 1, 0
    x1, y1 = 0, 1
    while b:
        q = a // b
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
        a, b = b, a % b
    return a, x0, y0 # return GCD and the Bézout coefficients

# given example
m = [2147483743, 2147483713, 2147483693, 2147483659, 2147483647, 2147483629]
a = [1246736738, 748761, 1829651881, 2008266397, 748030137, 1460049539]
print(CRT_residue(a, m))
# 703066055325632897509116263399480311

Another possibility is to use the Garner's algorithm.

Ellette