writeups/icecrax/asis2019/tonelli-shanks.py

   1 def legendre(a, p):
   2     return pow(a, (p - 1) // 2, p)
   3
   4 def tonelli(n, p):
   5     assert legendre(n, p) == 1, "not a square (mod p)"
   6     q = p - 1
   7     s = 0
   8     while q % 2 == 0:
   9         q //= 2
  10         s += 1
  11     if s == 1:
  12         return pow(n, (p + 1) // 4, p)
  13     for z in range(2, p):
  14         if p - 1 == legendre(z, p):
  15             break
  16     c = pow(z, q, p)
  17     r = pow(n, (q + 1) // 2, p)
  18     t = pow(n, q, p)
  19     m = s
  20     t2 = 0
  21     while (t - 1) % p != 0:
  22         t2 = (t * t) % p
  23         for i in range(1, m):
  24             if (t2 - 1) % p == 0:
  25                 break
  26             t2 = (t2 * t2) % p
  27         b = pow(c, 1 << (m - i - 1), p)
  28         r = (r * b) % p
  29         c = (b * b) % p
  30         t = (t * c) % p
  31         m = i
  32     return r
  33
  34 H = 381704527450191606347421195235742637659723827441243208291869156144963
  35 factors = [1504073, 20492753, 59833457464970183, 467795120187583723534280000348743236593] # all of them are coprime
  36
  37 # The Tonelli–Shanks algorithm is used in modular arithmetic to solve for
  38 # r in a congruence of the form r^2 ≡ n (mod p), where p is a prime number.
  39 # Note: it can be applied only if " Euler's criterion" is fullfiled or "Legendre symbol" = 1
  40 if __name__ == '__main__':
  41     ttest = [(H,factors[0]),(H,factors[1]), (H,factors[2]), (H,factors[3])]
  42     for n, p in ttest:
  43         print("=====================================================")
  44         print("Legendre symbol: ",  legendre(n, p))
  45         r = tonelli(n, p)
  46         assert (r * r - n) % p == 0
  47         print("n = %d p = %d" % (n, p))
  48         print("\t  roots : %d %d" % (r, p - r))
  49
  50
  51 # Results:
  52 # p0 = 1504073:
  53 # 399666
  54 # 1104407
  55
  56 # p1 = 20492753:
  57 # 7111848
  58 # 13380905
  59
  60 # p2 = 59833457464970183:
  61 # 34240854883018057
  62 # 25592602581952126
  63
  64 # p3 = 467795120187583723534280000348743236593:
  65 # 308269479959806774875048102517512730884
  66 # 159525640227776948659231897831230505709