```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 ``` ```# modular testing def gcd(m,n): while n <> 0: m, n = n, m % n return m def primes(n): ps, sieve = [], [True] * (n + 1) for p in range(2, n + 1): if sieve[p]: ps.append(p) for i in range(p * p, n + 1, p): sieve[i] = False return ps def isPrime(n): if n % 2 == 0: return n == 2 d = 3 while d * d <= n: if n % d == 0: return False d = d + 2 return True def inverse(x, m): a, b, u = 0, m, 1 while x > 0: q = b // x x, a, b, u = b % x, u, x, a - q * u if b == 1: return a % m raise ArithmeticError("must be coprime") def jacobi(a, p): a = a % p; t = 1 while a != 0: while a % 2 == 0: a = a / 2 if p % 8 in [3, 5]: t = -t a, p = p, a if a % 4 == 3 and p % 4 == 3: t = -t a = a % p return t if p == 1 else 0 def modSqrt(a, p): a = a % p if p % 4 == 3: x = pow(a, (p+1)/4, p) return x, p-x if p % 8 == 5: x = pow(a, (p+3)/8, p) c = pow(x, 2, p) if a != c: x = (x * pow(2, (p-1)/4, p)) % p return x, p-x d = 2 while jacobi(d, p) != -1: d += 1 t, s = p-1, 0 while t % 2 == 0: t, s = t / 2, s + 1 at, dt, m = pow(a, t, p), pow(d, t, p), 0 for i in range(0, s): if pow(at * pow(dt, m), pow(2, s-i-1), p) == p - 1: m = m + pow(2, i) x = (pow(dt, m) * pow(a, (t+1)/2, p)) % p return x, p-x def testPrime(n): ps = primes(n) for i in range(2,n): if i in ps: if not isPrime(i): print "prime not isPrime", i else: if isPrime(i): print "isPrime not prime", i testPrime(200) def testInverse(n): for m in range(2, n): for i in range(0, m): if gcd(i,m) == 1: j = inverse(i,m) if (i*j)%m <> 1: print "inverse fails", i, m testInverse(100) def testQuadRes(n): for p in primes(n)[1:]: for a in range(1,p): if jacobi(a,p) == 1: x, y = modSqrt(a,p) if (x*x)%p <> a: print "failure", a, p, x if (y*y)%p <> a: print "failure", a, p, y testQuadRes(100) def modSqrt(a, p): a = a % p if p % 4 == 3: x = pow(a, (p+1)/4, p) return x, p-x if p % 8 == 5: x = pow(a, (p+3)/8, p) c = pow(x, 2, p) if a != c: x = (x * pow(2, (p-1)/4, p)) % p return x, p-x d = 2 while jacobi(d, p) != -1: d += 1 t, s = p-1, 0 while t % 2 == 0: t, s = t / 2, s + 1 at, dt, m = pow(a, t, p), pow(d, t, p), 0 for i in range(0, s): if pow(at * pow(dt, m), pow(2, s-i-1), p) == p - 1: m = m + pow(2, i) x = (pow(dt, m/2) * pow(a, (t+1)/2, p)) % p return x, p-x testQuadRes(100) ```
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 ``` ```failure 2 17 1 failure 2 17 16 failure 4 17 1 failure 4 17 16 failure 8 17 1 failure 8 17 16 failure 9 17 1 failure 9 17 16 failure 13 17 1 failure 13 17 16 failure 15 17 1 failure 15 17 16 failure 16 17 1 failure 16 17 16 failure 2 41 31 failure 2 41 10 failure 4 41 18 failure 4 41 23 failure 5 41 23 failure 5 41 18 failure 8 41 25 failure 8 41 16 failure 9 41 40 failure 9 41 1 failure 20 41 4 failure 20 41 37 failure 21 41 4 failure 21 41 37 failure 23 41 10 failure 23 41 31 failure 25 41 37 failure 25 41 4 failure 31 41 16 failure 31 41 25 failure 32 41 40 failure 32 41 1 failure 33 41 25 failure 33 41 16 failure 36 41 23 failure 36 41 18 failure 39 41 31 failure 39 41 10 failure 40 41 1 failure 40 41 40 failure 3 73 64 failure 3 73 9 failure 6 73 4 failure 6 73 69 failure 9 73 8 failure 9 73 65 failure 12 73 55 failure 12 73 18 failure 18 73 37 failure 18 73 36 failure 19 73 32 failure 19 73 41 failure 23 73 16 failure 23 73 57 failure 24 73 8 failure 24 73 65 failure 25 73 37 failure 25 73 36 failure 27 73 1 failure 27 73 72 failure 35 73 2 failure 35 73 71 failure 36 73 16 failure 36 73 57 failure 38 73 2 failure 38 73 71 failure 41 73 55 failure 41 73 18 failure 46 73 1 failure 46 73 72 failure 48 73 37 failure 48 73 36 failure 49 73 8 failure 49 73 65 failure 50 73 16 failure 50 73 57 failure 54 73 32 failure 54 73 41 failure 57 73 4 failure 57 73 69 failure 61 73 55 failure 61 73 18 failure 65 73 64 failure 65 73 9 failure 67 73 4 failure 67 73 69 failure 69 73 2 failure 69 73 71 failure 70 73 64 failure 70 73 9 failure 71 73 32 failure 71 73 41 failure 72 73 1 failure 72 73 72 failure 5 89 9 failure 5 89 80 failure 9 89 53 failure 9 89 36 failure 10 89 42 failure 10 89 47 failure 11 89 73 failure 11 89 16 failure 17 89 69 failure 17 89 20 failure 18 89 10 failure 18 89 79 failure 20 89 18 failure 20 89 71 failure 21 89 40 failure 21 89 49 failure 22 89 44 failure 22 89 45 failure 25 89 81 failure 25 89 8 failure 34 89 55 failure 34 89 34 failure 36 89 17 failure 36 89 72 failure 40 89 84 failure 40 89 5 failure 42 89 68 failure 42 89 21 failure 44 89 57 failure 44 89 32 failure 47 89 21 failure 47 89 68 failure 49 89 5 failure 49 89 84 failure 50 89 22 failure 50 89 67 failure 53 89 72 failure 53 89 17 failure 55 89 34 failure 55 89 55 failure 57 89 11 failure 57 89 78 failure 68 89 49 failure 68 89 40 failure 69 89 71 failure 69 89 18 failure 71 89 79 failure 71 89 10 failure 72 89 20 failure 72 89 69 failure 73 89 85 failure 73 89 4 failure 79 89 47 failure 79 89 42 failure 80 89 36 failure 80 89 53 failure 81 89 50 failure 81 89 39 failure 84 89 80 failure 84 89 9 failure 85 89 87 failure 85 89 2 failure 87 89 25 failure 87 89 64 failure 88 89 88 failure 88 89 1 failure 2 97 49 failure 2 97 48 failure 3 97 65 failure 3 97 32 failure 4 97 73 failure 4 97 24 failure 6 97 81 failure 6 97 16 failure 8 97 85 failure 8 97 12 failure 9 97 54 failure 9 97 43 failure 11 97 53 failure 11 97 44 failure 12 97 89 failure 12 97 8 failure 16 97 91 failure 16 97 6 failure 18 97 27 failure 18 97 70 failure 22 97 75 failure 22 97 22 failure 24 97 93 failure 24 97 4 failure 25 97 66 failure 25 97 31 failure 27 97 18 failure 27 97 79 failure 31 97 72 failure 31 97 25 failure 32 97 94 failure 32 97 3 failure 33 97 50 failure 33 97 47 failure 36 97 62 failure 36 97 35 failure 43 97 88 failure 43 97 9 failure 44 97 86 failure 44 97 11 failure 47 97 64 failure 47 97 33 failure 48 97 95 failure 48 97 2 failure 49 97 2 failure 49 97 95 failure 50 97 33 failure 50 97 64 failure 53 97 11 failure 53 97 86 failure 54 97 9 failure 54 97 88 failure 62 97 36 failure 62 97 61 failure 64 97 47 failure 64 97 50 failure 65 97 3 failure 65 97 94 failure 66 97 25 failure 66 97 72 failure 70 97 79 failure 70 97 18 failure 72 97 31 failure 72 97 66 failure 73 97 4 failure 73 97 93 failure 75 97 22 failure 75 97 75 failure 79 97 70 failure 79 97 27 failure 81 97 6 failure 81 97 91 failure 85 97 8 failure 85 97 89 failure 86 97 44 failure 86 97 53 failure 88 97 43 failure 88 97 54 failure 89 97 12 failure 89 97 85 failure 91 97 16 failure 91 97 81 failure 93 97 24 failure 93 97 73 failure 94 97 32 failure 94 97 65 failure 95 97 48 failure 95 97 49 failure 96 97 96 failure 96 97 1 ```