| 118 | | |
| 119 | | def randint(minvalue, maxvalue): |
| 120 | | """Returns a random integer x with minvalue <= x <= maxvalue""" |
| 121 | | # Safety - get a lot of random data even if the range is fairly |
| 122 | | # small |
| 123 | | min_nbits = 32 |
| 124 | | |
| 125 | | # The range of the random numbers we need to generate |
| 126 | | range = maxvalue - minvalue |
| 127 | | |
| 128 | | # Which is this number of bytes |
| 129 | | rangebytes = ceil(log(range, 2) / 8) |
| 130 | | |
| 131 | | # Convert to bits, but make sure it's always at least min_nbits*2 |
| 132 | | rangebits = max(rangebytes * 8, min_nbits * 2) |
| 133 | | |
| 134 | | # Take a random number of bits between min_nbits and rangebits |
| 135 | | nbits = random.randint(min_nbits, rangebits) |
| 136 | | |
| 137 | | return (read_random_int(nbits) % range) + minvalue |
| 138 | | |
| 139 | | def fermat_little_theorem(p): |
| 140 | | """Returns 1 if p may be prime, and something else if p definitely |
| 141 | | is not prime""" |
| 142 | | |
| 143 | | a = randint(1, p-1) |
| 144 | | return fast_exponentiation(a, p-1, p) |
| 145 | | |
| 146 | | def jacobi(a, b): |
| 147 | | """Calculates the value of the Jacobi symbol (a/b) |
| 148 | | """ |
| 149 | | |
| 150 | | if a == 0: return 0 |
| 151 | | if a == 1: return 1 |
| 152 | | |
| 153 | | if a % 2 == 0: |
| 154 | | if (b**2-1)/8 % 2 == 0: |
| 155 | | return jacobi(a/2, b) |
| 156 | | return -jacobi(a/2, b) |
| 157 | | |
| 158 | | if (a-1) * (b-1) / 4 % 2 == 0: |
| 159 | | return jacobi(b % a, a) |
| 160 | | |
| 161 | | return -jacobi(b % a, a) |
| 162 | | |
| 163 | | def jacobi_witness(x, n): |
| 164 | | """Returns False if n is an Euler pseudo-prime with base x, and |
| 165 | | True otherwise. |
| 166 | | """ |
| 167 | | |
| 168 | | j = jacobi(x, n) % n |
| 169 | | f = fast_exponentiation(x, (n-1)/2, n) |
| 170 | | |
| 171 | | if j == f: return False |
| 172 | | return True |
| 173 | | |
| 174 | | def randomized_primality_testing(n, k): |
| 175 | | """Calculates whether n is composite (which is always correct) or |
| 176 | | prime (which is incorrect with error probability 2**-k) |
| 177 | | |
| 178 | | Returns False if the number if composite, and True if it's |
| 179 | | probably prime. |
| 180 | | """ |
| 181 | | |
| 182 | | q = 0.5 # Property of the jacobi_witness function |
| 183 | | |
| 184 | | # t = int(math.ceil(k / log(1/q, 2))) |
| 185 | | t = ceil(k / log(1/q, 2)) |
| 186 | | for i in range(t+1): |
| 187 | | x = randint(1, n-1) |
| 188 | | if jacobi_witness(x, n): return False |
| 189 | | |
| 190 | | return True |
| 191 | | |
| 192 | | def is_prime(number): |
| 193 | | """Returns True if the number is prime, and False otherwise. |
| 194 | | |
| 195 | | >>> is_prime(42) |
| 196 | | 0 |
| 197 | | >>> is_prime(41) |
| 198 | | 1 |
| 199 | | """ |
| 200 | | |
| 201 | | """ |
| 202 | | if not fermat_little_theorem(number) == 1: |
| 203 | | # Not prime, according to Fermat's little theorem |
| 204 | | return False |
| 205 | | """ |
| 206 | | |
| 207 | | if randomized_primality_testing(number, 5): |
| 208 | | # Prime, according to Jacobi |
| 209 | | return True |
| 210 | | |
| 211 | | # Not prime |
| 212 | | return False |
| 213 | | |
| 214 | | |
| 215 | | def getprime(nbits): |
| 216 | | """Returns a prime number of max. 'math.ceil(nbits/8)*8' bits. In |
| 217 | | other words: nbits is rounded up to whole bytes. |
| 218 | | |
| 219 | | >>> p = getprime(8) |
| 220 | | >>> is_prime(p-1) |
| 221 | | 0 |
| 222 | | >>> is_prime(p) |
| 223 | | 1 |
| 224 | | >>> is_prime(p+1) |
| 225 | | 0 |
| 226 | | """ |
| 227 | | |
| 228 | | nbytes = int(math.ceil(nbits/8)) |
| 229 | | |
| 230 | | while True: |
| 231 | | integer = read_random_int(nbits) |
| 232 | | |
| 233 | | # Make sure it's odd |
| 234 | | integer |= 1 |
| 235 | | |
| 236 | | # Test for primeness |
| 237 | | if is_prime(integer): break |
| 238 | | |
| 239 | | # Retry if not prime |
| 240 | | |
| 241 | | return integer |
| 242 | | |
| 243 | | def are_relatively_prime(a, b): |
| 244 | | """Returns True if a and b are relatively prime, and False if they |
| 245 | | are not. |
| 246 | | |
| 247 | | >>> are_relatively_prime(2, 3) |
| 248 | | 1 |
| 249 | | >>> are_relatively_prime(2, 4) |
| 250 | | 0 |
| 251 | | """ |
| 252 | | |
| 253 | | d = gcd(a, b) |
| 254 | | return (d == 1) |
| 255 | | |
| 256 | | def find_p_q(nbits): |
| 257 | | """Returns a tuple of two different primes of nbits bits""" |
| 258 | | |
| 259 | | print 'finding p' |
| 260 | | p = getprime(nbits) |
| 261 | | while True: |
| 262 | | print 'finding q' |
| 263 | | q = getprime(nbits) |
| 264 | | if not q == p: break |
| 265 | | |
| 266 | | return (p, q) |
| 267 | | |
| 268 | | def extended_euclid_gcd(a, b): |
| 269 | | """Returns a tuple (d, i, j) such that d = gcd(a, b) = ia + jb |
| 270 | | """ |
| 271 | | |
| 272 | | if b == 0: |
| 273 | | return (a, 1, 0) |
| 274 | | |
| 275 | | q = abs(a % b) |
| 276 | | r = long(a / b) |
| 277 | | (d, k, l) = extended_euclid_gcd(b, q) |
| 278 | | |
| 279 | | return (d, l, k - l*r) |
| 280 | | |
| 281 | | # Main function: calculate encryption and decryption keys |
| 282 | | def calculate_keys(p, q, nbits): |
| 283 | | """Calculates an encryption and a decryption key for p and q, and |
| 284 | | returns them as a tuple (e, d)""" |
| 285 | | |
| 286 | | n = p * q |
| 287 | | phi_n = (p-1) * (q-1) |
| 288 | | |
| 289 | | while True: |
| 290 | | print 'c' |
| 291 | | # Make sure e has enough bits so we ensure "wrapping" through |
| 292 | | # modulo n |
| 293 | | e = getprime(max(8, nbits/2)) |
| 294 | | if are_relatively_prime(e, n) and are_relatively_prime(e, phi_n): break |
| 295 | | |
| 296 | | (d, i, j) = extended_euclid_gcd(e, phi_n) |
| 297 | | |
| 298 | | if not d == 1: |
| 299 | | raise Exception("e (%d) and phi_n (%d) are not relatively prime" % (e, phi_n)) |
| 300 | | |
| 301 | | if not (e * i) % phi_n == 1: |
| 302 | | raise Exception("e (%d) and i (%d) are not mult. inv. modulo phi_n (%d)" % (e, i, phi_n)) |
| 303 | | |
| 304 | | return (e, i) |
| 305 | | |
| 306 | | |
| 307 | | def gen_keys(nbits): |
| 308 | | """Generate RSA keys of nbits bits. Returns (p, q, e, d). |
| 309 | | """ |
| 310 | | |
| 311 | | while True: |
| 312 | | print 'find' |
| 313 | | (p, q) = find_p_q(nbits) |
| 314 | | print 'calculate' |
| 315 | | (e, d) = calculate_keys(p, q, nbits) |
| 316 | | |
| 317 | | # For some reason, d is sometimes negative. We don't know how |
| 318 | | # to fix it (yet), so we keep trying until everything is shiny |
| 319 | | if d > 0: break |
| 320 | | |
| 321 | | return (p, q, e, d) |
| | 95 | |
| 603 | | |
| | 366 | def getRandomNumber(low, high): |
| | 367 | if low >= high: |
| | 368 | raise AssertionError() |
| | 369 | howManyBits = numBits(high) |
| | 370 | howManyBytes = numBytes(high) |
| | 371 | lastBits = howManyBits % 8 |
| | 372 | while 1: |
| | 373 | bytes = getRandomBytes(howManyBytes) |
| | 374 | if lastBits: |
| | 375 | bytes[0] = bytes[0] % (1 << lastBits) |
| | 376 | n = bytesToNumber(bytes) |
| | 377 | if n >= low and n < high: |
| | 378 | return n |
| | 379 | |
| | 380 | |
| | 381 | |
| | 382 | def numberToBytes(n): |
| | 383 | howManyBytes = numBytes(n) |
| | 384 | bytes = createByteArrayZeros(howManyBytes) |
| | 385 | for count in range(howManyBytes-1, -1, -1): |
| | 386 | bytes[count] = int(n % 256) |
| | 387 | n >>= 8 |
| | 388 | return bytes |
| | 389 | |
| | 390 | def stringToBytes(s): |
| | 391 | bytes = createByteArrayZeros(0) |
| | 392 | bytes.fromstring(s) |
| | 393 | return bytes |
| | 394 | |
| | 395 | |
| | 396 | import array |
| | 397 | def createByteArraySequence(seq): |
| | 398 | return array.array('B', seq) |
| | 399 | def createByteArrayZeros(howMany): |
| | 400 | return array.array('B', [0] * howMany) |
| | 401 | |
| | 402 | def bytesToNumber(bytes): |
| | 403 | total = 0L |
| | 404 | multiplier = 1L |
| | 405 | for count in range(len(bytes)-1, -1, -1): |
| | 406 | byte = bytes[count] |
| | 407 | total += multiplier * byte |
| | 408 | multiplier *= 256 |
| | 409 | return total |
| | 410 | |
| | 411 | import math |
| | 412 | def numBits(n): |
| | 413 | if n==0: |
| | 414 | return 0 |
| | 415 | s = "%x" % n |
| | 416 | return ((len(s)-1)*4) + \ |
| | 417 | {'0':0, '1':1, '2':2, '3':2, |
| | 418 | '4':3, '5':3, '6':3, '7':3, |
| | 419 | '8':4, '9':4, 'a':4, 'b':4, |
| | 420 | 'c':4, 'd':4, 'e':4, 'f':4, |
| | 421 | }[s[0]] |
| | 422 | return int(math.floor(math.log(n, 2))+1) |
| | 423 | |
| | 424 | def numBytes(n): |
| | 425 | if n==0: |
| | 426 | return 0 |
| | 427 | bits = numBits(n) |
| | 428 | return int(math.ceil(bits / 8.0)) |
| 624 | | if 1: |
| 625 | | import time |
| 626 | | (pub,priv) = gen_pubpriv_keys(1024) |
| 627 | | if 1: |
| 628 | | pass |
| 629 | | if 0: |
| 630 | | (pub,priv) = ({'e': 3197601411731116245564574341873883965700807663084732800379376818822337768214254275271880556652543127861042149268538435129289119875433847550131501030387489L, 'n': 10010369815382334240798530070718437854966721919010883362239947723044651469296597049454697875962523892167961291050690952736009970933207150773523631818070845272169487515533531012017986671583117870935704791444203091193928433685036216475710517809854134259468721246954372989146266743601614235731605829521314456547427727170490340220674730134652649494795255568521600458207445549064461533102227417295842548674337715761362209701734789853469920681184060067484771558465649908138716260602630770402224385695023184536234858539880642604697788322269392260439339765412460944589951571879796455505538134420317166479947358156062601693631L}, {'q': 71965600310892660429596898829670710653460148180418095123381449844052583795523012040229068564004895983101524115981347946055591335247903874583493068435665706247560257067715751690858970054084093796406191757154228889454862489531588573733962647071944624542539072044760521321849853532919544810806179365743581833777L, 'p': 139099372090795607726534110588106724714748422622835743921578772078079056922772905042826722454890797681407412421239584468016215098928247799729322355221381502409668675456356288320338371140589835466925402699800592201081119770108520423450701555386146286042071274331899104700802281153510654832285005494641183409903L, 'd': 2303141903135405741196770316647349478451046851656460345238126953851335354457956582261071364112927329063172609648667594153989314522741915259780604774082900187421379799944807777478001765690017634728386784574675440128718738783220812396390380020488828055736972630556491383856716073773079586029351343732817861978611326631064317705218919924700761528865528832464733884127079042865050020476303293874373222849912706661252250304706695160160822827622086731394008802874568358702895903468129589361385276702313798198182238681153517515581052056557552995405472004415453043300567819425403784410663519000867938462388941313183838315329L}) |
| 631 | | if 0: |
| 632 | | #full |
| 633 | | t = time.time() |
| 634 | | message = 'serial '*5 |
| 635 | | #print 'cleartext ', message |
| 636 | | cypher = encrypt(message,pub) |
| 637 | | #print 'cyphertext: ',cypher |
| 638 | | #print 'decrypted', decrypt(cypher,priv) |
| 639 | | decrypt(cypher,priv) |
| 640 | | signed = sign(message,priv) |
| 641 | | #print 'signed', signed |
| 642 | | #print 'verified', message == verify(signed,pub) |
| 643 | | unblinder = getUnblinder(pub['n']) |
| 644 | | blinder = pow(invMod(unblinder, pub['n']), pub['e'],pub['n']) |
| 645 | | blinded = blind(message,blinder,pub) |
| 646 | | #print 'blinded', blinded |
| 647 | | #signedblind = sign(blinded,priv) |
| 648 | | signedblind = encrypt_int(blinded, priv['d'], priv['p']*priv['q']) |
| 649 | | #print 'signedblind', signedblind |
| 650 | | #unblinded = unblind(signedblind,unblinder,pub) |
| 651 | | unblinded = (signedblind * unblinder) % pub['n'] |
| 652 | | #print 'unblinded', unblinded |
| 653 | | #print 'verified', message == verify(unblinded,pub) |
| 654 | | print time.time() - t |
| 655 | | |
| 656 | | |
| 657 | | print '=' * 40 |
| 658 | | #no blinding |
| 659 | | t = time.time() |
| 660 | | message = 'serial '*5 |
| 661 | | #print 'cleartext ', message |
| 662 | | cypher = encrypt(message,pub) |
| 663 | | #print 'cyphertext: ',cypher |
| 664 | | #print 'decrypted', decrypt(cypher,priv) |
| 665 | | decrypt(cypher,priv) |
| 666 | | signed = sign(message,priv) |
| 667 | | #print 'signed', signed |
| 668 | | #print 'verified', message == verify(signed,pub) |
| 669 | | print time.time() - t |
| 670 | | |
| 671 | | |
| 672 | | print '=' * 40 |
| 673 | | times = [] |
| 674 | | #blinding |
| | 449 | import time |
| | 450 | t = time.time() |
| | 451 | #(pub,priv) = gen_pubpriv_keys(1024) |
| | 452 | (pub,priv) = (dummypub,dummypriv) |
| | 453 | print '=' * 40 |
| | 454 | times = [] |
| | 455 | #blinding |
| | 456 | #message = 'f'*65 |
| | 457 | message = 'c3ab8ff13720e8ad9047dd39466b3c8974e592c2fa383d4a3960714caef0c4f2' |
| | 458 | #print 'cleartext ', message |
| | 459 | unblinder = getUnblinder(pub['n']) |
| | 460 | blinder = pow(invMod(unblinder, pub['n']), pub['e'],pub['n']) |
| | 461 | times.append(time.time() - t) |
| | 462 | t = time.time() |
| | 463 | |
| | 464 | blinded = blind(message,blinder,pub) |
| | 465 | times.append(time.time() - t) |
| | 466 | t = time.time() |
| | 467 | |
| | 468 | signedblind = encrypt_int(blinded, priv['d'], priv['p']*priv['q']) |
| | 469 | times.append(time.time() - t) |
| | 470 | t = time.time() |
| | 471 | |
| | 472 | unblinded = (signedblind * unblinder) % pub['n'] |
| | 473 | times.append(time.time() - t) |
| | 474 | t = time.time() |
| | 475 | |
| | 476 | print 'verifyied', message == verify(unblinded,pub) |
| | 477 | times.append(time.time() - t) |
| | 478 | print sum(times) - times[2] |
| | 479 | |
| | 480 | if 0: |
| | 481 | #full |