Sage for Number Theory : A Primes Primer
system:sage

<h1><span style="font-family: arial,helvetica,sans-serif;">Sage Number Theory: A Primes Primer</span></h1>
<p><span style="font-family: arial,helvetica,sans-serif;">Erik Jacobson</span></p>
<p><span style="font-family: arial,helvetica,sans-serif;">Sage Days 13, March 1, 2009, Version 0.1</span></p>
<h2><span style="font-family: arial,helvetica,sans-serif;">Introduction<br /></span></h2>
<p><span style="font-family: arial,helvetica,sans-serif;">This primer is not intended to teach Sage or number theory but instead to provide examples introducing some of what is possible with Sage used for number theory.&nbsp; It is aimed at secondary math teachers and undergraduates who want a functional knowledge of Sage without the overhead of deep understanding.&nbsp; The documentation and tutorials at sagemath.org are a good place to look next; those interested in the programming side of things might check out python.org, the website for the language in which much of Sage is written.</span></p>
<p><span style="font-family: arial,helvetica,sans-serif;">This primer is divided into two sections: </span><span style="font-family: arial,helvetica,sans-serif;">"primality, seive of Eratosthenes, and the Euclidean algorithm" and </span><span style="font-family: arial,helvetica,sans-serif;">"Counting Primes, Plotting Primes, and Finding Sequences of Primes." The primality section explores modular division, factorization, and primality testing in Sage as well as implementing the seive of Eratosthenes and the Euclidean algorithm.&nbsp; The next section demonstrates the prime_pi() function and 2D-plotting functions in Sage, and provides examples exploring prime density and arithmetic sequences of primes.<br /></span></p>
<h2><span style="font-family: arial,helvetica,sans-serif;">Primality, Seive of Eratosthenes, and the Euclidean Algorithm<br /></span></h2>
<p><span style="font-family: arial,helvetica,sans-serif;">To use a Sage Notebook, type an expression you want to evaluate into the cell and press <span style="font-family: courier new,courier;">&lt;SHIFT&gt;&lt;ENTER&gt;</span> or click the <span style="font-family: courier new,courier;"><span style="text-decoration: underline;">evaluate</span></span> button.&nbsp; This worksheet is interactive--try evaluating the cell below; you can modify it if you wish.</span></p>
<p><span style="font-family: arial,helvetica,sans-serif;">Modular division is implemented in two ways in Sage: using the binary <span style="font-family: courier new,courier;">%</span> operator or the <span style="font-family: courier new,courier;">.mod()</span> method. </span><span style="font-family: arial,helvetica,sans-serif;">To display the value of an expresion, use the <span style="font-family: courier new,courier;">print</span> command followed by the expressions seperated by commas.&nbsp; If you are evaluating a single expression, the </span><span style="font-family: arial,helvetica,sans-serif;"><span style="font-family: courier new,courier;">print</span></span><span style="font-family: arial,helvetica,sans-serif;"> command is not necessary.</span></p>

{{{id=19|
print 24 % 5, 24.mod(5)
print 24 % 7, 24 % 2
///
}}}

<h3>Factoring<br /></h3>
<p>To factor a number in Sage, use the function factor() which can be called as a method of an integer or be passed an integer parameter.</p>

{{{id=22|
print 1278, " factors as ", 1278.factor()
print 1242, " factors as ", factor(1242)
print 59, "factors as ", factor(59)
///

1278  factors as  2 * 3^2 * 71
1242  factors as  2 * 3^3 * 23
59 factors as  59
}}}

<p>In the integers, a number is prime if and only if it is irreducible (it has only trivial factors: 1 and -1).&nbsp; The next example checks the first 100 numbers for primality using the function factor().&nbsp; The statement factor(i)[0][0] returns the smallest (and possibly only) factor of of i.</p>

{{{id=20|
for i in range(2,53):
    if i == factor(i)[0][0]:
        print i, " is prime."
///

2  is prime.
3  is prime.
5  is prime.
7  is prime.
11  is prime.
13  is prime.
17  is prime.
19  is prime.
23  is prime.
29  is prime.
31  is prime.
37  is prime.
41  is prime.
43  is prime.
47  is prime.
}}}

