1#! /usr/bin/env python
2#
3# Provide some simple capabilities from number theory.
4#
5# Version of 2008.11.14.
6#
7# Written in 2005 and 2006 by Peter Pearson and placed in the public domain.
8# Revision history:
9#   2008.11.14: Use pow(base, exponent, modulus) for modular_exp.
10#               Make gcd and lcm accept arbitrarly many arguments.
11
12from __future__ import division
13
14from six import integer_types, PY2
15from six.moves import reduce
16
17try:
18    xrange
19except NameError:
20    xrange = range
21try:
22    from gmpy2 import powmod
23
24    GMPY2 = True
25    GMPY = False
26except ImportError:
27    GMPY2 = False
28    try:
29        from gmpy import mpz
30
31        GMPY = True
32    except ImportError:
33        GMPY = False
34
35import math
36import warnings
37
38
39class Error(Exception):
40    """Base class for exceptions in this module."""
41
42    pass
43
44
45class SquareRootError(Error):
46    pass
47
48
49class NegativeExponentError(Error):
50    pass
51
52
53def modular_exp(base, exponent, modulus):  # pragma: no cover
54    """Raise base to exponent, reducing by modulus"""
55    # deprecated in 0.14
56    warnings.warn(
57        "Function is unused in library code. If you use this code, "
58        "change to pow() builtin.",
59        DeprecationWarning,
60    )
61    if exponent < 0:
62        raise NegativeExponentError(
63            "Negative exponents (%d) not allowed" % exponent
64        )
65    return pow(base, exponent, modulus)
66
67
68def polynomial_reduce_mod(poly, polymod, p):
69    """Reduce poly by polymod, integer arithmetic modulo p.
70
71    Polynomials are represented as lists of coefficients
72    of increasing powers of x."""
73
74    # This module has been tested only by extensive use
75    # in calculating modular square roots.
76
77    # Just to make this easy, require a monic polynomial:
78    assert polymod[-1] == 1
79
80    assert len(polymod) > 1
81
82    while len(poly) >= len(polymod):
83        if poly[-1] != 0:
84            for i in xrange(2, len(polymod) + 1):
85                poly[-i] = (poly[-i] - poly[-1] * polymod[-i]) % p
86        poly = poly[0:-1]
87
88    return poly
89
90
91def polynomial_multiply_mod(m1, m2, polymod, p):
92    """Polynomial multiplication modulo a polynomial over ints mod p.
93
94    Polynomials are represented as lists of coefficients
95    of increasing powers of x."""
96
97    # This is just a seat-of-the-pants implementation.
98
99    # This module has been tested only by extensive use
100    # in calculating modular square roots.
101
102    # Initialize the product to zero:
103
104    prod = (len(m1) + len(m2) - 1) * [0]
105
106    # Add together all the cross-terms:
107
108    for i in xrange(len(m1)):
109        for j in xrange(len(m2)):
110            prod[i + j] = (prod[i + j] + m1[i] * m2[j]) % p
111
112    return polynomial_reduce_mod(prod, polymod, p)
113
114
115def polynomial_exp_mod(base, exponent, polymod, p):
116    """Polynomial exponentiation modulo a polynomial over ints mod p.
117
118    Polynomials are represented as lists of coefficients
119    of increasing powers of x."""
120
121    # Based on the Handbook of Applied Cryptography, algorithm 2.227.
122
123    # This module has been tested only by extensive use
124    # in calculating modular square roots.
125
126    assert exponent < p
127
128    if exponent == 0:
129        return [1]
130
131    G = base
132    k = exponent
133    if k % 2 == 1:
134        s = G
135    else:
136        s = [1]
137
138    while k > 1:
139        k = k // 2
140        G = polynomial_multiply_mod(G, G, polymod, p)
141        if k % 2 == 1:
142            s = polynomial_multiply_mod(G, s, polymod, p)
143
144    return s
145
146
147def jacobi(a, n):
148    """Jacobi symbol"""
149
150    # Based on the Handbook of Applied Cryptography (HAC), algorithm 2.149.
151
152    # This function has been tested by comparison with a small
153    # table printed in HAC, and by extensive use in calculating
154    # modular square roots.
155
156    assert n >= 3
157    assert n % 2 == 1
158    a = a % n
159    if a == 0:
160        return 0
161    if a == 1:
162        return 1
163    a1, e = a, 0
164    while a1 % 2 == 0:
165        a1, e = a1 // 2, e + 1
166    if e % 2 == 0 or n % 8 == 1 or n % 8 == 7:
167        s = 1
168    else:
169        s = -1
170    if a1 == 1:
171        return s
172    if n % 4 == 3 and a1 % 4 == 3:
173        s = -s
174    return s * jacobi(n % a1, a1)
175
176
177def square_root_mod_prime(a, p):
178    """Modular square root of a, mod p, p prime."""
179
180    # Based on the Handbook of Applied Cryptography, algorithms 3.34 to 3.39.
181
182    # This module has been tested for all values in [0,p-1] for
183    # every prime p from 3 to 1229.
184
185    assert 0 <= a < p
186    assert 1 < p
187
188    if a == 0:
189        return 0
190    if p == 2:
191        return a
192
193    jac = jacobi(a, p)
194    if jac == -1:
195        raise SquareRootError("%d has no square root modulo %d" % (a, p))
196
197    if p % 4 == 3:
198        return pow(a, (p + 1) // 4, p)
199
200    if p % 8 == 5:
201        d = pow(a, (p - 1) // 4, p)
202        if d == 1:
203            return pow(a, (p + 3) // 8, p)
204        if d == p - 1:
205            return (2 * a * pow(4 * a, (p - 5) // 8, p)) % p
206        raise RuntimeError("Shouldn't get here.")
207
208    if PY2:
209        # xrange on python2 can take integers representable as C long only
210        range_top = min(0x7FFFFFFF, p)
211    else:
212        range_top = p
213    for b in xrange(2, range_top):
214        if jacobi(b * b - 4 * a, p) == -1:
215            f = (a, -b, 1)
216            ff = polynomial_exp_mod((0, 1), (p + 1) // 2, f, p)
217            assert ff[1] == 0
218            return ff[0]
219    raise RuntimeError("No b found.")
220
221
222if GMPY2:
223
224    def inverse_mod(a, m):
225        """Inverse of a mod m."""
226        if a == 0:
227            return 0
228        return powmod(a, -1, m)
229
230
231elif GMPY:
232
233    def inverse_mod(a, m):
234        """Inverse of a mod m."""
235        # while libgmp likely does support inverses modulo, it is accessible
236        # only using the native `pow()` function, and `pow()` sanity checks
237        # the parameters before passing them on to underlying implementation
238        # on Python2
239        if a == 0:
240            return 0
241        a = mpz(a)
242        m = mpz(m)
243
244        lm, hm = mpz(1), mpz(0)
245        low, high = a % m, m
246        while low > 1:
247            r = high // low
248            lm, low, hm, high = hm - lm * r, high - low * r, lm, low
249
250        return lm % m
251
252
253else:
254
255    def inverse_mod(a, m):
256        """Inverse of a mod m."""
257
258        if a == 0:
259            return 0
260
261        lm, hm = 1, 0
262        low, high = a % m, m
263        while low > 1:
264            r = high // low
265            lm, low, hm, high = hm - lm * r, high - low * r, lm, low
266
267        return lm % m
268
269
270try:
271    gcd2 = math.gcd
272except AttributeError:
273
274    def gcd2(a, b):
275        """Greatest common divisor using Euclid's algorithm."""
276        while a:
277            a, b = b % a, a
278        return b
279
280
281def gcd(*a):
282    """Greatest common divisor.
283
284    Usage: gcd([ 2, 4, 6 ])
285    or:    gcd(2, 4, 6)
286    """
287
288    if len(a) > 1:
289        return reduce(gcd2, a)
290    if hasattr(a[0], "__iter__"):
291        return reduce(gcd2, a[0])
292    return a[0]
293
294
295def lcm2(a, b):
296    """Least common multiple of two integers."""
297
298    return (a * b) // gcd(a, b)
299
300
301def lcm(*a):
302    """Least common multiple.
303
304    Usage: lcm([ 3, 4, 5 ])
305    or:    lcm(3, 4, 5)
306    """
307
308    if len(a) > 1:
309        return reduce(lcm2, a)
310    if hasattr(a[0], "__iter__"):
311        return reduce(lcm2, a[0])
312    return a[0]
313
314
315def factorization(n):
316    """Decompose n into a list of (prime,exponent) pairs."""
317
318    assert isinstance(n, integer_types)
319
320    if n < 2:
321        return []
322
323    result = []
324    d = 2
325
326    # Test the small primes:
327
328    for d in smallprimes:
329        if d > n:
330            break
331        q, r = divmod(n, d)
332        if r == 0:
333            count = 1
334            while d <= n:
335                n = q
336                q, r = divmod(n, d)
337                if r != 0:
338                    break
339                count = count + 1
340            result.append((d, count))
341
342    # If n is still greater than the last of our small primes,
343    # it may require further work:
344
345    if n > smallprimes[-1]:
346        if is_prime(n):  # If what's left is prime, it's easy:
347            result.append((n, 1))
348        else:  # Ugh. Search stupidly for a divisor:
349            d = smallprimes[-1]
350            while 1:
351                d = d + 2  # Try the next divisor.
352                q, r = divmod(n, d)
353                if q < d:  # n < d*d means we're done, n = 1 or prime.
354                    break
355                if r == 0:  # d divides n. How many times?
356                    count = 1
357                    n = q
358                    while d <= n:  # As long as d might still divide n,
359                        q, r = divmod(n, d)  # see if it does.
360                        if r != 0:
361                            break
362                        n = q  # It does. Reduce n, increase count.
363                        count = count + 1
364                    result.append((d, count))
365            if n > 1:
366                result.append((n, 1))
367
368    return result
369
370
371def phi(n):  # pragma: no cover
372    """Return the Euler totient function of n."""
373    # deprecated in 0.14
374    warnings.warn(
375        "Function is unused by library code. If you use this code, "
376        "please open an issue in "
377        "https://github.com/warner/python-ecdsa",
378        DeprecationWarning,
379    )
380
381    assert isinstance(n, integer_types)
382
383    if n < 3:
384        return 1
385
386    result = 1
387    ff = factorization(n)
388    for f in ff:
389        e = f[1]
390        if e > 1:
391            result = result * f[0] ** (e - 1) * (f[0] - 1)
392        else:
393            result = result * (f[0] - 1)
394    return result
395
396
397def carmichael(n):  # pragma: no cover
398    """Return Carmichael function of n.
399
400    Carmichael(n) is the smallest integer x such that
401    m**x = 1 mod n for all m relatively prime to n.
402    """
403    # deprecated in 0.14
404    warnings.warn(
405        "Function is unused by library code. If you use this code, "
406        "please open an issue in "
407        "https://github.com/warner/python-ecdsa",
408        DeprecationWarning,
409    )
410
411    return carmichael_of_factorized(factorization(n))
412
413
414def carmichael_of_factorized(f_list):  # pragma: no cover
415    """Return the Carmichael function of a number that is
416    represented as a list of (prime,exponent) pairs.
417    """
418    # deprecated in 0.14
419    warnings.warn(
420        "Function is unused by library code. If you use this code, "
421        "please open an issue in "
422        "https://github.com/warner/python-ecdsa",
423        DeprecationWarning,
424    )
425
426    if len(f_list) < 1:
427        return 1
428
429    result = carmichael_of_ppower(f_list[0])
430    for i in xrange(1, len(f_list)):
431        result = lcm(result, carmichael_of_ppower(f_list[i]))
432
433    return result
434
435
436def carmichael_of_ppower(pp):  # pragma: no cover
437    """Carmichael function of the given power of the given prime."""
438    # deprecated in 0.14
439    warnings.warn(
440        "Function is unused by library code. If you use this code, "
441        "please open an issue in "
442        "https://github.com/warner/python-ecdsa",
443        DeprecationWarning,
444    )
445
446    p, a = pp
447    if p == 2 and a > 2:
448        return 2 ** (a - 2)
449    else:
450        return (p - 1) * p ** (a - 1)
451
452
453def order_mod(x, m):  # pragma: no cover
454    """Return the order of x in the multiplicative group mod m."""
455    # deprecated in 0.14
456    warnings.warn(
457        "Function is unused by library code. If you use this code, "
458        "please open an issue in "
459        "https://github.com/warner/python-ecdsa",
460        DeprecationWarning,
461    )
462
463    # Warning: this implementation is not very clever, and will
464    # take a long time if m is very large.
465
466    if m <= 1:
467        return 0
468
469    assert gcd(x, m) == 1
470
471    z = x
472    result = 1
473    while z != 1:
474        z = (z * x) % m
475        result = result + 1
476    return result
477
478
479def largest_factor_relatively_prime(a, b):  # pragma: no cover
480    """Return the largest factor of a relatively prime to b."""
481    # deprecated in 0.14
482    warnings.warn(
483        "Function is unused by library code. If you use this code, "
484        "please open an issue in "
485        "https://github.com/warner/python-ecdsa",
486        DeprecationWarning,
487    )
488
489    while 1:
490        d = gcd(a, b)
491        if d <= 1:
492            break
493        b = d
494        while 1:
495            q, r = divmod(a, d)
496            if r > 0:
497                break
498            a = q
499    return a
500
501
502def kinda_order_mod(x, m):  # pragma: no cover
503    """Return the order of x in the multiplicative group mod m',
504    where m' is the largest factor of m relatively prime to x.
505    """
506    # deprecated in 0.14
507    warnings.warn(
508        "Function is unused by library code. If you use this code, "
509        "please open an issue in "
510        "https://github.com/warner/python-ecdsa",
511        DeprecationWarning,
512    )
513
514    return order_mod(x, largest_factor_relatively_prime(m, x))
515
516
517def is_prime(n):
518    """Return True if x is prime, False otherwise.
519
520    We use the Miller-Rabin test, as given in Menezes et al. p. 138.
521    This test is not exact: there are composite values n for which
522    it returns True.
523
524    In testing the odd numbers from 10000001 to 19999999,
525    about 66 composites got past the first test,
526    5 got past the second test, and none got past the third.
527    Since factors of 2, 3, 5, 7, and 11 were detected during
528    preliminary screening, the number of numbers tested by
529    Miller-Rabin was (19999999 - 10000001)*(2/3)*(4/5)*(6/7)
530    = 4.57 million.
531    """
532
533    # (This is used to study the risk of false positives:)
534    global miller_rabin_test_count
535
536    miller_rabin_test_count = 0
537
538    if n <= smallprimes[-1]:
539        if n in smallprimes:
540            return True
541        else:
542            return False
543
544    if gcd(n, 2 * 3 * 5 * 7 * 11) != 1:
545        return False
546
547    # Choose a number of iterations sufficient to reduce the
548    # probability of accepting a composite below 2**-80
549    # (from Menezes et al. Table 4.4):
550
551    t = 40
552    n_bits = 1 + int(math.log(n, 2))
553    for k, tt in (
554        (100, 27),
555        (150, 18),
556        (200, 15),
557        (250, 12),
558        (300, 9),
559        (350, 8),
560        (400, 7),
561        (450, 6),
562        (550, 5),
563        (650, 4),
564        (850, 3),
565        (1300, 2),
566    ):
567        if n_bits < k:
568            break
569        t = tt
570
571    # Run the test t times:
572
573    s = 0
574    r = n - 1
575    while (r % 2) == 0:
576        s = s + 1
577        r = r // 2
578    for i in xrange(t):
579        a = smallprimes[i]
580        y = pow(a, r, n)
581        if y != 1 and y != n - 1:
582            j = 1
583            while j <= s - 1 and y != n - 1:
584                y = pow(y, 2, n)
585                if y == 1:
586                    miller_rabin_test_count = i + 1
587                    return False
588                j = j + 1
589            if y != n - 1:
590                miller_rabin_test_count = i + 1
591                return False
592    return True
593
594
595def next_prime(starting_value):
596    """Return the smallest prime larger than the starting value."""
597
598    if starting_value < 2:
599        return 2
600    result = (starting_value + 1) | 1
601    while not is_prime(result):
602        result = result + 2
603    return result
604
605
606smallprimes = [
607    2,
608    3,
609    5,
610    7,
611    11,
612    13,
613    17,
614    19,
615    23,
616    29,
617    31,
618    37,
619    41,
620    43,
621    47,
622    53,
623    59,
624    61,
625    67,
626    71,
627    73,
628    79,
629    83,
630    89,
631    97,
632    101,
633    103,
634    107,
635    109,
636    113,
637    127,
638    131,
639    137,
640    139,
641    149,
642    151,
643    157,
644    163,
645    167,
646    173,
647    179,
648    181,
649    191,
650    193,
651    197,
652    199,
653    211,
654    223,
655    227,
656    229,
657    233,
658    239,
659    241,
660    251,
661    257,
662    263,
663    269,
664    271,
665    277,
666    281,
667    283,
668    293,
669    307,
670    311,
671    313,
672    317,
673    331,
674    337,
675    347,
676    349,
677    353,
678    359,
679    367,
680    373,
681    379,
682    383,
683    389,
684    397,
685    401,
686    409,
687    419,
688    421,
689    431,
690    433,
691    439,
692    443,
693    449,
694    457,
695    461,
696    463,
697    467,
698    479,
699    487,
700    491,
701    499,
702    503,
703    509,
704    521,
705    523,
706    541,
707    547,
708    557,
709    563,
710    569,
711    571,
712    577,
713    587,
714    593,
715    599,
716    601,
717    607,
718    613,
719    617,
720    619,
721    631,
722    641,
723    643,
724    647,
725    653,
726    659,
727    661,
728    673,
729    677,
730    683,
731    691,
732    701,
733    709,
734    719,
735    727,
736    733,
737    739,
738    743,
739    751,
740    757,
741    761,
742    769,
743    773,
744    787,
745    797,
746    809,
747    811,
748    821,
749    823,
750    827,
751    829,
752    839,
753    853,
754    857,
755    859,
756    863,
757    877,
758    881,
759    883,
760    887,
761    907,
762    911,
763    919,
764    929,
765    937,
766    941,
767    947,
768    953,
769    967,
770    971,
771    977,
772    983,
773    991,
774    997,
775    1009,
776    1013,
777    1019,
778    1021,
779    1031,
780    1033,
781    1039,
782    1049,
783    1051,
784    1061,
785    1063,
786    1069,
787    1087,
788    1091,
789    1093,
790    1097,
791    1103,
792    1109,
793    1117,
794    1123,
795    1129,
796    1151,
797    1153,
798    1163,
799    1171,
800    1181,
801    1187,
802    1193,
803    1201,
804    1213,
805    1217,
806    1223,
807    1229,
808]
809
810miller_rabin_test_count = 0
811