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

4 4
3 0
}}}

<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 N 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|
@interact
def _(N=(2,(2..200)), M=(100,(100..500))):
    for i in range(N, max(N,M)):
        if i == factor(i)[0][0]:
            print i, " is prime."
///

<html><!--notruncate--><div padding=6 id='div-interact-20'> <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-20' class='ui-slider ui-slider-3' style='margin:0px;'><span class='ui-slider-handle'></span></div>
    	</td><td><font color='black' id='slider-N-20-lbl'></font></td></tr></table><script>(function(){ var values = ["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","54","55","56","57","58","59","60","61","62","63","64","65","66","67","68","69","70","71","72","73","74","75","76","77","78","79","80","81","82","83","84","85","86","87","88","89","90","91","92","93","94","95","96","97","98","99","100","101","102","103","104","105","106","107","108","109","110","111","112","113","114","115","116","117","118","119","120","121","122","123","124","125","126","127","128","129","130","131","132","133","134","135","136","137","138","139","140","141","142","143","144","145","146","147","148","149","150","151","152","153","154","155","156","157","158","159","160","161","162","163","164","165","166","167","168","169","170","171","172","173","174","175","176","177","178","179","180","181","182","183","184","185","186","187","188","189","190","191","192","193","194","195","196","197","198","199","200"]; setTimeout(function() {
    $('#slider-N-20').slider({
    	stepping: 1, min: 0, max: 198, startValue: 0,
    	change: function (e,ui) { var position = ui.value; if(values!=null) $('#slider-N-20-lbl').text(values[position]); interact(20, "sage.server.notebook.interact.update(20, \"N\", 10, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals());sage.server.notebook.interact.recompute(20)"); },
    	slide: function(e,ui) { if(values!=null) $('#slider-N-20-lbl').text(values[ui.value]); }
    });
    if(values != null) $('#slider-N-20-lbl').text(values[$('#slider-N-20').slider('value')]);
    }, 1); })();</script></td></tr>
<tr><td align=right><font color="black">M&nbsp;</font></td><td><table><tr><td>
    	<div id='slider-M-20' class='ui-slider ui-slider-3' style='margin:0px;'><span class='ui-slider-handle'></span></div>
    	</td><td><font color='black' id='slider-M-20-lbl'></font></td></tr></table><script>(function(){ var values = ["100","101","102","103","104","105","106","107","108","109","110","111","112","113","114","115","116","117","118","119","120","121","122","123","124","125","126","127","128","129","130","131","132","133","134","135","136","137","138","139","140","141","142","143","144","145","146","147","148","149","150","151","152","153","154","155","156","157","158","159","160","161","162","163","164","165","166","167","168","169","170","171","172","173","174","175","176","177","178","179","180","181","182","183","184","185","186","187","188","189","190","191","192","193","194","195","196","197","198","199","200","201","202","203","204","205","206","207","208","209","210","211","212","213","214","215","216","217","218","219","220","221","222","223","224","225","226","227","228","229","230","231","232","233","234","235","236","237","238","239","240","241","242","243","244","245","246","247","248","249","250","251","252","253","254","255","256","257","258","259","260","261","262","263","264","265","266","267","268","269","270","271","272","273","274","275","276","277","278","279","280","281","282","283","284","285","286","287","288","289","290","291","292","293","294","295","296","297","298","299","300","301","302","303","304","305","306","307","308","309","310","311","312","313","314","315","316","317","318","319","320","321","322","323","324","325","326","327","328","329","330","331","332","333","334","335","336","337","338","339","340","341","342","343","344","345","346","347","348","349","350","351","352","353","354","355","356","357","358","359","360","361","362","363","364","365","366","367","368","369","370","371","372","373","374","375","376","377","378","379","380","381","382","383","384","385","386","387","388","389","390","391","392","393","394","395","396","397","398","399","400","401","402","403","404","405","406","407","408","409","410","411","412","413","414","415","416","417","418","419","420","421","422","423","424","425","426","427","428","429","430","431","432","433","434","435","436","437","438","439","440","441","442","443","444","445","446","447","448","449","450","451","452","453","454","455","456","457","458","459","460","461","462","463","464","465","466","467","468","469","470","471","472","473","474","475","476","477","478","479","480","481","482","483","484","485","486","487","488","489","490","491","492","493","494","495","496","497","498","499","500"]; setTimeout(function() {
    $('#slider-M-20').slider({
    	stepping: 1, min: 0, max: 400, startValue: 0,
    	change: function (e,ui) { var position = ui.value; if(values!=null) $('#slider-M-20-lbl').text(values[position]); interact(20, "sage.server.notebook.interact.update(20, \"M\", 11, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals());sage.server.notebook.interact.recompute(20)"); },
    	slide: function(e,ui) { if(values!=null) $('#slider-M-20-lbl').text(values[ui.value]); }
    });
    if(values != null) $('#slider-M-20-lbl').text(values[$('#slider-M-20').slider('value')]);
    }, 1); })();</script></td></tr>
