414 - 2010-01-25: powers, primality testing

850 days ago by wstein

Writing a number in binary
def write_in_binary(m): digits = [] while m: if m%2: digits.append(1) else: digits.append(0) m = m//2 return digits 
       
write_in_binary(521) 
        
[1, 0, 0, 1, 0, 0, 0, 0, 0, 1]
[1, 0, 0, 1, 0, 0, 0, 0, 0, 1]
1 + 0*2 + 0*2^2 + 2^3 + 0*2^4 + 0*2^5 + 0*2^6 + 0*2^7 + 0*2^8 + 2^9 
        
521
521
Mod(7,100)^25 
        
7
7
 
       

Problem: Compute 3^{521} \pmod{10000}.
m = 521 n = 10000 a = 3 
       
m 
        
521
521
m.digits(2) 
        
[1, 0, 0, 1, 0, 0, 0, 0, 0, 1]
[1, 0, 0, 1, 0, 0, 0, 0, 0, 1]
m.str(2) 
        
'1000001001'
'1000001001'
m == 2^9 + 2^3 + 1 
        
True
True
abar = Mod(a,n) 
       
abar^(2^9 + 2^3 + 1) 
        
3203
3203
Mod(a,n)^m 
        
3203
3203
abar * ((abar^2)^2)^2 * (((((((((abar^2)^2)^2)^2)^2)^2)^2)^2)^2) 
        
3203
3203
 
       

Exponentiation Mod n Calculator
%auto @interact def _(a = 2, m = 4, n = 100): print "%s = (%s)_2 (in binary)"%(m, m.str(2)) print "%s^%s = %s (mod %s)"%(a, m, Mod(a,n)^m, n) 
       
a 
m 
n 

Click to the left again to hide and once more to show the dynamic interactive window
 
       

Our First Primality Test

If 2^{n-1} \equiv 1 \pmod{n} then n has a chance of being prime.
# Is 323 prime? Nope! Mod(2,323)^322 
        
157
157
# Is 341 prime? Maybe! But, not for sure. Mod(2, 341)^340 
        
1
1
# What about a big odd random number? # Nope, but this doesn't give any info about a factorization! set_random_seed(0) n = 2*ZZ.random_element(10^50)+1 print n time print Mod(2,n)^(n-1) 
        
172693837716124418522092982979366532296314630498967
79829587767292159558640084729953131077372528003242
Time: CPU 0.00 s, Wall: 0.00 s
172693837716124418522092982979366532296314630498967
79829587767292159558640084729953131077372528003242
Time: CPU 0.00 s, Wall: 0.00 s
time factor(n) 
        
3 * 449862907000167013519 * 127960344532295706107782898131
Time: CPU 1.70 s, Wall: 1.70 s
3 * 449862907000167013519 * 127960344532295706107782898131
Time: CPU 1.70 s, Wall: 1.70 s
%auto def maybe_prime(n): return 2.powermod(n-1,n) == 1 
       
maybe_prime(2011) 
        
True
True
%auto @interact def _(up_to=(100,110,..,30000)): fails = 0 v = [3..up_to] print "Fails for n =", for n in v: if maybe_prime(n) != is_prime(n): print n, fails += 1 print "\n\nPrimality test using 2 fails %s percent of the time"%(100*float(fails) / len(v)) 
       
up_to 

Click to the left again to hide and once more to show the dynamic interactive window
 
       

Idle Question:  Is there an elementary argument that the percentage of failures go to 0 as n\to\infty?
 
       

Finding a Charmichael Number
%auto def is_charmichael_number(n): if is_prime(n): return False for a in [2..n]: if gcd(a,n) == 1: if a.powermod(n-1,n) != 1: return False return True 
       
for n in [2..7000]: if is_charmichael_number(n): print "n = %s"%n 
        
n = 561
n = 1105
n = 1729
n = 2465
n = 2821
n = 6601
n = 561
n = 1105
n = 1729
n = 2465
n = 2821
n = 6601
12^3 + 1^3 
        
1729
1729
 
       

Miller-Rabin

An extremely powerful probabilistic primality test.  Every call increases the chances of a correct conclusion.
%auto def miller_rabin_round(n): k = valuation(n-1, 2) m = (n-1)//2^k a = ZZ.random_element(x=2,y=n-1) b = a.powermod(m,n) if b == 1 or b == n-1: return True for r in [1..k-1]: if b.powermod(2^r,n) == n-1: return True return False def miller_rabin(n, rounds=10): for i in range(rounds): if not miller_rabin_round(n): return False return True 
       
%auto @interact def _(nmax=1000, rounds=[1..5], retry=['Go!']): print "Failures up to %s: "%nmax, k = 0 for n in [5..nmax]: if miller_rabin(n,rounds) != is_prime(n): print n, k += 1 if k == 0: print "None" 
       
nmax 
rounds 
retry 

Click to the left again to hide and once more to show the dynamic interactive window