[top] [TitleIndex] [WordIndex]

# 2008/480a/schedule/2008-04-18/worksheet

Math 480 -- 2008-04-18: Prime and Factorization system:sage

<h1>Prime Numbers</h1>

A <b>prime number</b> is a positive integer divisible by only itself and 1.

id=8|
for n in [1..20]:
print '%-4s%-7s'%(n, is_prime(n))
///
1   False
2   True
3   True
4   False
5   True
6   False
7   True
8   False
9   False
10  False
11  True
12  False
13  True
14  False
15  False
16  False
17  True
18  False
19  True
20  False


id=10|
@interact
def pr(max=(10..1000)):
v = prime_range(max)
print v
show(points([(x,0) for x in v], pointsize=10), xmin=0, ymin=-0.1, ymax=0.1, figsize=[10,2])
///

<html><div padding=6 id='div-interact-10'> <table width=800px height=400px bgcolor='#c5c5c5'
cellpadding=15><tr><td bgcolor='#f9f9f9' valign=top align=left><table><tr><td align=right><font color="black">max </font></td><td><div id='slider-max-10' class='ui-slider-1' style='padding:0px;margin:0px;'><span class='ui-slider-handle'></span></div><script>setTimeout(function() { $('#slider-max-10').slider({ stepping: 1, minValue: 0, maxValue: 990, startValue: 0, change: function () { var position = Math.ceil($('#slider-max-10').slider('value')); interact(10, "sage.server.notebook.interact.update(10, \"max\", 68, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals())"); }
});}, 1);</script></td></tr>
</table><div id='cell-interact-10'><?__SAGE__START>
<table border=0 bgcolor='#white' width=100% height=100%>
<tr><td bgcolor=white align=left valign=top><pre><?__SAGE__TEXT></pre></td></tr>
<tr><td  align=left valign=top><?__SAGE__HTML></td></tr>
</table><?__SAGE__END></div></td>
</tr></table></div>
</html>


Finding the next prime after a given prime is quite fast

id=11|
time next_prime(10^100, proof=False)
///
10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000267
CPU time: 0.01 s,  Wall time: 0.01 s


id=12|
time next_prime(10^100, proof=True)
///
10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000267
CPU time: 0.28 s,  Wall time: 0.33 s


id=14|
# Let's make a personal prime''
time next_prime(101010101010101010101010101010101010101010101010)
///
101010101010101010101010101010101010101010101039
Time: CPU 0.02 s, Wall: 0.02 s


<br><br>By the way, it was an open problem for a very long time to give a deterministic polynomial time algorithm for primality. This was done a few years ago finally. See <a href="http://en.wikipedia.org/wiki/AKS_primality_test">http://en.wikipedia.org/wiki/AKS_primality_test</a>.

<br><br>This is not a function in Sage though, since there are <i>better in practice</i> nondeterministic primality tests, especially for numbers that aren't too big.

<br><br>The distribution of prime numbers is a problem that has fascinated people for a very long time.

<br><br>But what the heck does <i>distribution</i> mean?! There are many possible interpretations.

<br><br>One way to measure the distribution is to simply count the primes up to : $$• \pi(x) = \#\{\text{primes } p : p \leq x\}.$$

id=17|
@interact
def primecount(m=(1..1000)):
n = 100*m
print "Primes up to %s"%n
show(plot(prime_pi, 1, n), xmin=0)
///

<html><div padding=6 id='div-interact-17'> <table width=800px height=400px bgcolor='#c5c5c5'
cellpadding=15><tr><td bgcolor='#f9f9f9' valign=top align=left><table><tr><td align=right><font color="black">m </font></td><td><div id='slider-m-17' class='ui-slider-1' style='padding:0px;margin:0px;'><span class='ui-slider-handle'></span></div><script>setTimeout(function() { $('#slider-m-17').slider({ stepping: 1, minValue: 0, maxValue: 999, startValue: 0, change: function () { var position = Math.ceil($('#slider-m-17').slider('value')); interact(17, "sage.server.notebook.interact.update(17, \"m\", 69, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals())"); }
});}, 1);</script></td></tr>
</table><div id='cell-interact-17'><?__SAGE__START>
<table border=0 bgcolor='#white' width=100% height=100%>
<tr><td bgcolor=white align=left valign=top><pre><?__SAGE__TEXT></pre></td></tr>
<tr><td  align=left valign=top><?__SAGE__HTML></td></tr>
</table><?__SAGE__END></div></td>
</tr></table></div>
</html>


<br><br> That curve above looks pretty smooth. It's natural to wonder what it is. Maybe it's just really complicated since there isn't a known easy way to find primes?!