<h3>Sieve of Eratosthenes</h3>
<p>Factoring each number takes a long time and an ancient (and faster) way to find the primes between 2 and 53 is the Sieve of Eratosthenes, attrbuted to an ancient Greek living around 200 BC.&nbsp; The idea is to start with all the numbers between 2 and 53 and one by one, remove the multiples of each prime.</p>

<p><span style="font-family: arial,helvetica,sans-serif;"><br /> </span></p>

{{{id=34|
def seive(a_known_prime, possible_primes):
    print possible_primes
    checked_primes = possible_primes
    for i in possible_primes:     
      if (i%a_known_prime == 0) and (i != a_known_prime):
          checked_primes.remove(i)
    print
    print "Remove mulitples of ", a_known_prime
    print
    possible_primes = checked_primes
    if possible_primes.index(a_known_prime)<len(possible_primes)-1:
        return seive(possible_primes[possible_primes.index(a_known_prime)+1], possible_primes)
    else:
        return possible_primes
          
possible_primes = list(range(2, 53+1))
seive(2,possible_primes)
///
}}}

<p>Sage has a built in function that checks the primality of an integer called is_prime().&nbsp; The previous example can be rewritten with this code as follows.</p>

{{{id=31|
for i in range(2,53):
    if is_prime(i):
        print i, " is prime."
///

[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53]

Remove mulitples of  2

[2, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 43, 45, 47, 49, 51, 53]

Remove mulitples of  3

[2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, 41, 43, 47, 49, 53]

Remove mulitples of  5

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53]

Remove mulitples of  7

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  11

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  13

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  17

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  19

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  23

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  29

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  31

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  37

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  41

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  43

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  47

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]

Remove mulitples of  53

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53]
}}}

<h3><span style="font-family: arial,helvetica,sans-serif;">Euclidean Algorithm</span></h3>
<p><span style="font-family: arial,helvetica,sans-serif;">The Euclidean Algorithm finds the greatest common factor of two integers by dividing successively dividing each divisor by the remainder of the previous division.&nbsp; The example below implements the Euclidean algorithm using recursion, a function which calls itself.</span></p>

{{{id=37|
def EuclidAlg(a,b):
    r = a % b
    if r == 0:
        return
    else:
        EuclidAlg(b,r)
    print r,"=",a,"-",b,"*",(a-r)/b

print "EuclidAlg(121,105):"
EuclidAlg(121,105)
print "\n","EuclidAlg(67890,12345):"
EuclidAlg(67890,12345)
///

EuclidAlg(121,105):
1 = 7 - 2 * 3
2 = 9 - 7 * 1
7 = 16 - 9 * 1
9 = 105 - 16 * 6
16 = 121 - 105 * 1

EuclidAlg(67890,12345):
15 = 12345 - 6165 * 2
6165 = 67890 - 12345 * 5
}}}

<p>Sage has the built-in function gcd() that finds the greatest common factor of two or more numbers.</p>

{{{id=38|
print gcd(121,105)
print "\n",gcd(67890,12345)
///

1

15
}}}

<h2><span style="font-family: arial,helvetica,sans-serif;">Counting Primes, Plotting Primes, and Finding Sequences of Primes</span></h2>
<h3>Ploting the Prime Number Theorem</h3>
<p>That there are an infinite number of primes was known to ancient Greeks.&nbsp; The function primes_pi() in Sage returns the number of primes below a given parameter.&nbsp; The example below is an interactive plot comparing primes_pi() and its asymtotic bound x/log(x)-1.&nbsp; The @interact command allows you to create dynamic plots.&nbsp; After @interact, define a function which takes the user input.&nbsp; In this case we define a slider N, initialized to 100 that takes on values between 2 and 2000.&nbsp; The next line is html and latex expressing the title of the plot in mathematical type.&nbsp; The next line shows the plot of the prime_pi() function superimposed on a plot of x/log(x) -1 in blue.&nbsp; Both are graphed on the domain bounded above by N.</p>
<p>This example due to William Stein.</p>

