Twierdzenie o resztach wydaje się dość proste, ale to wrażenie szybko mija, kiedy człowiek zacznie się wczytywać w zawiłe matematyczne symbole na Wikipedii. Muszę sobie tu kilka ważnych rzeczy zanotować na przyszłość.

Definicja

Układ kongruencji:

x = a1 (mod n1) 
x = a2 (mod n2) 
x = a3 (mod n3) 
... 
x = an (mod nk) (1)

przy założeniu, że n1, n2, ..., nn są względnie pierwsze, spełnia dokładnie jedna liczba x taka, że:

1 <= x <= n1*n2*n3*...*nk

Rozwiązanie

Założenia

Jeśli przyjmiemy:

M = n1*n2*n3*...*nk

oraz zdefiniujemy:

Mi = M / ni

To wtedy na podstawie faktu, że Mi oraz ni są względnie pierwsze, można zastosować rozszerzony algortym Euklidesa do znalezienia współczynników fi, gi (i <= k), że:

fi*ni + gi*Mi = 1 (2)

Rozwiązanie teretyczne

Wtedy x (rozwiązanie układu kongruencji) wynosi:

x = a1*g1*M1 + a2*g2*M2 + ... + ak*gk*Mk

Jak znaleźć współczynniki g1,g2, ..., gn? Z równania (2):

gi*Mi = (-fi)*ni + 1

a więc prawa strona, tak samo jak lewa (gi*Mi) przy dzieleniu przez ni da resztę 1:

gi*Mi = 1 (mod ni)

Aby więc znaleźć gi, należy zastosować algorytm poszukiwania odwrotności modulo dla wartości Mi w arytmetyce modulo ni.

Czego szukam?

Szukam więc takiej liczby g, że iloczyn Mi*g da resztę 1 w dzieleniu przez ni. Zbadam więc jaka wartość g (z przedziału od 1 do ni-1) da resztę 1 (mod ni). Oto algortym w modularInverse(g, ni) w Pythonie:

1
2
3
4
5
6

def modularInverse(g, mod): 
  gmod = g % mod 
  for x in range(1, mod): 
    if ((gmod * x) % mod == 1): 
      return x 
  return 1

Rozwiązanie praktyczne w Pythonie

Rozwiązaniem układu kongruencji (1) będzie zatem następujący kod (rezultat sumowania jest obliczony modulo M, ponieważ szukamy jego najmniejszej wartości; spełnia je także każda wielokrotność M, która daje taką samą wartość modulo M):

1
2
3
4
5
6
7
8

def crt(congruences):
    M = reduce(lambda acc, val: acc * val[0], congruences, 1)

    ssum = 0
    for (ni, ai) in congruences:
        Mi = floor(M / ni)
        ssum += ai * modularInverse(Mi, ni) * Mi
    return ssum % M