<br><br> In fact, there is a <b>one hundred thousand dollar (100,000)</b> cash prize for finding a prime with at least ten million (10,000,000) decimal digits. See <a href="http://w2.eff.org/awards/coop-prime-rules.php">http://w2.eff.org/awards/coop-prime-rules.php</a>. (Note: a 50,000 dollar award was given in 2000 for a million digit prime.)

<br><br> Surprisingly $$• \pi(x) \sim x/\log(x)$$

id=21|
@interact
def primecount(m=(1..1000)):
n = 100*m; var('x')
print "Primes up to %s"%n
show(plot(prime_pi, 2, n) + plot(x/log(x), 4, n, color='red') , xmin=0, figsize=[9,3])
///

<html><div padding=6 id='div-interact-21'> <table width=800px height=400px bgcolor='#c5c5c5'
cellpadding=15><tr><td bgcolor='#f9f9f9' valign=top align=left><table><tr><td align=right><font color="black">m </font></td><td><div id='slider-m-21' class='ui-slider-1' style='padding:0px;margin:0px;'><span class='ui-slider-handle'></span></div><script>setTimeout(function() { $('#slider-m-21').slider({ stepping: 1, minValue: 0, maxValue: 999, startValue: 0, change: function () { var position = Math.ceil($('#slider-m-21').slider('value')); interact(21, "sage.server.notebook.interact.update(21, \"m\", 70, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals())"); }
});}, 1);</script></td></tr>
</table><div id='cell-interact-21'><?__SAGE__START>
<table border=0 bgcolor='#white' width=100% height=100%>
<tr><td bgcolor=white align=left valign=top><pre><?__SAGE__TEXT></pre></td></tr>
<tr><td  align=left valign=top><?__SAGE__HTML></td></tr>
</table><?__SAGE__END></div></td>
</tr></table></div>
</html>


<br><br> In fact it is a <b>theorem</b> (the prime number theorem) that $$• \lim_{x\to\infty} \pi(x) / (x/\log(x)) = 1.$$

<br><br> Actually isn't <i>that</i> good, though it looks superb in the picture above...

id=28|
prime_pi(10^7)
///
664579


id=29|
float(10^7/log(10^7))
///
620420.68843321688


<br><br> Let $$• \text{Li}(x) = \int_{2}^x \frac{1}{\log(t)} dt$$ This is just the area under the curve of the graph of , which is somehow a more precise way of <i>counting</i>.

id=35|
@interact
def primecount(m=(1..1000)):
n = 100*m; var('x')
print "Primes up to %s"%n
show(plot(prime_pi, 2, n) + plot(Li, 4, n, color='red') , xmin=0, figsize=[9,3])
///

<html><div padding=6 id='div-interact-35'> <table width=800px height=400px bgcolor='#c5c5c5'
cellpadding=15><tr><td bgcolor='#f9f9f9' valign=top align=left><table><tr><td align=right><font color="black">m </font></td><td><div id='slider-m-35' class='ui-slider-1' style='padding:0px;margin:0px;'><span class='ui-slider-handle'></span></div><script>setTimeout(function() { $('#slider-m-35').slider({ stepping: 1, minValue: 0, maxValue: 999, startValue: 0, change: function () { var position = Math.ceil($('#slider-m-35').slider('value')); interact(35, "sage.server.notebook.interact.update(35, \"m\", 71, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals())"); }
});}, 1);</script></td></tr>
</table><div id='cell-interact-35'><?__SAGE__START>
<table border=0 bgcolor='#white' width=100% height=100%>
<tr><td bgcolor=white align=left valign=top><pre><?__SAGE__TEXT></pre></td></tr>
<tr><td  align=left valign=top><?__SAGE__HTML></td></tr>
</table><?__SAGE__END></div></td>
</tr></table></div>
</html>


<br><br> This is <i>really good!</i>

id=34|
prime_pi(10^7)
///
664579


id=33|
Li(10^7)
///
664917.35989


<br><br> It is a <b>conjecture</b> (the <i>Riemann Hypothesis</i>!) that for all , we have $$• |\pi(x) - \text{Li}(x)| \leq \sqrt{x}\cdot \log(x).$$ That is <i>square root accurate</i>, which is amazingly good.

<br><br> Oh, by the way, you get a million dollars if you can prove the above inequality... <a href="http://www.claymath.org/millennium/">http://www.claymath.org/millennium/</a>

id=39|

///


id=38|

///


<br><br> <b>THEOREM:</b> Every positive integer can be written uniquely as product of primes. Proof: Surprisingly not trivial. Basically use the Euclidean algorithm to prove a special divisibility property, namely "if divides , then divides a or divides ". The rest is straightforward.