{{{id=36|

///
}}}

{{{id=16|
@interact
def _(N=(100,(2..2000))):
    html("<font color='red'>$\pi(x)$</font> and <font color='blue'>$x/(\log(x)-1)$</font> for $x < %s$"%N)
    show(plot(prime_pi, 0, N, rgbcolor='red') + plot(x/(log(x)-1), 5, N, rgbcolor='blue'))
///

<html><!--notruncate--><div padding=6 id='div-interact-16'> <table width=800px height=20px bgcolor='#c5c5c5'
                 cellpadding=15><tr><td bgcolor='#f9f9f9' valign=top align=left><table><tr><td align=right><font color="black">N&nbsp;</font></td><td><table><tr><td>
    	<div id='slider-N-16' class='ui-slider ui-slider-3' style='margin:0px;'><span class='ui-slider-handle'></span></div>
    	</td><td><font color='black' id='slider-N-16-lbl'></font></td></tr></table><script>(function(){ var values = ["2","6","10","14","18","22","26","30","34","38","42","46","50","54","58","62","66","70","74","78","82","86","90","94","98","102","106","110","114","118","122","126","130","134","138","142","146","150","154","158","162","166","170","174","178","182","186","190","194","198","202","206","210","214","218","222","226","230","234","238","242","246","250","254","258","262","266","270","274","278","282","286","290","294","298","302","306","310","314","318","322","326","330","334","338","342","346","350","354","358","362","366","370","374","378","382","386","390","394","398","402","406","410","414","418","422","426","430","434","438","442","446","450","454","458","462","466","470","474","478","482","486","490","494","498","502","506","510","514","518","522","526","530","534","538","542","546","550","554","558","562","566","570","574","578","582","586","590","594","598","602","606","610","614","618","622","626","630","634","638","642","646","650","654","658","662","666","670","674","678","682","686","690","694","698","702","706","710","714","718","722","726","730","734","738","742","746","750","754","758","762","766","770","774","778","782","786","790","794","798","802","806","810","814","818","822","826","830","834","838","842","846","850","854","858","862","866","870","874","878","882","886","890","894","898","902","906","910","914","918","922","926","930","934","938","942","946","950","954","958","962","966","970","974","978","982","986","990","994","998","1003","1007","1011","1015","1019","1023","1027","1031","1035","1039","1043","1047","1051","1055","1059","1063","1067","1071","1075","1079","1083","1087","1091","1095","1099","1103","1107","1111","1115","1119","1123","1127","1131","1135","1139","1143","1147","1151","1155","1159","1163","1167","1171","1175","1179","1183","1187","1191","1195","1199","1203","1207","1211","1215","1219","1223","1227","1231","1235","1239","1243","1247","1251","1255","1259","1263","1267","1271","1275","1279","1283","1287","1291","1295","1299","1303","1307","1311","1315","1319","1323","1327","1331","1335","1339","1343","1347","1351","1355","1359","1363","1367","1371","1375","1379","1383","1387","1391","1395","1399","1403","1407","1411","1415","1419","1423","1427","1431","1435","1439","1443","1447","1451","1455","1459","1463","1467","1471","1475","1479","1483","1487","1491","1495","1499","1503","1507","1511","1515","1519","1523","1527","1531","1535","1539","1543","1547","1551","1555","1559","1563","1567","1571","1575","1579","1583","1587","1591","1595","1599","1603","1607","1611","1615","1619","1623","1627","1631","1635","1639","1643","1647","1651","1655","1659","1663","1667","1671","1675","1679","1683","1687","1691","1695","1699","1703","1707","1711","1715","1719","1723","1727","1731","1735","1739","1743","1747","1751","1755","1759","1763","1767","1771","1775","1779","1783","1787","1791","1795","1799","1803","1807","1811","1815","1819","1823","1827","1831","1835","1839","1843","1847","1851","1855","1859","1863","1867","1871","1875","1879","1883","1887","1891","1895","1899","1903","1907","1911","1915","1919","1923","1927","1931","1935","1939","1943","1947","1951","1955","1959","1963","1967","1971","1975","1979","1983","1987","1991","1995","2000"]; setTimeout(function() {
    $('#slider-N-16').slider({
    	stepping: 1, min: 0, max: 499, startValue: 24,
    	change: function (e,ui) { var position = ui.value; if(values!=null) $('#slider-N-16-lbl').text(values[position]); interact(16, "sage.server.notebook.interact.update(16, \"N\", 1, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals());sage.server.notebook.interact.recompute(16)"); },
    	slide: function(e,ui) { if(values!=null) $('#slider-N-16-lbl').text(values[ui.value]); }
    });
    if(values != null) $('#slider-N-16-lbl').text(values[$('#slider-N-16').slider('value')]);
    }, 1); })();</script></td></tr>
</table><div id='cell-interact-16'><?__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>
}}}

