sys/src/libauthsrv/msqrt.mp

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