id=42|
import random
def ftree(rows, v, i, F):
if len(v) > 0: # add a row to g at the ith level.
rows.append(v)
w = []
for i in range(len(v)):
k, _, _ = v[i]
if k is None or is_prime(k):
w.append((None,None,None))
else:
d = random.choice(divisors(k)[1:-1])
w.append((d,k,i))
e = k//d
if e == 1:
w.append((None,None))
else:
w.append((e,k,i))
if len(w) > len(v):
ftree(rows, w, i+1, F)
def draw_ftree(rows,font):
g = Graphics()
for i in range(len(rows)):
cur = rows[i]
for j in range(len(cur)):
e, f, k = cur[j]
if not e is None:
if is_prime(e):
c = (1,0,0)
else:
c = (0,0,.4)
g += text(str(e), (j*2-len(cur),-i), fontsize=font, rgbcolor=c)
if not k is None and not f is None:
g += line([(j*2-len(cur),-i), ((k*2)-len(rows[i-1]),-i+1)],
alpha=0.5)
return g

@interact
def factor_tree(n=100, font=(10, (8..20)), redraw=['Redraw']):
n = Integer(n)
rows = []
v = [(n,None,0)]
ftree(rows, v, 0, factor(n))
show(draw_ftree(rows, font), axes=False)
///

<html><div padding=6 id='div-interact-42'> <table width=800px height=400px bgcolor='#c5c5c5'
cellpadding=15><tr><td bgcolor='#f9f9f9' valign=top align=left><table><tr><td align=right><font color="black">n </font></td><td><input type='text' value='100' width=200px onchange='interact(42, "sage.server.notebook.interact.update(42, \"n\", 72, sage.server.notebook.interact.standard_b64decode(\""+encode64(this.value)+"\"), globals())")'></input></td></tr>
<tr><td align=right><font color="black">font </font></td><td><div id='slider-font-42' class='ui-slider-1' style='padding:0px;margin:0px;'><span class='ui-slider-handle'></span></div><script>setTimeout(function() { $('#slider-font-42').slider({ stepping: 1, minValue: 0, maxValue: 12, startValue: 2, change: function () { var position = Math.ceil($('#slider-font-42').slider('value')); interact(42, "sage.server.notebook.interact.update(42, \"font\", 73, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals())"); }
});}, 1);</script></td></tr>
<tr><td align=right><font color="black">redraw </font></td><td><table style="border:1px solid #dfdfdf;background-color:#efefef">
<tr><td><button style='' value='0' onclick='interact(42, "sage.server.notebook.interact.update(42, \"redraw\", 74, sage.server.notebook.interact.standard_b64decode(\""+encode64(this.value)+"\"), globals())")'>Redraw</button>
</td></tr></table></td></tr>
</table><div id='cell-interact-42'><?__SAGE__START>
<table border=0 bgcolor='#white' width=100% height=100%>
<tr><td bgcolor=white align=left valign=top><pre><?__SAGE__TEXT></pre></td></tr>
<tr><td  align=left valign=top><?__SAGE__HTML></td></tr>
</table><?__SAGE__END></div></td>
</tr></table></div>
</html>


id=44|

///


id=48|
n = next_prime(10^28) * next_prime(10^29)
///


id=45|
time factor(n)
///
10000000000000000000000000331 * 100000000000000000000000000319
CPU time: 10.84 s,  Wall time: 11.67 s


Integer factorization in Sage:

<ol>

• <li> The <tt>factor</tt> command. This calls PARI's factor command, which implements the standard modern factoring algorithms including the elliptic curve method and the quadratic sieve. <li> The <tt>qsieve</tt> command. This calls Sage's (Bill Hart's) quadratic sieve implementation of the quadratic sieve. The sieve reduces factoring to an linear algebra problem. This implementation is better than the one in PARI. <li> The <tt>ecm</tt> command. This calls Sage's GMP-ECM (Paul Zimmerman), which is better than the ECM implementation in PARI.

</ol>

<br><br> The Sage factor command never calls qsieve or ecm. It should. Fixing this would be a <i>superb student project</i>.

id=47|
time qsieve(n)
///
([10000000000000000000000000331, 100000000000000000000000000319], '')
CPU time: 0.01 s,  Wall time: 6.83 s


id=43|
n = next_prime(10^14) * next_prime(10^50)
///


id=40|
time factor(n)
///
100000000000031 * 100000000000000000000000000000000000000000000000151
CPU time: 3.81 s,  Wall time: 3.95 s


id=50|
time ecm.factor(n)
///
[100000000000031, 100000000000000000000000000000000000000000000000151]
CPU time: 0.01 s,  Wall time: 0.45 s


<br><br> Parallel computation is extremely relevant fast integer factorization, though it's not used at all above. See the <tt>DistributedFactor</tt> command in Sage...

id=53|

///


id=51|

///


id=52|

///


id=49|

///


2013-05-11 18:32