lib/msqr.py

   1 # from http://eli.thegreenplace.net/2009/03/07/computing-modular-square-roots-in-python/
   2
   3 def modular_sqrt(a, p):
   4     """ Find a quadratic residue (mod p) of 'a'. p
   5     must be an odd prime.
   6
   7     Solve the congruence of the form:
   8     x^2 = a (mod p)
   9     And returns x. Note that p - x is also a root.
  10
  11     0 is returned is no square root exists for
  12     these a and p.
  13
  14     The Tonelli-Shanks algorithm is used (except
  15     for some simple cases in which the solution
  16     is known from an identity). This algorithm
  17     runs in polynomial time (unless the
  18     generalized Riemann hypothesis is false).
  19     """
  20     # Simple cases
  21     #
  22     if legendre_symbol(a, p) != 1:
  23         return 0
  24     elif a == 0:
  25         return 0
  26     elif p == 2:
  27         return p
  28     elif p % 4 == 3:
  29         return pow(a, (p + 1) / 4, p)
  30
  31     # Partition p-1 to s * 2^e for an odd s (i.e.
  32     # reduce all the powers of 2 from p-1)
  33     #
  34     s = p - 1
  35     e = 0
  36     while s % 2 == 0:
  37         s /= 2
  38         e += 1
  39
  40     # Find some 'n' with a legendre symbol n|p = -1.
  41     # Shouldn't take long.
  42     #
  43     n = 2
  44     while legendre_symbol(n, p) != -1:
  45         n += 1
  46
  47     # Here be dragons!
  48     # Read the paper "Square roots from 1; 24, 51,
  49     # 10 to Dan Shanks" by Ezra Brown for more
  50     # information
  51     #
  52
  53     # x is a guess of the square root that gets better
  54     # with each iteration.
  55     # b is the "fudge factor" - by how much we're off
  56     # with the guess. The invariant x^2 = ab (mod p)
  57     # is maintained throughout the loop.
  58     # g is used for successive powers of n to update
  59     # both a and b
  60     # r is the exponent - decreases with each update
  61     #
  62     x = pow(a, (s + 1) / 2, p)
  63     b = pow(a, s, p)
  64     g = pow(n, s, p)
  65     r = e
  66
  67     while True:
  68         t = b
  69         m = 0
  70         for m in xrange(r):
  71             if t == 1:
  72                 break
  73             t = pow(t, 2, p)
  74
  75         if m == 0:
  76             return x
  77
  78         gs = pow(g, 2 ** (r - m - 1), p)
  79         g = (gs * gs) % p
  80         x = (x * gs) % p
  81         b = (b * g) % p
  82         r = m
  83
  84 def legendre_symbol(a, p):
  85     """ Compute the Legendre symbol a|p using
  86     Euler's criterion. p is a prime, a is
  87     relatively prime to p (if p divides
  88     a, then a|p = 0)
  89
  90     Returns 1 if a has a square root modulo
  91     p, -1 otherwise.
  92     """
  93     ls = pow(a, (p - 1) / 2, p)
  94     return -1 if ls == p - 1 else ls