| 1 | ###################################################################################### |
|---|
| 2 | # NUMBTHY.PY |
|---|
| 3 | # Basic Number Theory functions implemented in Python |
|---|
| 4 | # Note: Currently requires Python 2.x (uses +=, %= and other 2.x-isms) |
|---|
| 5 | # Note: Currently requires Python 2.3 (uses implicit long integers - could be back-ported though) |
|---|
| 6 | # Author: Robert Campbell, <campbell@math.umbc.edu> |
|---|
| 7 | # Date: 27 April, 2007 |
|---|
| 8 | # Version 0.41 |
|---|
| 9 | ###################################################################################### |
|---|
| 10 | """Basic number theory functions. |
|---|
| 11 | Functions implemented are: |
|---|
| 12 | gcd(a,b) - Compute the greatest common divisor of a and b. |
|---|
| 13 | xgcd(a,b) - Find [g,x,y] such that g=gcd(a,b) and g = ax + by. |
|---|
| 14 | powmod(b,e,n) - Compute b^e mod n efficiently. |
|---|
| 15 | isprime(n) - Test whether n is prime using a variety of pseudoprime tests. |
|---|
| 16 | isprimeF(n,b) - Test whether n is prime or a Fermat pseudoprime to base b. |
|---|
| 17 | isprimeE(n,b) - Test whether n is prime or an Euler pseudoprime to base b. |
|---|
| 18 | eulerphi(n) - Computer Euler's Phi function of n - the number of integers strictly less than n which are coprime to n. |
|---|
| 19 | carmichaellambda(n) - Computer Carmichael's Lambda function of n - the smallest exponent e such that b**e = 1 for all b coprime to n. |
|---|
| 20 | factor(n) - Find a factor of n using a variety of methods. |
|---|
| 21 | factors(n) - Return a sorted list of the prime factors of n. |
|---|
| 22 | factorPR(n) - Find a factor of n using the Pollard Rho method |
|---|
| 23 | isprimitive(g,n) - Test whether g is primitive - generates the group of units mod n. |
|---|
| 24 | A list of the functions implemented in numbthy is printed by the command info().""" |
|---|
| 25 | |
|---|
| 26 | |
|---|
| 27 | Version = 'NUMBTHY.PY, version 0.1, 10 Jan, 2003, by Robert Campbell, campbell@math.umbc.edu' |
|---|
| 28 | |
|---|
| 29 | import math # Use sqrt, floor |
|---|
| 30 | |
|---|
| 31 | def gcd(a,b): |
|---|
| 32 | """gcd(a,b) returns the greatest common divisor of the integers a and b.""" |
|---|
| 33 | if a == 0: |
|---|
| 34 | return b |
|---|
| 35 | return abs(gcd(b % a, a)) |
|---|
| 36 | |
|---|
| 37 | def powmod(b,e,n): |
|---|
| 38 | """powmod(b,e,n) computes the eth power of b mod n. |
|---|
| 39 | (Actually, this is not needed, as pow(b,e,n) does the same thing for positive integers. |
|---|
| 40 | This will be useful in future for non-integers or inverses.""" |
|---|
| 41 | accum = 1; i = 0; bpow2 = b |
|---|
| 42 | while ((e>>i)>0): |
|---|
| 43 | if((e>>i) & 1): |
|---|
| 44 | accum = (accum*bpow2) % n |
|---|
| 45 | bpow2 = (bpow2*bpow2) % n |
|---|
| 46 | i+=1 |
|---|
| 47 | return accum |
|---|
| 48 | |
|---|
| 49 | def xgcd(a,b): |
|---|
| 50 | """xgcd(a,b) returns a list of form [g,x,y], where g is gcd(a,b) and |
|---|
| 51 | x,y satisfy the equation g = ax + by.""" |
|---|
| 52 | a1=1; b1=0; a2=0; b2=1; aneg=1; bneg=1; swap = False |
|---|
| 53 | if(a < 0): |
|---|
| 54 | a = -a; aneg=-1 |
|---|
| 55 | if(b < 0): |
|---|
| 56 | b = -b; bneg=-1 |
|---|
| 57 | if(b > a): |
|---|
| 58 | swap = True |
|---|
| 59 | [a,b] = [b,a] |
|---|
| 60 | while (1): |
|---|
| 61 | quot = -(a / b) |
|---|
| 62 | a = a % b |
|---|
| 63 | a1 = a1 + quot*a2; b1 = b1 + quot*b2 |
|---|
| 64 | if(a == 0): |
|---|
| 65 | if(swap): |
|---|
| 66 | return [b, b2*bneg, a2*aneg] |
|---|
| 67 | else: |
|---|
| 68 | return [b, a2*aneg, b2*bneg] |
|---|
| 69 | quot = -(b / a) |
|---|
| 70 | b = b % a; |
|---|
| 71 | a2 = a2 + quot*a1; b2 = b2 + quot*b1 |
|---|
| 72 | if(b == 0): |
|---|
| 73 | if(swap): |
|---|
| 74 | return [a, b1*bneg, a1*aneg] |
|---|
| 75 | else: |
|---|
| 76 | return [a, a1*aneg, b1*bneg] |
|---|
| 77 | |
|---|
| 78 | def isprime(n): |
|---|
| 79 | """isprime(n) - Test whether n is prime using a variety of pseudoprime tests.""" |
|---|
| 80 | if (n in [2,3,5,7,11,13,17,19,23,29]): return True |
|---|
| 81 | return isprimeE(n,2) and isprimeE(n,3) and isprimeE(n,5) |
|---|
| 82 | |
|---|
| 83 | def isprimeF(n,b): |
|---|
| 84 | """isprimeF(n) - Test whether n is prime or a Fermat pseudoprime to base b.""" |
|---|
| 85 | return (pow(b,n-1,n) == 1) |
|---|
| 86 | |
|---|
| 87 | def isprimeE(n,b): |
|---|
| 88 | """isprimeE(n) - Test whether n is prime or an Euler pseudoprime to base b.""" |
|---|
| 89 | if (not isprimeF(n,b)): return False |
|---|
| 90 | r = n-1 |
|---|
| 91 | while (r % 2 == 0): r /= 2 |
|---|
| 92 | c = pow(b,r,n) |
|---|
| 93 | if (c == 1): return True |
|---|
| 94 | while (1): |
|---|
| 95 | if (c == 1): return False |
|---|
| 96 | if (c == n-1): return True |
|---|
| 97 | c = pow(c,2,n) |
|---|
| 98 | |
|---|
| 99 | def factor(n): |
|---|
| 100 | """factor(n) - Find a prime factor of n using a variety of methods.""" |
|---|
| 101 | if (isprime(n)): return n |
|---|
| 102 | for fact in [2,3,5,7,11,13,17,19,23,29]: |
|---|
| 103 | if n%fact == 0: return fact |
|---|
| 104 | return factorPR(n) # Needs work - no guarantee that a prime factor will be returned |
|---|
| 105 | |
|---|
| 106 | def factors(n): |
|---|
| 107 | """factors(n) - Return a sorted list of the prime factors of n.""" |
|---|
| 108 | if (isprime(n)): |
|---|
| 109 | return [n] |
|---|
| 110 | fact = factor(n) |
|---|
| 111 | if (fact == 1): return "Unable to factor "+str(n) |
|---|
| 112 | facts = factors(n/fact) + factors(fact) |
|---|
| 113 | facts.sort() |
|---|
| 114 | return facts |
|---|
| 115 | |
|---|
| 116 | def factorPR(n): |
|---|
| 117 | """factorPR(n) - Find a factor of n using the Pollard Rho method. |
|---|
| 118 | Note: This method will occasionally fail.""" |
|---|
| 119 | for slow in [2,3,4,6]: |
|---|
| 120 | numsteps=2*math.floor(math.sqrt(math.sqrt(n))); fast=slow; i=1 |
|---|
| 121 | while i<numsteps: |
|---|
| 122 | slow = (slow*slow + 1) % n |
|---|
| 123 | i = i + 1 |
|---|
| 124 | fast = (fast*fast + 1) % n |
|---|
| 125 | fast = (fast*fast + 1) % n |
|---|
| 126 | g = gcd(fast-slow,n) |
|---|
| 127 | if (g != 1): |
|---|
| 128 | if (g == n): |
|---|
| 129 | break |
|---|
| 130 | else: |
|---|
| 131 | return g |
|---|
| 132 | return 1 |
|---|
| 133 | |
|---|
| 134 | def eulerphi(n): |
|---|
| 135 | """eulerphi(n) - Computer Euler's Phi function of n - the number of integers |
|---|
| 136 | strictly less than n which are coprime to n. Otherwise defined as the order |
|---|
| 137 | of the group of integers mod n.""" |
|---|
| 138 | thefactors = factors(n) |
|---|
| 139 | thefactors.sort() |
|---|
| 140 | phi = 1 |
|---|
| 141 | oldfact = 1 |
|---|
| 142 | for fact in thefactors: |
|---|
| 143 | if fact==oldfact: |
|---|
| 144 | phi = phi*fact |
|---|
| 145 | else: |
|---|
| 146 | phi = phi*(fact-1) |
|---|
| 147 | oldfact = fact |
|---|
| 148 | return phi |
|---|
| 149 | |
|---|
| 150 | def carmichaellambda(n): |
|---|
| 151 | """carmichaellambda(n) - Computer Carmichael's Lambda function |
|---|
| 152 | of n - the smallest exponent e such that b**e = 1 for all b coprime to n. |
|---|
| 153 | Otherwise defined as the exponent of the group of integers mod n.""" |
|---|
| 154 | thefactors = factors(n) |
|---|
| 155 | thefactors.sort() |
|---|
| 156 | thefactors += [0] # Mark the end of the list of factors |
|---|
| 157 | carlambda = 1 # The Carmichael Lambda function of n |
|---|
| 158 | carlambda_comp = 1 # The Carmichael Lambda function of the component p**e |
|---|
| 159 | oldfact = 1 |
|---|
| 160 | for fact in thefactors: |
|---|
| 161 | if fact==oldfact: |
|---|
| 162 | carlambda_comp = (carlambda_comp*fact) |
|---|
| 163 | else: |
|---|
| 164 | if ((oldfact == 2) and (carlambda_comp >= 4)): carlambda_comp /= 2 # Z_(2**e) is not cyclic for e>=3 |
|---|
| 165 | if carlambda == 1: |
|---|
| 166 | carlambda = carlambda_comp |
|---|
| 167 | else: |
|---|
| 168 | carlambda = (carlambda * carlambda_comp)/gcd(carlambda,carlambda_comp) |
|---|
| 169 | carlambda_comp = fact-1 |
|---|
| 170 | oldfact = fact |
|---|
| 171 | return carlambda |
|---|
| 172 | |
|---|
| 173 | def isprimitive(g,n): |
|---|
| 174 | """isprimitive(g,n) - Test whether g is primitive - generates the group of units mod n.""" |
|---|
| 175 | if gcd(g,n) != 1: return False # Not in the group of units |
|---|
| 176 | order = eulerphi(n) |
|---|
| 177 | if carmichaellambda(n) != order: return False # Group of units isn't cyclic |
|---|
| 178 | orderfacts = factors(order) |
|---|
| 179 | oldfact = 1 |
|---|
| 180 | for fact in orderfacts: |
|---|
| 181 | if fact!=oldfact: |
|---|
| 182 | if pow(g,order/fact,n) == 1: return False |
|---|
| 183 | oldfact = fact |
|---|
| 184 | return True |
|---|
| 185 | |
|---|
| 186 | def info(): |
|---|
| 187 | """Return information about the module""" |
|---|
| 188 | print locals() |
|---|
| 189 | |
|---|
| 190 | # import numbthy |
|---|
| 191 | # reload(numbthy) |
|---|
| 192 | # print numbthy.__doc__ |
|---|
| 193 | |
|---|
