Writing a number in binary
{{{id=46|
def write_in_binary(m):
digits = []
while m:
if m%2: digits.append(1)
else: digits.append(0)
m = m//2
return digits
///
}}}
{{{id=45|
write_in_binary(521)
///
[1, 0, 0, 1, 0, 0, 0, 0, 0, 1]
}}}
{{{id=47|
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
}}}
{{{id=48|
Mod(7,100)^25
///
7
}}}
{{{id=44|
///
}}}
Problem: Compute $3^{521} \pmod{10000}$.
{{{id=1|
m = 521
n = 10000
a = 3
///
}}}
{{{id=6|
m
///
521
}}}
{{{id=3|
m.digits(2)
///
[1, 0, 0, 1, 0, 0, 0, 0, 0, 1]
}}}
{{{id=8|
m.str(2)
///
'1000001001'
}}}
{{{id=9|
m == 2^9 + 2^3 + 1
///
True
}}}
{{{id=10|
abar = Mod(a,n)
///
}}}
{{{id=11|
abar^(2^9 + 2^3 + 1)
///
3203
}}}
{{{id=4|
Mod(a,n)^m
///
3203
}}}
{{{id=13|
abar * ((abar^2)^2)^2 * (((((((((abar^2)^2)^2)^2)^2)^2)^2)^2)^2)
///
3203
}}}
{{{id=17|
///
}}}
Exponentiation Mod n Calculator
{{{id=16|
%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)
///
}}}
{{{id=15|
///
}}}
Our First Primality Test
If $2^{n-1} \equiv 1 \pmod{n}$ then $n$ has a chance of being prime.
{{{id=33|
# Is 323 prime? Nope!
Mod(2,323)^322
///
157
}}}
{{{id=32|
# Is 341 prime? Maybe! But, not for sure.
Mod(2, 341)^340
///
1
}}}
{{{id=35|
# 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
}}}
{{{id=34|
time factor(n)
///
3 * 449862907000167013519 * 127960344532295706107782898131
Time: CPU 1.70 s, Wall: 1.70 s
}}}
{{{id=28|
%auto
def maybe_prime(n):
return 2.powermod(n-1,n) == 1
///
}}}
{{{id=31|
maybe_prime(2011)
///
True
}}}
{{{id=25|
%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))
///
}}}
{{{id=29|
///
}}}
Idle Question: Is there an elementary argument that the percentage of failures go to 0 as $n\to\infty$?
{{{id=24|
///
}}}
Finding a Charmichael Number
{{{id=14|
%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
///
}}}
{{{id=21|
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
^C
^C
}}}
{{{id=22|
12^3 + 1^3
///
1729
}}}
{{{id=36|
///
}}}
Miller-Rabin
An extremely powerful probabilistic primality test. Every call increases the chances of a correct conclusion.
{{{id=37|
%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
///
}}}
{{{id=23|
%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"
///
}}}
{{{id=40|
///
}}}
Lucas-Lehnmer test for Primality of Mersenne Numbers
{{{id=39|
%auto
def is_prime_lucas_lehmer(p):
s = Mod(4, 2^p - 1)
for i in range(3, p+1):
s = s^2 - 2
return s == 0
///
}}}
Find the Mersenne primes $2^p-1$ for $p<2000$.
{{{id=42|
%time
for p in prime_range(2000):
if is_prime_lucas_lehmer(p):
print p, len((2^p-1).digits())
///
3 1
5 2
7 3
13 4
17 6
19 6
31 10
61 19
89 27
107 33
127 39
521 157
607 183
1279 386
CPU time: 3.41 s, Wall time: 3.47 s
}}}
{{{id=43|
///
}}}