</table><div id='cell-interact-20'><?__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>Sieve of Eratosthenes</h3>
<p>Factoring each number takes a long time and an ancient (and faster) way to find the primes below a number N is called the Sieve of Eratosthenes, attributed to an ancient Greek living around 200 BC.&nbsp; The idea is to start with a list of all the numbers below N, considered to be possibly prime and one by one, remove the multiples of each known prime.&nbsp; Whatever is left is prime since it is not a multiple of some smaller prime.&nbsp;</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 p in possible_primes:     
      if (p%a_known_prime == 0) and (p != a_known_prime):
          checked_primes.remove(p)
    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
          
@interact
def _(N=(54,(5..200))):
    possible_primes = list(range(2, N))
    list_o_primes = seive(2,possible_primes)
    print "The primes less than ",N," are:", list_o_primes
///

<html><!--notruncate--><div padding=6 id='div-interact-34'> <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-34' class='ui-slider ui-slider-3' style='margin:0px;'><span class='ui-slider-handle'></span></div>
    	</td><td><font color='black' id='slider-N-34-lbl'></font></td></tr></table><script>(function(){ var values = ["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","54","55","56","57","58","59","60","61","62","63","64","65","66","67","68","69","70","71","72","73","74","75","76","77","78","79","80","81","82","83","84","85","86","87","88","89","90","91","92","93","94","95","96","97","98","99","100","101","102","103","104","105","106","107","108","109","110","111","112","113","114","115","116","117","118","119","120","121","122","123","124","125","126","127","128","129","130","131","132","133","134","135","136","137","138","139","140","141","142","143","144","145","146","147","148","149","150","151","152","153","154","155","156","157","158","159","160","161","162","163","164","165","166","167","168","169","170","171","172","173","174","175","176","177","178","179","180","181","182","183","184","185","186","187","188","189","190","191","192","193","194","195","196","197","198","199","200"]; setTimeout(function() {
    $('#slider-N-34').slider({
    	stepping: 1, min: 0, max: 195, startValue: 49,
    	change: function (e,ui) { var position = ui.value; if(values!=null) $('#slider-N-34-lbl').text(values[position]); interact(34, "sage.server.notebook.interact.update(34, \"N\", 18, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals());sage.server.notebook.interact.recompute(34)"); },
    	slide: function(e,ui) { if(values!=null) $('#slider-N-34-lbl').text(values[ui.value]); }
    });
    if(values != null) $('#slider-N-34-lbl').text(values[$('#slider-N-34').slider('value')]);
    }, 1); })();</script></td></tr>
</table><div id='cell-interact-34'><?__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>
}}}

<p>This example illustrates a recursive function. (Each call of seive removes the multiples of one prime; if there are more prime multiples to remove, <em>seive calls itself</em> with the next prime and the list of remaining possible primes.) It also illustrates the binary logical operators and and or.&nbsp; Inside the for loop, we check if both p is divisible by a_known_prime and also not equal to it (!= means not equal).&nbsp; In addition, there are some features of Sage lists to highlight.&nbsp; Notice that possible_primes is a list built from the range function.&nbsp; Since it is a list, we can dynamically alter its elements. For example, we can .remove() an element from the list, find the .index() of a list, and find its length using the len() function.&nbsp; The function append(), not used in this example, adds a new element to the end of a list.</p>
<p>&nbsp;</p>
<p>Sage has a built in function that checks the primality of an integer called is_prime().&nbsp; The previous example can be rewritten as:</p>

