[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 $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 $$

$$ 

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 $$

$$

<br><br> Actually $x/\log(x)$ 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 $$

$$ This is just the area under the curve of the graph of $1/log(t)$, 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 $x\geq 2.01$, we have $$

$$ 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 $p$ divides $a*b$, then $p$ divides a or $p$ divides $b$". 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>

</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