1417-Writing a number in binary.sagews

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*2 + 0*2^2 + 2^3 + 0*2^4 + 0*2^5 + 0*2^6 + 0*2^7 + 0*2^8 + 2^9

521

Mod(7,100)^25

7

Problem: Compute $3^{521} \pmod{10000}$.

m = 521
n = 10000
a = 3

m

521

m.digits(2)

[1, 0, 0, 1, 0, 0, 0, 0, 0, 1]

m.str(2)

'1000001001'

m == 2^9 + 2^3 + 1

True

abar = Mod(a,n)

abar^(2^9 + 2^3 + 1)

3203

Mod(a,n)^m

3203

abar * ((abar^2)^2)^2 * (((((((((abar^2)^2)^2)^2)^2)^2)^2)^2)^2)

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  [removed] [removed] [removed] [removed]

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

# Is 341 prime?  Maybe!  But, not for sure.
Mod(2, 341)^340

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

time factor(n)

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

%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  [removed] [removed] [removed] [removed]

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

12^3 + 1^3

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  [removed] [removed] [removed] [removed]