{{{id=31|
for i in range(2,100):
    if is_prime(i):
        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><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(105,121)
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
105 = 105 - 121 * 0

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))):
    var('x')
    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\", 3, 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 Progressions 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
///

1289 420 [1709, 2129, 2549, 2969, 3389]
1321 210 [1531, 1741, 1951, 2161, 2371]
7841 330 [8171, 8501, 8831, 9161, 9491]
1973 90 [2063, 2153, 2243, 2333, 2423]
179 210 [389, 599, 809, 1019, 1229]
1039 210 [1249, 1459, 1669, 1879, 2089]
5179 300 [5479, 5779, 6079, 6379, 6679]
2857 330 [3187, 3517, 3847, 4177, 4507]
1301 210 [1511, 1721, 1931, 2141, 2351]
5273 210 [5483, 5693, 5903, 6113, 6323]
1289 420 [1709, 2129, 2549, 2969, 3389]
4021 420 [4441, 4861, 5281, 5701, 6121]
389 210 [599, 809, 1019, 1229, 1439]
7321 420 [7741, 8161, 8581, 9001, 9421]
2239 420 [2659, 3079, 3499, 3919, 4339]
3079 420 [3499, 3919, 4339, 4759, 5179]
3499 210 [3709, 3919, 4129, 4339, 4549]
359 30 [389, 419, 449, 479, 509]
281 120 [401, 521, 641, 761, 881]
8101 210 [8311, 8521, 8731, 8941, 9151]
1669 360 [2029, 2389, 2749, 3109, 3469]
1091 210 [1301, 1511, 1721, 1931, 2141]
3919 210 [4129, 4339, 4549, 4759, 4969]
619 210 [829, 1039, 1249, 1459, 1669]
8101 210 [8311, 8521, 8731, 8941, 9151]
1039 210 [1249, 1459, 1669, 1879, 2089]
4603 210 [4813, 5023, 5233, 5443, 5653]
1453 210 [1663, 1873, 2083, 2293, 2503]
11 60 [71, 131, 191, 251, 311]
11 480 [491, 971, 1451, 1931, 2411]
577 360 [937, 1297, 1657, 2017, 2377]
6079 150 [6229, 6379, 6529, 6679, 6829]
491 330 [821, 1151, 1481, 1811, 2141]
4637 150 [4787, 4937, 5087, 5237, 5387]
83 300 [383, 683, 983, 1283, 1583]
1663 210 [1873, 2083, 2293, 2503, 2713]
199 210 [409, 619, 829, 1039, 1249]
2339 120 [2459, 2579, 2699, 2819, 2939]
587 300 [887, 1187, 1487, 1787, 2087]
619 210 [829, 1039, 1249, 1459, 1669]
4877 210 [5087, 5297, 5507, 5717, 5927]
743 420 [1163, 1583, 2003, 2423, 2843]
2129 330 [2459, 2789, 3119, 3449, 3779]
7 30 [37, 67, 97, 127, 157]
7 150 [157, 307, 457, 607, 757]
5273 210 [5483, 5693, 5903, 6113, 6323]
2221 30 [2251, 2281, 2311, 2341, 2371]
257 210 [467, 677, 887, 1097, 1307]
1301 210 [1511, 1721, 1931, 2141, 2351]
2207 240 [2447, 2687, 2927, 3167, 3407]
1451 450 [1901, 2351, 2801, 3251, 3701]
2347 420 [2767, 3187, 3607, 4027, 4447]
2953 210 [3163, 3373, 3583, 3793, 4003]
5153 330 [5483, 5813, 6143, 6473, 6803]
}}}

{{{id=39|

///
}}}