number_theory.primes_diophantine_equations_1.0
system:sage

<h2 style="text-align: center;"><span class="cmr-17">Introduction to Number Theory in SAGE: A Primer</span></h2>
<h3 style="text-align: center;">Erik Jacobson<br /></h3>
<h4 style="text-align: center;">University of Georgia</h4>
<p class="noindent" style="text-align: center;"><br /> <span class="cmbx-12">Version 1.0</span><br /> <span class="cmbx-12">June 21, 2009</span></p>
<p>"To a large degree, number theory is an experimental science."- John William Scott Cassels</p>
<h2><span style="font-family: book antiqua,palatino;">Introduction<br /></span></h2>
<p><span style="font-family: book antiqua,palatino;">Number theory is perhaps the mathematical field that most makes use of inductive reasoning; Gauss, for example, was said to have made his famous conjecture concerning the limiting behavior of $\pi(n)$ after studying values of $\frac{\pi(n) - \pi(n-1000)}{1000}$ for large n. If you are interested in number theory, you will be interested in the data computers can provide.<br /></span></p>
<p><span style="font-family: book antiqua,palatino;"> This goal of this interactive worksheet is not to teach either Sage or number theory; </span><span style="font-family: book antiqua,palatino;">for these purposes, </span><span style="font-family: book antiqua,palatino;">there already exist Sage tutorials (http://www.sagemath.org/doc/tutorial/) and number theory resources (http://math.uga.edu/~pete/numbertheory2009.html).&nbsp; Rather the purpose is to provide substantive examples that motivate the use of Sage for number theoretic explorations and illustrate some of what is possible when Sage is used for number theory.&nbsp; Our goal is not deep understanding but a piqued curiosity.</span><span style="font-family: book antiqua,palatino;">If we succeed, the details will soon be essential: sagemath.org and python.org are excellent resources.</span></p>
<p><span style="font-family: book antiqua,palatino;">The primer is divided into two sections, each with two subsections.&nbsp; In the section entitled <em>Primes</em>, we introduce modular division, factorization, primality testing and demonstrate implementations of the sieve of Eratosthenes and the Euclidean algorithm. In addition, illustrations of the </span><span style="font-family: book antiqua,palatino;">built-in $\pi(n)$ function introduce</span><span style="font-family: book antiqua,palatino;"> 2D plotting in Sage using plot() and the matplotlib library. The section concludes with code that randomly searches for (short) arithmetic progressions of primes.</span><span style="font-family: arial,helvetica,sans-serif;"><br /></span></p>
<p><span style="font-family: book antiqua,palatino;">In the section entitled <em>Diophantine Equations</em>, we first explore linear Diophantine equations, solving them using the extended gcd algorithm. Next we consider diagonal quadratic forms, performing a computational analysis of sums of two and three squares and implementing a formula using Kronecker symbols to calculate sums of an arbitrary number of squares.&nbsp; We conclude with </span><span style="font-family: book antiqua,palatino;">Pythagorean triples, demonstrating </span><span style="font-family: book antiqua,palatino;">an interactive triple-generating function and a plot of all triples with two elements less than 500.<br /></span></p>
<h3><span style="font-family: book antiqua,palatino;">A Partial List of Topics<br /></span></h3>
<p><span style="font-family: book antiqua,palatino;"><span style="text-decoration: underline;">Programing and Python</span></span></p>
<p><span style="font-family: book antiqua,palatino;">function definition. def <em>name</em>(<em>arg1</em>, <em>arg2</em>=<em>initial_value</em>, ..., <em>argN</em>):</span></p>
<p><span style="font-family: book antiqua,palatino;">for-loops. for <em>var</em> in <em>sequence</em>:</span></p>
<p><span style="font-family: book antiqua,palatino;">in-line list creation. (<em>expression_x</em> for <em>x</em> in <em>list</em>)</span></p>
<p><span style="font-family: book antiqua,palatino;">interactive function definition. @<em>interact</em> before function definition, vars determine the controls</span></p>
<p><span style="font-family: book antiqua,palatino;">Sage enumeration. (a..b)</span></p>
<p><span style="font-family: book antiqua,palatino;">type coercion. QQ(.125) ==&gt; Sage rational number with value 1/8; also str(15) =&gt; "15"</span></p>
<p><span style="font-family: book antiqua,palatino;">adding graphic objects for show(). show(plot()+plot()) ==&gt; the plots are shown on the same set of axes</span></p>
<p><span style="font-family: book antiqua,palatino;"><span style="text-decoration: underline;">Sage functions:</span></span></p>
<p><span style="font-family: book antiqua,palatino;">range([a],b [, d]), append(elt), sum(list),</span></p>
<p><span style="font-family: book antiqua,palatino;">sqrt(num), is_square(int),</span></p>
<p><span style="font-family: book antiqua,palatino;">kronecker_symbol(a,b), lcm(m,n), gcd(m,n), crt(a,b,m,n), n.divides(m), divisors(n)</span></p>
<p><span style="font-family: book antiqua,palatino;">vars("<em>new_var</em>"), min(list), max(list), solve(<em>list_of_equations</em>, <em>var1</em>, <em>var2</em>)</span></p>
<p><span style="font-family: book antiqua,palatino;">QQ.numerator(), .denominator()</span></p>
<p><span style="font-family: book antiqua,palatino;">point(<em>list_ordered_pairs</em>, <em>keywords</em>), line(<em>list_ordered_pairs</em>, <em>keywords</em>), plot(<em>graphic_object, x_min, x_max, keywords</em>)</span></p>
<p><span style="font-family: book antiqua,palatino;">show(..., <em>aspect_ratio</em> = 1), rgbcolor = hue(<em>between0&amp;1</em>)</span></p>
<p>&nbsp;</p>
<p><span style="font-family: book antiqua,palatino;">And now we begin...</span></p>
<p><span style="font-family: book antiqua,palatino;">First some tips for the notebook interface: &lt;shift-enter&gt; evaluates a cell (output is displayed bellow), &lt;control-;&gt; splits a cell, &lt;control-backspace&gt; joins cells, &lt;backspace&gt; deletes an empty cell, &lt;click&gt; the blue line to insert a new cell, and &lt;shift-click&gt; the blue line to insert text.</span></p>
<p><span style="font-family: book antiqua,palatino;">Next, we briefly explore some useful Sage functions</span>. <span style="font-family: book antiqua,palatino;">The for-loop command </span><span style="font-family: book antiqua,palatino;"><span style="font-family: courier new,courier;">for -var- in -list-</span> iterates the variable through each element of the list; in this case</span><span style="font-family: book antiqua,palatino;"><span style="font-family: courier new,courier;">(2..4)</span> returns a list of integers from 2 to 4 inclusive</span><span style="font-family: book antiqua,palatino;">. The command </span><span style="font-family: book antiqua,palatino;"><span style="font-family: courier new,courier;">print</span> displays output, each comma inserts a space, "\t" a tab, and "\n" a new line (not used here).&nbsp; The command </span><span style="font-family: book antiqua,palatino;"><span style="font-family: courier new,courier;">gcd(m, n)</span> returns the greatest common divisor of m and n and </span><span style="font-family: courier new,courier;">lcm(m,n)</span> returns their least common multiple.<span style="font-family: book antiqua,palatino;"> </span></p>

{{{id=49|
for m in (2..4): 
    n = m+4
    print m, n, "\t",gcd(m, n), lcm(m,n)
///
2 6 	2 6
3 7 	1 21
4 8 	4 8
}}}

<p><span style="font-family: book antiqua,palatino;">When gcd(m,n) = d, you can find s and t such that d = s*m + t*n using <span style="font-family: courier new,courier;">xgcd(m,n)</span>.&nbsp; Each format operator <span style="font-family: courier new,courier;">%s</span> holds the place of a successive element in the list <span style="font-family: courier new,courier;">(1, 3, -2)</span> returned by <span style="font-family: courier new,courier;">xgcd(5,7)</span>. The variable m divided by n has remainder <span style="font-family: courier new,courier;">m % n</span>, and if m | n then <span style="font-family: courier new,courier;">m.divides(n)</span> returns <span style="font-family: courier new,courier;">True</span>. Here <span style="font-family: courier new,courier;">str(-expr-)</span> coerces the expression into a string for concatenation (+ only makes sense if the things being added are both numbers or, in this case, both strings). The binomial coefficient $\binom{m}{n}$ is given by binomial(m,n).</span></p>

{{{id=74|
print "%s = %s*5 + %s*7"%xgcd(5,7)
print "5*"+str(17%5)+" = 17 so 5 | 17 is "+str(5.divides(17))
print "17 choose 5 is", binomial(17,5) #this is a comment: Sage ignores all after hash
///
1 = 3*5 + -2*7
5*2 = 17 so 5 | 17 is False
17 choose 5 is 6188
}}}

<h2><span style="font-family: book antiqua,palatino;">Primes<br /></span></h2>
<h3><span style="font-family: book antiqua,palatino;">Factoring<br /></span></h3>
<p><span style="font-family: book antiqua,palatino;">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.</span></p>

{{{id=75|
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><span style="font-family: book antiqua,palatino;">In the integers, a number p is prime if and only if it is irreducible (it has only trivial factors: 1, -1, p, -p).&nbsp; The next example checks the numbers n between N and M for primality using the function factor().&nbsp; The statement factor(i)[0][0] returns the smallest (and possibly only) factor of of i. If this factor is n, then p is prime.<br /></span></p>

{{{id=76|
@interact
def _(N=(18,(1..200)), M=(38,(1..200))):
    for n in (min(N,M)..max(N,M)):
        if n == factor(n)[0][0]:
            print n, "is prime."
///
}}}

<h3><span style="font-family: book antiqua,palatino;">Sieve of Eratosthenes</span></h3>
<p><span style="font-family: book antiqua,palatino;">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;</span></p>

{{{id=86|
def sieve(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 "Check for and remove any multiples of ", a_known_prime
    print
    possible_primes = checked_primes
    if a_known_prime < possible_primes[len(possible_primes)-1]/2:        
        return sieve(possible_primes[possible_primes.index(a_known_prime)+1], possible_primes)
    else:
        return possible_primes
          
@interact
def _(N=(38,(5..200))):
    possible_primes = list(range(2, N))
    list_o_primes = sieve(2,possible_primes)
    print "The primes less than ",N," are:", list_o_primes
///
}}}

<p><span style="font-family: book antiqua,palatino;">This example illustrates a recursive function. (Each call of sieve removes the multiples of one prime; if there are more prime multiples to remove, sieve<em> calls itself</em></span> 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><span style="font-family: book antiqua,palatino;">Sage has a built in function that checks the primality of an integer called is_prime().&nbsp; The example from above can be rewritten as:</span></p>

{{{id=87|
@interact
def _(N=(18,(1..200)), M=(38,(1..200))):
    for n in (min(N,M)..max(N,M)):
        if is_prime(n):
            print n, "is prime."
///
}}}

<h3><span style="font-family: book antiqua,palatino;">The Euclidean Algorithm</span></h3>
<p><span style="font-family: book antiqua,palatino;">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.</span></p>

{{{id=88|
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><span style="font-family: book antiqua,palatino;">Sage has the built-in function gcd() that finds the greatest common factor of two or more numbers. Notice "\n" which introduces a new line in the <span style="font-family: courier new,courier;">print</span> output.<br /></span></p>

{{{id=89|
print "The gcd of 121 and 105 is", gcd(121,105),"\nThe gcd of 67890 and 12345 is",gcd(67890,12345)
///
The gcd of 121 and 105 is 1 
The gcd of 67890 and 12345 is 15
}}}

<h3><span style="font-family: book antiqua,palatino;">Counting Primes: Ploting the Pi Function</span></h3>
<p><span style="font-family: book antiqua,palatino;">That there are an infinite number of primes was known to ancient Greeks.&nbsp; The function <span style="font-family: courier new,courier;">primes_pi()</span> in Sage returns the number of primes below a given parameter.&nbsp; The example below is an interactive plot comparing <span style="font-family: courier new,courier;">primes_pi()</span> and its asymptotic bound x/log(x)-1.&nbsp; The <span style="font-family: courier new,courier;">@interact</span> command allows you to create dynamic plots.&nbsp; After <span style="font-family: courier new,courier;">@interact</span>, define a function that 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 <span style="font-family: courier new,courier;">prime_pi() </span>function superimposed on a plot of x/log(x) -1 in blue.&nbsp; Both are graphed on the domain bounded above by N.</span></p>

{{{id=90|
# due to William Stein.
@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'))
///
}}}

<h3><span style="font-family: book antiqua,palatino;">Plotting the Density of Primes</span></h3>
<p><span style="font-family: book antiqua,palatino;">Here we wish to visualize the density of primes over the interval (0, 20000),&nbsp; another perspective on </span><span style="font-family: book antiqua,palatino;"><span style="font-family: courier new,courier;">prime_pi()</span></span><span style="font-family: book antiqua,palatino;"> function plotted above.&nbsp; Aparanently, when Gauss was 15, he studied the slope of these segments before making his famous conjecture about the asymptotic behavior of <span style="font-family: courier new,courier;">prime_pi()</span>. </span><span style="font-family: book antiqua,palatino;">You can use the matplotlib (http://matplotlib.sourceforge.net) library if you need more control over plotting than the <span style="font-family: courier new,courier;">plot()</span> function provides. In Sage, the commands for matplotlib can be used only after importing them.&nbsp; This examples illustrates some functionality of matplotlib. </span></p>

{{{id=77|
x = map(lambda x: x*10^3+10^3,range(20))
y = [prime_pi(10^3)]+map(lambda j: (prime_pi(x[j])-prime_pi(x[j-1])), range(1,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^3 from 0 to 2*10^4')
ax.set_xlabel('Upper bound of intervals of width 4*10^3')
ax.set_ylabel('Number of Primes')
ax.plot(x,y, 'go-', linewidth=1.0 )
canvas.print_figure('ex1_linear.png')
///
}}}

<p><span style="font-family: book antiqua,palatino;">The code above stores an 20 term arithmetic progression of integers with a constant difference of 10^3 in the list x, the x-coordinates.&nbsp; Then the number of primes in each interval of width 10^3 is stored in the list y, the y coordinates. The axis lables are set with the by the <span style="font-family: courier new,courier;">ax.set_xlabel()</span> and <span style="font-family: courier new,courier;">ax.set_ylabel()</span> commands. The title is set with <span style="font-family: courier new,courier;">ax.set_title()</span>.&nbsp; The number of grid lines is set dependent on the lists x and y using <span style="font-family: courier new,courier;">MaxNLocator(-num-)</span> passed to <span style="font-family: courier new,courier;">ax.xaxis.set_major_locator()</span> and <span style="font-family: courier new,courier;">ax.yaxis.set_major_locator()</span>, respectively.</span></p>
<h3><span style="font-family: book antiqua,palatino;">Arithmetic Progressions of Primes<br /></span></h3>
<p><span style="font-family: book antiqua,palatino;">Arithmetic progressions of primes of any length are conjectured to exist, but finding them is difficult.&nbsp; As of March 2009, the longest known progression of primes has 25 terms.&nbsp; Shorter sequences are (luckily) easy to find. The idea in the code below is to pick a random prime and search for arithmetic progressions starting at this prime with differences between 1 and say, 200.&nbsp; </span></p>
<p><span style="font-family: book antiqua,palatino;">The code stores the length of the desired sequence in <span style="font-family: courier new,courier;">seq_length</span>.&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 200 times, using the for loop indexed by i.&nbsp; The variable <span style="font-family: courier new,courier;">j</span> stores a random index into the list of primes <span style="font-family: courier new,courier;">P</span>.&nbsp; For each possible difference, scale the list from 1 to </span><span style="font-family: book antiqua,palatino;"><span style="font-family: courier new,courier;">seq_length</span></span><span style="font-family: book antiqua,palatino;"> by<span style="font-family: courier new,courier;"> k</span> and shift by the randomly chosen prime <span style="font-family: courier new,courier;">P[j]</span>, then check to see if all the numbers in this list are primes.&nbsp; Finally, print the initial prime, the difference, and the sequence of primes. If we find no primes, run it again!<br /></span></p>

{{{id=78|
print "prime \tdiff. \t[progression]","\n-----\t-----\t-------------"

seq_length = 5
P=prime_range(10^4)
for i in range(1,200):
    j = int(random()*len(P))
    for k in range(1,200):
        arithmetic_sequence = map(lambda x: x*k+P[j], range(1, seq_length+1))
        if set(arithmetic_sequence).issubset(set(P)):
            print P[j],"\t", k,"\t", arithmetic_sequence
///
prime 	diff. 	[progression] 
-----	-----	-------------
1013 	180 	[1193, 1373, 1553, 1733, 1913]
157 	150 	[307, 457, 607, 757, 907]
2351 	90 	[2441, 2531, 2621, 2711, 2801]
4637 	150 	[4787, 4937, 5087, 5237, 5387]
}}}

<h2>Diophantine Equations</h2>
<p><span style="font-family: book antiqua,palatino;">
<h3 class="line874">Linear Diophantine Equations</h3>
<p class="line874">We solve the equation 2345x + 4567y = 1 using the extended gcd function, xgcd(). First, we check if 2345 and 4567 are relatively prime.</p>
</span></p>

{{{id=51|
a = 2345
b = 4567
d = gcd(2345, 4567); print d
///
1
}}}

<p><span style="font-family: book antiqua,palatino;">Next, we find s and t such that 2345*s + 4567*t = 1 using <span style="font-family: courier new,courier;">xgcd()</span>.</span></p>

{{{id=55|
s,t = xgcd(a,b)[1:] # [1:] returns the list starting at the 2nd value, index 1
print "2345*"+str(s)+" + 4567*"+str(t)+" = 1"
///
2345*1708 + 4567*-877 = 1
}}}

<p><span style="font-family: book antiqua,palatino;"><span style="text-decoration: underline;">Proposition</span> (stated w/o proof): The solutions of a linear Diophantine equation ax + by = c are x = x0 + (b/d)t, y = y0 + (a/d)t, t in Z, where d = gcd(a,b) if and only if d | c.</span></p>
<p><span style="font-family: book antiqua,palatino;">Following the proposition, we generate and plot solution pairs for t from -9 to 9.</span></p>

{{{id=57|
x0 = s # you must run the preceding cell to initialize s, t
y0 = t
solns = []
for t in (-3..3):
    solns.append((x0+b*t, y0-a*t))

print solns
show(point(solns, rgbcolor = 'purple'))
///
[(-11993, 6158), (-7426, 3813), (-2859, 1468), (1708, -877), (6275, -3222), (10842, -5567), (15409, -7912)]
}}}

<p><span style="font-family: book antiqua,palatino;">Suppose you wish to find k solutions of n linear Diophantine equations. Even for n&gt;5, this becomes tedious using the method outlined above.&nbsp; To conclude this section, we generalize our prior work in a new user-defined Sage function solve_Diophantine(), use the function to generate solutions to a list of randomly-generated linear Diophantine equations, plotting the solutions on the same axes.</span></p>

{{{id=52|
#generate a random list of coefficients of linear Diophantine equations
coefs_list = []
for i in (10..13):
    coefs_list.append((random_int_upto(100)+1, random_int_upto(100)+1, i))

#if possible, return k (default 5) solutions for ax + by = c, where cfs = (a,b,c)
def solve_Diophantine(cfs, k=5): 
    solns = []
    a,b,c = cfs
    d,x0,y0 = xgcd(a,b)
    if d.divides(c):
        for t in (-ceil(k/2)..ceil(k/2)):
            solns.append((x0+(b/d)*t, y0-(a/d)*t))
    else:
        solns.append(None)
    return solns


#for coefficient set, generate solutions, store in soln_list
soln_list = []
for coefs in coefs_list:
    soln = solve_Diophantine(coefs)
    if soln != [None]:
        soln_list.append((soln,coefs))

#plot origin and the plot label
p = point((0,0))
t = text("      Linear Diophantine Equations",(0,0), \
                                                rgbcolor = 'purple', \
                                                vertical_alignment='bottom', \
                                                horizontal_alignment='left')

#plot each set of solution points and an equation label for each
n = len(soln_list)
for i in range(n):
    p += point(soln_list[i][0], hue = (i/Integer(n)))
    t += text( " %s*x + %s*y = %s"%soln_list[i][1], soln_list[i][0][0], \
                                                hue =(i/Integer(n)), \
                                                vertical_alignment='bottom', \
                                                horizontal_alignment='left')
show(p+t)
///
}}}

<h3><span style="font-family: book antiqua,palatino;">Congruences Mod <em>m</em></span><br /></h3>
<p><span style="font-family: book antiqua,palatino;">From linear Diophantine equations, we digress to the equivalent problem of solving congruences mod m, and use Sage to approach the problem from this angle.&nbsp; Consider again $2345x + 4567y = 1$ and note that it is equivalent to the statement $2345x \equiv 1 \pmod{4567}$. From before, we have a list of 7 solutions to this equation, and we verify that each one satisfies the modular congruence.<br /></span></p>

{{{id=58|
solns = [(-11993, 7044), (-7426, 4699), (-2859, 2354), (1708, 9), (6275, -2336),
(10842, -4681), (15409, -7026)]
for soln in solns: print Mod(2345*soln[0], 4567)==1
///
True
True
True
True
True
True
True
}}}

<p class="line874"><span style="font-family: book antiqua,palatino;">Sucess!</span></p>
<p class="line874"><span style="font-family: book antiqua,palatino;">Now we briefly discuss solving systems of congruence using the Chinese remainder theorem, implemented in Sage as crt(a,b,m,b) which returns x such that $x \equiv ~ a \pmod{m}$ and&nbsp; $x \equiv ~ b \pmod{n}$.&nbsp; Note that this function is not defined unless a and b are coprime.</span></p>
<p class="line874"><span style="font-family: book antiqua,palatino;"><br />Suppose that in addition to $2345x \equiv 1 \pmod{4567}$ we have the condition $789x \equiv 1 \pmod{1213}$. What values of x satisfy both conditions?&nbsp; Since 4567 and 1213 are prime, all elements 2345 and 789 are units and we can solve both congruences for x by computing the inverse of 2345 and 789 (mod 4567 and 1213, respectively), using inverse_mod().&nbsp; Then we compute a solution using crt() and check the solution using Mod().<br /></span></p>

{{{id=60|
a,b = inverse_mod(2345,4567),inverse_mod(789,1213)
print a,"*2345 (mod 4567) =",Mod(a*2345,4567)," and ", b,"*789 (mod 1213) =",Mod(b*789,1213)
x =crt(a,b,4567,1213)
print x,"=",a,"(mod 4567) and",x,"=",b,"(mod 1213)?"
print Mod(x,4567)==a, " and ", Mod(x,1213)==b, "!"
///
1708 *2345 (mod 4567) = 1  and  884 *789 (mod 1213) = 1
-1313357884 = 1708 (mod 4567) and -1313357884 = 884 (mod 1213)?
True  and  True !
}}}

<h3><span style="font-family: book antiqua,palatino;">A Few Diagonal Quadratic Forms: Sums of <em>k</em>-Squares</span></h3>
<p><span style="font-family: book antiqua,palatino;">(adapted from "Experimental Number Theory" by Fernando Rodriquez Villegas, section 1.4)</span></p>
<p><span style="font-family: book antiqua,palatino;">The simplest diagonal quadratic forms correspond to the identity matrices of various dimension.&nbsp; Here we examine all sums of one(?), two, three, and four squares:</span></p>
<p><span style="font-family: book antiqua,palatino;">$\begin{align}n &amp;= x_1^2 \\ n &amp;= x_1^2 + x_2^2 \\ n &amp;= x_1^2 + x_2^2 + x_3^2 \\ n &amp;= x_1^2 + x_2^2 + x_3^2 + x_4^2 \\ \end{align}$</span></p>
<p><span style="font-family: book antiqua,palatino;">using Sage to both enumerate solutions for small n and find the number of solutions for arbitrary n.</span></p>
<p><span style="font-family: book antiqua,palatino;">The brute force method of finding solutions is straightforward.&nbsp; Certainly $x_i  \leq \sqrt{n}$, and so it is a simple matter to run through possible values for $x_i$ and check for solutions. Let $k$ be the number of squares in each sum.</span></p>
<p><span style="font-family: book antiqua,palatino;">When k = 1, it is certainly sufficient (although a trifle ridiculous) to check if the square of any integer less than or equal to $\sqrt{n}$ is $n$.&nbsp; After demonstrating this idea, hereafter we use the Sage command is_square() which has the same functionality.</span></p>

{{{id=8|
def sum1square(n): #the purpose here is pedagogical, hereafter we use is_square() instead
    for x in range(1, sqrt(n)+1):
        if n == x^2:
            return x
    return None

print "n:"," sum1square(n)", "  is_square(n,True)"
print "-- ", "-------------", "  ----------------"
for n in range(1,10): print str(n),"  ",sum1square(n), "            ",is_square(n,True)
///
n:  sum1square(n)   is_square(n,True)
--  -------------   ----------------
1    1              (True, 1)
2    None              (False, None)
3    None              (False, None)
4    2              (True, 2)
5    None              (False, None)
6    None              (False, None)
7    None              (False, None)
8    None              (False, None)
9    3              (True, 3)
}}}

<p><span style="font-family: book antiqua,palatino;">When k = 2, we must iterate two variables, $ 0 \leq x_1 \leq \sqrt{n}$ and $0 \leq x_2 \leq \sqrt{n}$.&nbsp; Note however, that the method outlined above neglects to exploit the solution symmetries obtained by permuting and negating coordinates.&nbsp; Further restricting the search to solutions where each $x_i &lt; x_j$ if and only if $i&lt;j$ reduces the necessary number of computations by a factor of $2^k k!$.&nbsp; The variables x1 and x2 are iterated over the values provided by Python command range(a, b) which returns a list of integers x such that $a \leq x &lt; b$.&nbsp; The code below implements this algorithm in the case $k = 2$, printing all solutions for the given n.</span></p>

{{{id=11|
def sum2squares_v1(n):
    solns = []
    
    for x1 in range(0,sqrt(n)+1):
        for x2 in range(x1,sqrt(n)+1):
            if n == x1^2 + x2^2:
                 solns.append((x1, x2))
    return solns

for n in range(1,10): print "n: "+str(n), sum2squares_v1(n)
///
n: 1 [(0, 1)]
n: 2 [(1, 1)]
n: 3 []
n: 4 [(0, 2)]
n: 5 [(1, 2)]
n: 6 []
n: 7 []
n: 8 [(2, 2)]
n: 9 [(0, 3)]
}}}

<p><span style="font-family: book antiqua,palatino;">But we can do better.&nbsp; Once x1 is fixed, there is only one value of x2 that it makes sense to check, namely $\sqrt{n-x1^2}$.&nbsp; This means we can get away with just one for-loop.&nbsp; In addition, if $n \equiv 3 \pmod{4}$ then n is not a sum of two squares (why?), and we screen for this case before beginning the computation.&nbsp;</span></p>

{{{id=13|
def sum2squares_v2(n):
    solns = []    
    if n >=0 and n%4!=3:
        for x1 in range(0,sqrt(floor(n/2))+1):
            soln = is_square(n-x1^2,True)
            if soln[0]:
                 solns.append((x1, soln[1]))
    return solns

for n in range(1,10): print "n: "+str(n), sum2squares_v2(n)
///
n: 1 [(0, 1)]
n: 2 [(1, 1)]
n: 3 []
n: 4 [(0, 2)]
n: 5 [(1, 2)]
n: 6 []
n: 7 []
n: 8 [(2, 2)]
n: 9 [(0, 3)]
}}}

<p><span style="font-family: book antiqua,palatino;">The case for k = 3 is similar.&nbsp; We can screen out all $n \equiv 7 \pmod{8}$ (why?) and add another for loop.</span></p>

{{{id=14|
def sum3squares(n):
    solns = []    
    if n >=0 and n%8!=7:
        for x1 in range(0,floor(sqrt(n/3))+1):
            n1 = n - x1^2
            for x2 in range(x1, floor(sqrt(n1/2))+1):
                soln = is_square(n1-x2^2,True)
                if soln[0]:
                     solns.append((x1, x2, soln[1]))
        return solns

for n in range(10): print "n: "+str(n), sum3squares(n)
///
n: 0 [(0, 0, 0)]
n: 1 [(0, 0, 1)]
n: 2 [(0, 1, 1)]
n: 3 [(1, 1, 1)]
n: 4 [(0, 0, 2)]
n: 5 [(0, 1, 2)]
n: 6 [(1, 1, 2)]
n: 7 None
n: 8 [(0, 2, 2)]
n: 9 [(0, 0, 3), (1, 2, 2)]
}}}

<p><span style="font-family: book antiqua,palatino;">We leave the implementation of sum4squares() to the reader.&nbsp; By Lagrange, every integer is a sum of 4 squares, so we need not go further.&nbsp;</span></p>
<p><span style="font-family: book antiqua,palatino;">Turning now to related problem, let $r_k(n)$ be the number of solutions to the equation $ n =\sum_i^k x_i^2$.&nbsp; By mixing Boolean and standard arithmetic, we can modify the existing sum2squares_v2() function to return $r_k(n)$, and count the solutions we ignored by exploiting symmetry.&nbsp;</span></p>
<p><span style="font-family: book antiqua,palatino;">We can permute the sign of $x_i$ unless it is zero and we can permute the order of $x_1$ and $x_2$ unless they are equal. Conveniently in Sage, you can treat False as the number 0 and True as the number 1. Thus we have that each sum2squares_v2() solution (x1, x2) represents ((x1&gt;0)+1)*(x2&gt;0)+1)*((x1!=x2)+1) unique solutions (see r2_v1()).</span></p>
<p><span style="font-family: book antiqua,palatino;">We can also implement the explicit formula</span></p>
<p><span style="font-family: book antiqua,palatino;">$r_2(n) = 4\sum_{d|n}{\chi (d)}$</span></p>
<p><span style="font-family: book antiqua,palatino;">by taking four times the sum of the Kronecker symbol (-4, d) over the divisors of n (see r2_v2()).</span></p>

{{{id=16|
def r2_v1(n):
    count = 0    
    if n >=0 and n%4!=3:
        for x1 in range(0,sqrt(floor(n/2))+1):
            soln = is_square(n-x1^2,True)
            if soln[0]:
                 x2 = soln[1]
                 count = count + ((x1>0)+1)*((x2>0)+1)*((x1!=x2)+1)
    return count

def r2_v2(n):
    return 4*sum(kronecker_symbol(-4,d) for d in divisors(n))

print "n","v1 ", "v2"
for n in range(1,10): 
    print str(n), r2_v1(n)," ", r2_v2(n)
for N in range(int(10^16), int(10^16+10)): print j, r2_v2(N)
///
n v1  v2
1 4   4
2 4   4
3 0   0
4 4   4
5 8   8
6 0   0
7 0   0
8 4   4
9 4   4
751 68
751 128
751 0
751 0
751 32
751 0
751 0
751 0
751 8
751 32
}}}

<h3 class="line874"><span style="font-family: book antiqua,palatino;">Pythagorean Triples</span></h3>
<p class="line874"><span style="font-family: book antiqua,palatino;">A problem that is similar to the problem of finding all pairs of two squares that sum to a given number is the enumeration of Pythagorean triples. We want to find all n, x1, x2 such that</span></p>
<p class="line874"><span style="font-family: book antiqua,palatino;">$ n^2 = x_1^2 + x_2^2$</span></p>
<p class="line874"><span style="font-family: book antiqua,palatino;">The brute force strategy? Check pairs of squares to see if the sum is square. We begin by modifying sum2squares_v2 to generate Pythagorean triples with <em>x1</em> and <em>x2</em> less than the variable <em>max</em>.</span></p>

{{{id=25|
def pythag_triples(max):
    triples = []
    
    for x1 in range(1,max+1):
        for x2 in range(1,max+1):
            if is_square(x1^2 + x2^2):
                 triples.append((sqrt(x1^2+x2^2), x1, x2))
    return triples

print pythag_triples(16)
///
[(5, 3, 4), (5, 4, 3), (13, 5, 12), (10, 6, 8), (10, 8, 6), (17, 8, 15), (15, 9, 12), (13, 12, 5), (15, 12, 9), (20, 12, 16), (17, 15, 8), (20, 16, 12)]
}}}

<p><span style="font-family: book antiqua,palatino;">But even among these first three few examples, one observes some triples are related.&nbsp; For example, (5,4,3) is related to (5,3,4)&nbsp; by permutation, and each element of (5, 4, 3) is doubled in (10, 8, 6).&nbsp; Algebraic reasoning shows that any multiple of a triple is a triple, so how to enumerate the generating triples, the list of triples that are "unique up to multiples," becomes the next significant question.&nbsp;</span></p>
<p><span style="font-family: book antiqua,palatino;">One idea is to divide each triple by n, then each generating triple corresponds to the rational point (x1/n, x2/n) on the unit circle. It is a fact that a point on the circle has rational coordinates if and only if it is the second intersection point of a line passing through (-1,0) with rational slope.&nbsp; Becuase of the symmetry of the circle, it suffices to examine such points in the first quadrant and these points are the intersections of the line $y=\frac{a}{b}(x+1)x, ~ a,b \in \mathbb{Z}^+$ and the unit circle.</span></p>
<p><span style="font-family: book antiqua,palatino;">The interactive function below displays the line, the corresponding triple for each intersection, and a triangle similar to the one represented by the Pythagorean triple.</span></p>

{{{id=33|
m =100
@interact
def pythag_intersections(a = (1..m),b=(m/2..m)):
    var("x,y")
    circ = sqrt(1-x^2)
    slope = min((a/b),1)
    lin = slope*(x+1)
    
    point = solve([1 == x^2+y^2, y==lin], x, y,solution_dict = True)[1]
    numer_a, numer_b, denom = QQ(point[x]).numerator(),QQ(point[y]).numerator(),QQ(point[x]).denominator()
    
    print "The slope of the line is "+str(slope)+".\n"
    print point, "\t==>\t ("+str(denom)+", "+str(numer_a)+", "+str(numer_b)+")" 
    
    c = plot(circ, x, (-1, 1),rgbcolor='green')
    l = plot(lin, x, (-1,1), rgbcolor='purple')
    t = line([[0,0],[point[x],point[y]],[point[x],0],[0,0]], rgbcolor='red')
    show(c+l+t,aspect_ratio=1)
///
}}}

<p><span style="font-family: book antiqua,palatino;">Next we plot the Pythagorean triples with a and b less than 1000 in the first quadrant using point2d().&nbsp; We generate list of points (a,b) using a modified version of pythag_triples(). The hue records the sum of coordinates of each point mod 4.&nbsp; First, notice the symmetry.&nbsp; Question: why do the coordinate-sum moduli of points on the two most distinct lines in the center of the graph cycle in the opposite order than the moduli of points on the slightly less distinct lines above and below?</span></p>

{{{id=44|
def pythag_triples_v2(max):
    points = []   
    for x1 in range(1,max+1):
        for x2 in range(1,max+1):
            if is_square(x1^2 + x2^2):
                 points.append((x1, x2))
    return points
    
pts = point((0,0))
for pair in pythag_triples_v2(500):
    pts+=point(pair,rgbcolor = hue((sum(pair)%4)/4))
show(pts, aspect_ratio = 1)
///
}}}

<p><span style="font-family: book antiqua,palatino;">We conclude with a discussion of the following<span style="text-decoration: underline;"> </span></span></p>
<p><span style="font-family: book antiqua,palatino;"><span style="text-decoration: underline;">Proposition</span> (w/o proof): Let p be a prime number congruent to 1 mod 4. There exists a right triangle with integer sides such that the length of the hypotenuse is p.</span></p>
<p><span style="font-family: book antiqua,palatino;">We exploit the fact that these primes have factorizations in Z[i ] = $\{a + bi | a, b \in \mathbb{Z}\}$.&nbsp; First, we illustrate this by using Sage to factor 13 in Z[i].&nbsp; Our strategy is to construct the Gaussian integers, and use the built-in <span style="font-family: courier new,courier;">factor()</span> function, stipulating further that the factors be given in generator-reduced form by using <span style="font-family: courier new,courier;">gen_reduced()</span>.</span></p>

{{{id=65|
var('x')
f = x^2 +1
K.<I> = NumberField(f)

print "K is a",K
factors_13 = K.factor(13)[0][0].gens_reduced()[0],K.factor(13)[1][0].gens_reduced()[0]
print "The factors of 13 in K are: \n",factors_13[0] , " and ", factors_13[1]
///
K is a Number Field in I with defining polynomial x^2 + 1
The factors of 13 in K are: 
-3*I - 2  and  3*I - 2
}}}

<p><span style="font-family: book antiqua,palatino;">Next, we square a factor of 13, taking the square root of part, real and imaginary, and so generate a Pythangorean triple.</span></p>

{{{id=63|
q = factors_13[0]^2
a,b = abs(q[0]),abs(q[1])
print str(a)+"^2 +", str(b)+"^2 = 13^2"
///
5^2 + 12^2 = 13^2
}}}

<p><span style="font-family: book antiqua,palatino;"> We generalize, calculating triples for the first 7 primes congruent to 1 mod 4 and plotting the triangles that correspond to these triples. The command <span style="font-family: courier new,courier;">filter(-bool func-, -list-)</span> returns a list with the elements of list for which the function is <span style="font-family: courier new,courier;">True</span>.&nbsp; The command <span style="font-family: courier new,courier;">lambda -var-: expr(-var-) </span>creates an inline function that takes one variable (it might be a list!) and returns an expression. We initialize an empty Graphics object so that the <span style="font-family: courier new,courier;">+=</span> works on the first pass (+= means add, then set equal to the sum).<br /></span></p>

{{{id=43|
# adapted from Alex Clemesha
max = 100
primes_list = filter(lambda x: x%4==1, primes(max))[:7]

g = Graphics() #empty Graphics object
for p in primes_list:
    qq = K.factor(p)[0][0].gens_reduced()[0]^2
    a,b = abs(qq[0]),abs(qq[1])
    lv = [[0, 0], [a, 0], [a, b], [0, 0]]
    l = line(lv, hue=(p/max))
    s = "$(%s,\ %s,\ %s)$"%(a, b, sqrt(a^2 + b^2))
    t = text(s, (a, b+1), fontsize=12,hue=(p/max))
    g+=(l+t)

show(g, aspect_ratio = 1)
///
}}}

<p style="text-align: center;"><span style="font-family: book antiqua,palatino;">~ finis ~</span></p>

{{{id=91|

///
}}}