<h3>Plotting the density of primes</h3>
<p>Sage uses the matplotlib (http://matplotlib.sourceforge.net) library for its plotting needs and if one requires more control over plotting than the plot() function provides, the capabilities of matplotlib can be used directly.&nbsp; While a complete explanation of how matplotlib works is beyond the scope of this primer, this examples illustrates some functionality of matplotlib. Here we wish to plot the density of primes over the interval (0, 10^20).&nbsp; This is another perspective on prime_pi() function plotted above.</p>

{{{id=9|
x = map(lambda x:x*10^7+10^7,range(20))
y = map(lambda j: prime_pi(x[j])-prime_pi(x[j]-10^7+1), range(len(x)))

from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.figure import Figure
from matplotlib.ticker import *

fig = Figure()
canvas = FigureCanvasAgg(fig)
ax = fig.add_subplot(111)
ax.xaxis.set_major_formatter( FormatStrFormatter( '%d' ))
ax.xaxis.set_major_locator( MaxNLocator(5) )

ax.yaxis.set_major_formatter( FormatStrFormatter( '%d' ))
ax.yaxis.set_major_locator( MaxNLocator(20) )
ax.yaxis.grid(True, linestyle='-', which='minor')

ax.grid(True, linestyle='-', linewidth=.5)
ax.set_title('Number of Primes in intervals of width 10^7 from 0 to 10^20')
ax.set_xlabel('Upper bound of intervals of width 10^7')
ax.set_ylabel('Number of Primes')
ax.plot(x,y, 'go-', linewidth=1.0 )
canvas.print_figure('ex1_linear.png')
///
}}}

<p>The code above stores an 20 term arithmetic sequence of integers with a constant difference of 10^7 in the list x, the x-coordinates.&nbsp; Then the number of primes in each interval of width 10^7 is stored in the list y, the y coordinates. The axis lables are set with the by the ax.set_xlabel() and ax.set_ylabel() commands. The title is set with ax.set_title().&nbsp; The number of grid lines is set dependent on the lists x and y using MaxNLocator(number) passed to ax.xaxis.set_major_locator() and ax.yaxis.set_major_locator(), respectively.</p>
<h3>Arithmetic Sequences of Primes<br /></h3>
<p>Arithmetic sequences of primes are common but enumerating all of them is non-trivial.&nbsp; Shorter sequences are (luckily) easy to find. The idea is to pick a random prime and search for arithemetic sequences starting at this prime with differences between 1 and say, 500.&nbsp; The code below stores the length of the desired sequence in seq_length.&nbsp; P stores a list of primes between 0 and 10^5. &nbsp; This procedure is unlikely to produce results for each guess, so try 500 times with the for loop indexed by i.&nbsp; the variable j stores a random index into the list of primes, P.&nbsp; For each possible difference, scale the list from 1 to seq_lenth by k and shift by the randomly chosen prime P[j], and check to see if all the numbers in this list are primes.&nbsp; Then print the initial prime, the difference, and the sequence of primes.</p>

{{{id=7|
seq_length = 5
P=prime_range(10^4)
for i in range(1,500):
    j = int(random()*len(P))
    for k in range(1,500):
        arithmetic_sequence = map(lambda x: x*k+P[j], range(1, seq_length+1))
        if set(arithmetic_sequence).issubset(set(P)):
            print P[j], k, arithmetic_sequence
///

331 240 [571, 811, 1051, 1291, 1531]
881 210 [1091, 1301, 1511, 1721, 1931]
6863 420 [7283, 7703, 8123, 8543, 8963]
6323 240 [6563, 6803, 7043, 7283, 7523]
8101 210 [8311, 8521, 8731, 8941, 9151]
6869 240 [7109, 7349, 7589, 7829, 8069]
8101 210 [8311, 8521, 8731, 8941, 9151]
2659 420 [3079, 3499, 3919, 4339, 4759]
7457 30 [7487, 7517, 7547, 7577, 7607]
4129 210 [4339, 4549, 4759, 4969, 5179]
2039 420 [2459, 2879, 3299, 3719, 4139]
41 420 [461, 881, 1301, 1721, 2141]
1381 450 [1831, 2281, 2731, 3181, 3631]
829 210 [1039, 1249, 1459, 1669, 1879]
5273 210 [5483, 5693, 5903, 6113, 6323]
1201 270 [1471, 1741, 2011, 2281, 2551]
2339 120 [2459, 2579, 2699, 2819, 2939]
227 360 [587, 947, 1307, 1667, 2027]
6863 420 [7283, 7703, 8123, 8543, 8963]
8761 210 [8971, 9181, 9391, 9601, 9811]
641 60 [701, 761, 821, 881, 941]
1301 210 [1511, 1721, 1931, 2141, 2351]
2953 210 [3163, 3373, 3583, 3793, 4003]
5273 210 [5483, 5693, 5903, 6113, 6323]
3607 450 [4057, 4507, 4957, 5407, 5857]
331 240 [571, 811, 1051, 1291, 1531]
397 180 [577, 757, 937, 1117, 1297]
5651 450 [6101, 6551, 7001, 7451, 7901]
743 420 [1163, 1583, 2003, 2423, 2843]
6491 210 [6701, 6911, 7121, 7331, 7541]
5227 210 [5437, 5647, 5857, 6067, 6277]
43 420 [463, 883, 1303, 1723, 2143]
2351 90 [2441, 2531, 2621, 2711, 2801]
5779 420 [6199, 6619, 7039, 7459, 7879]
157 150 [307, 457, 607, 757, 907]
1039 210 [1249, 1459, 1669, 1879, 2089]
2437 120 [2557, 2677, 2797, 2917, 3037]
7717 150 [7867, 8017, 8167, 8317, 8467]
5651 450 [6101, 6551, 7001, 7451, 7901]
2833 420 [3253, 3673, 4093, 4513, 4933]
5021 240 [5261, 5501, 5741, 5981, 6221]
227 360 [587, 947, 1307, 1667, 2027]
4877 210 [5087, 5297, 5507, 5717, 5927]
227 360 [587, 947, 1307, 1667, 2027]
3607 450 [4057, 4507, 4957, 5407, 5857]
8467 270 [8737, 9007, 9277, 9547, 9817]
6949 420 [7369, 7789, 8209, 8629, 9049]
4129 210 [4339, 4549, 4759, 4969, 5179]
2857 330 [3187, 3517, 3847, 4177, 4507]
4129 210 [4339, 4549, 4759, 4969, 5179]
6863 420 [7283, 7703, 8123, 8543, 8963]
13 90 [103, 193, 283, 373, 463]
13 210 [223, 433, 643, 853, 1063]
641 60 [701, 761, 821, 881, 941]
4637 150 [4787, 4937, 5087, 5237, 5387]
389 210 [599, 809, 1019, 1229, 1439]
5081 90 [5171, 5261, 5351, 5441, 5531]
5273 210 [5483, 5693, 5903, 6113, 6323]
41 420 [461, 881, 1301, 1721, 2141]
613 420 [1033, 1453, 1873, 2293, 2713]
199 210 [409, 619, 829, 1039, 1249]
}}}