sd15: number theory in sage
system:sage

<h1 style="text-align: center;"><span style="font-size: xx-large;"><span style="color: #0000ff;">Number Theory in Sage</span></span></h1>
<h2 style="text-align: center;"><span style="font-size: x-large;"><span style="color: #0000ff;">Craig Citro</span></span></h2>
<h2 style="text-align: center;"><span style="font-size: large;"><span style="color: #0000ff;"><span style="font-size: small;"><span style="font-size: x-large;">Sage Days 15</span><br /></span></span></span></h2>

{{{id=34|
prime_range(100)
///
}}}

{{{id=39|
len(prime_range(100))
///
}}}

{{{id=254|
nth_prime(100)
///
}}}

<h3>Sage also has fast code for computing $\pi(x)$, the prime-counting function.</h3>

<p>Sage also has fast code for computing $\pi(x)$, the prime-counting function. The algorithm in Sage (implemented by Andrew Ohana) can compute for $x$ up to about $x=2^{40}$. For contrast, Mathematica has one of (or probably <strong>the</strong>) best implementation of this among math software systems; they can compute this up to about $x=2^{45}$, and use an asymptotically better algorithm.</p>
<p>On the other hand, Fredrik Johansson is looking to implement an algorithm of Meissel, Lehmer, Lagarias, Miller, and Odlyzko here at SD15 -- so maybe we'll have better things to say in just a few days ...</p>

{{{id=33|
prime_pi(100)
///
}}}

{{{id=38|
time prime_pi(10**9)
///
}}}

{{{id=37|
plot(prime_pi, xmin=0, xmax=100)
///
}}}

{{{id=36|
plot(prime_pi, xmin=0, xmax=1000)
///
}}}

<p>The Prime Number Theorem (roughly) states that</p>
<p>$$ \pi(x) \sim \frac{x}{\log(x)}. $$</p>
<p>Here we graph them and compare $\ldots$</p>

{{{id=42|
@interact
def print_graph(max=input_box(1000)):
    f = lambda x: x/log(x)
    show(plot(f, xmin=50, xmax=max, color='red') + plot(prime_pi, xmin=50, xmax=max))
    return
///
}}}

<p>In fact, one can see that a better approximation is given by</p>
<p>$$ \pi(x) \sim \frac{x}{\log(x)-1}. $$</p>

{{{id=41|
@interact
def print_graph(max=input_box(1000)):
    f = lambda x: x/(log(x)-1)
    show(plot(f, xmin=50, xmax=max, color='red') + plot(prime_pi, xmin=50, xmax=max))
    return
///
}}}

{{{id=40|

///
}}}

<p>Of course, Sage has a number of useful facilities for factoring. Several of these aren't "turned on" by default, but recent work by Robert Bradshaw, Tom Boothby, and others means that this is imminent $\ldots$</p>

{{{id=45|
factor(2009)
///
}}}

{{{id=44|
factor(5162009)
///
}}}

{{{id=43|
n = next_prime(10^20) * next_prime(10^22)
print len(n.digits())
time factor(n)
///
}}}

{{{id=47|
factor(1234567654321234567)
///
}}}

{{{id=46|
is_prime(144169)
///
}}}

{{{id=2|

///
}}}

<h2 style="text-align: center;">Euclidean Algorithm and GCD</h2>

<p>Sage has two important functions for finding greatest common divisors: $\operatorname{gcd}$ and $\operatorname{xgcd}$.</p>
<ul>
<li>
<p>$\operatorname{gcd}(x,y)$ simply returns the greatest common divisor</p>
</li>
<li>$\operatorname{xgcd}(x,y)$ returns a pair of integers $a$, $b$ such that $$ \operatorname{gcd(x,y)} = xa + yb.$$</li>
</ul>

{{{id=1|
gcd(1819, 963)
///
}}}

{{{id=48|
xgcd(1819, 963)
///
}}}

{{{id=49|
-1 * 1819 + 2 * 963
///
}}}

{{{id=50|
@interact
def f(x=input_box('517'), y=input_box('799')):
    g, a, b = xgcd(x,y)
    print "The gcd of %s and %s is %s."%(x,y,g)
    print "One can write %s = (%s)(%s) + (%s)(%s)."%(g,x,a,y,b)
///
}}}

{{{id=30|
gcd(12780049, 1673450759)
///
}}}

{{{id=27|
@interact
def euclidean_alg(x=input_box(215441), y=input_box(46189)):
    if x < y:
        tmp = y
        y = x
        x = tmp
    orig_x = x
    orig_y = y
    n = len(x.str())                    
    print "Performing Euclidean algorithm with %s and %s."%(x,y)
    print
    r = y
    while r:
        q, r = x.quo_rem(y)
        if r:
            print ' ' * (n - len(x.str())), "%s = (%s)(%s) + %s"%(x, q, y, r)
        else:
            g = y
            print ' ' * (n - len(x.str())), "%s = (%s)(%s)"%(x, q, y)
        x = y
        y = r
    print 
    print "The gcd of %s and %s is %s."%(orig_x, orig_y, g)
///
}}}

{{{id=22|

///
}}}

<p>Similarly, Sage has a number of useful functions for using the Chinese Remainder Theorem (CRT). Recall that the CRT states that given a collection of integers $n_1$, $n_2$, $\ldots$, $n_k$, with $(n_i,n_j) = 1$ for all $i\neq j$, we have an isomorphism</p>
<p>$$ \mathbb{Z}/n \simeq \mathbb{Z}/n_1 \times \mathbb{Z}/n_2 \times \cdots \times \mathbb{Z}/n_k $$</p>
<p>where $n = n_1\cdot n_2\cdot \cdots \cdot n_k$.</p>

{{{id=21|
CRT(1,4,3,5)
///
}}}

{{{id=20|
19 % 3, 19 % 5
///
}}}

{{{id=19|
CRT_list([1,2,3],[2,3,5])
///
}}}

<p>Of course, usually one wants to explicitly construct the map from the product of the $\mathbb{Z}/n_i$ back to $\mathbb{Z}/n$. The command <span style="font-family: courier new,courier;">CRT_basis</span> accomplishes this.</p>

{{{id=18|
CRT_basis([2,3,5])
///
}}}

{{{id=17|
(15 * 1 + 10 * 2 + 6 * 3) % (2*3*5)
///
}}}

{{{id=16|
@interact
def crt_test(moduli=input_box([2,3,5]), residues=input_box([1,2,4])):
    print "The map from", ' '.join(['Z/%s x'%n for n in moduli[:-1]]), 'Z/%s'%moduli[-1],
    print "to Z/%s is given by:"%prod(moduli)
    
    crt_ls = CRT_basis(moduli)
    n = len(moduli)
    
    print
    print "(" + ', '.join(['x_%s'%i for i in range(1,n+1)]) + ') |--->',
    print ' + '.join(['%s*x_%s'%(crt_ls[i],i+1) for i in range(n)])

    print
    print "Under this map, (" + ', '.join([str(x) for x in residues]) + ") maps to",
    print "%s."%sum([x*y for x,y in zip(residues, crt_ls)])%prod(moduli)
///
}}}

{{{id=15|

///
}}}

<h2 style="text-align: center;">Modular Arithmetic</h2>

{{{id=14|
F = Integers(101)
F
///
}}}

{{{id=13|
F(105)
///
}}}

{{{id=12|
F(1)*102
///
}}}

{{{id=11|
F(50) + 100
///
}}}

{{{id=95|
F.multiplicative_group_is_cyclic()
///
}}}

{{{id=10|
F.unit_gens()
///
}}}

{{{id=9|
F.unit_group_exponent()
///
}}}

{{{id=132|
euler_phi(101)
///
}}}

{{{id=96|
F.quadratic_nonresidue()
///
}}}

{{{id=108|
GF(101)
///
}}}

{{{id=110|
F2.<a> = GF(101^2)
///
}}}

{{{id=111|
a.minpoly()
///
}}}

{{{id=112|
a.trace()
///
}}}

{{{id=107|
a.norm()
///
}}}

{{{id=8|

///
}}}

{{{id=56|
R = Integers(800)
R
///
}}}

{{{id=94|
R.multiplicative_group_is_cyclic()
///
}}}

{{{id=58|
factor(800)
///
}}}

{{{id=57|
R.unit_gens()
///
}}}

{{{id=59|
[ x.multiplicative_order() for x in R.unit_gens() ]
///
}}}

{{{id=136|
prod([ x.multiplicative_order() for x in R.unit_gens() ])
///
}}}

{{{id=135|
euler_phi(800)
///
}}}

{{{id=55|
R.unit_group_exponent()
///
}}}

{{{id=102|

///
}}}

<h2 style="text-align: center;">Quadratic Reciprocity</h2>

{{{id=101|
F = GF(7)
///
}}}

{{{id=100|
F(3).is_square()
///
}}}

{{{id=99|
F(2).is_square()
///
}}}

{{{id=118|
[ (x, x.is_square()) for x in F ]
///
}}}

{{{id=120|

///
}}}

{{{id=119|
F = GF(691)
///
}}}

{{{id=123|
F(2).is_square()
///
}}}

{{{id=122|
F(-2).is_square()
///
}}}

{{{id=117|

///
}}}

{{{id=131|
kronecker_symbol(691,2)
///
}}}

{{{id=130|
kronecker_symbol(691,-2)
///
}}}

{{{id=129|
kronecker_symbol(691,-1)
///
}}}

{{{id=128|

///
}}}

<p>Here the "upside-down Kronecker character" for $p$ is the character that sends $$x\mapsto\left(\frac{x}{p}\right).$$ One usually abbreviates this to $$\left(\frac{\cdot}{p}\right)$$ hence "upside-down."</p>

{{{id=116|
chi = kronecker_character_upside_down(691)
///
}}}

{{{id=115|
chi(2)
///
}}}

{{{id=125|
chi(-2)
///
}}}

{{{id=124|
chi(-1)
///
}}}

{{{id=53|

///
}}}

<h2 style="text-align: center;">Dirichlet Characters</h2>

<p>Recall that a Dirichlet character modulo an integer $N$ is simply a character of the group $(\mathbb{Z}/N)^\times$.</p>

{{{id=93|
G = DirichletGroup(49)
G
///
}}}

{{{id=92|
euler_phi(49)
///
}}}

{{{id=91|
G.gens()
///
}}}

{{{id=90|
chi = G.gens()[0]
chi
///
}}}

{{{id=89|
chi(17)
///
}}}

<p>Here the Gauss sum of a Dirichlet character $\chi$ of conductor $n$ is given by</p>
<p>$$ G(\chi) = \sum_{k=0}^n \chi(k)\zeta_n^k $$</p>
<p>for a fixed choice of a primitive $n^{\text{th}}$ root of unity $\zeta_n$ (e.g. $\zeta_n = e^{2\pi i / n}$).</p>

{{{id=88|
chi.gauss_sum()
///
}}}

{{{id=87|
lcm(49,euler_phi(49))
///
}}}

{{{id=85|
chi.gauss_sum_numerical()
///
}}}

{{{id=86|
chi.conductor()
///
}}}

{{{id=84|
(chi^7).conductor()
///
}}}

{{{id=83|
psi = chi^7
psi
///
}}}

{{{id=81|
psi.minimize_base_ring()
///
}}}

{{{id=174|
psi.bernoulli(1)
///
}}}

{{{id=80|

///
}}}

<h2 style="text-align: center;">Fibonacci Numbers</h2>

{{{id=143|
fibonacci(100)
///
}}}

{{{id=142|
list(fibonacci_sequence(15))
///
}}}

{{{id=141|
time n = fibonacci(5*10^6)
///
}}}

{{{id=140|
time len(n.digits())
///
}}}

{{{id=160|

///
}}}

{{{id=159|
fib = lambda n: reduce(lambda x,y: (x[1], x[2], x[1]+x[2]), xrange(n), (-1,1,0))[2]
///
}}}

{{{id=158|
[ fib(i) for i in range(15) ]
///
}}}

{{{id=139|

///
}}}

<h2 style="text-align: center;">Sloane's OEIS</h2>

{{{id=168|
sloane_find([2,3,5,7,11,13])
///
}}}

{{{id=167|
len(sloane_find([2,3,5,7,11,13]))
///
}}}

{{{id=166|
n = 7.38905609893065
ls = list(str(n)) ; ls.remove('.')
ls
///
}}}

{{{id=207|
sloane_find(ls)
///
}}}

{{{id=206|
@interact
def use_sloane(n=input_box(pi)):
    ls = [ x for x in list(str(n.n())) if x.isdigit() ]
    while ls[-1] == '0':
        ls.pop()
    print "Searching Sloane database for %s ..."%ls
    print
    print sloane_find(ls)
///
}}}

{{{id=210|

///
}}}

<h2 style="text-align: center;">Continued Fractions</h2>

{{{id=173|
cf = continued_fraction(100/691)
///
}}}

{{{id=172|
cf
///
}}}

{{{id=272|
QQ(cf)
///
}}}

{{{id=171|
show(cf)
///
}}}

{{{id=170|
cf.convergents()
///
}}}

{{{id=169|
cf_e = continued_fraction(e)
///
}}}

{{{id=186|
show(cf_e)
///
}}}

{{{id=191|
cf_e.convergents()
///
}}}

{{{id=185|
cf_pi = continued_fraction(pi)
///
}}}

{{{id=184|
cf_pi.convergents()
///
}}}

{{{id=190|
cf_pi.convergent(1)
///
}}}

{{{id=189|
cf_pi.convergent(3)
///
}}}

{{{id=188|
cf_pi.convergent(100)
///
}}}

{{{id=192|

///
}}}

{{{id=187|
cf_pi = continued_fraction(pi,bits=1000)
///
}}}

{{{id=183|
cf_pi.convergent(100)
///
}}}

{{{id=182|
@interact
def cf_conv(number=input_box(pi),convergent=input_box(10)):
    curr_prec=100
    cf = continued_fraction(number, bits=curr_prec)
    while len(cf) < convergent+20:
        curr_prec += 50
        cf = continued_fraction(number, bits=curr_prec)
    x = cf.convergent(convergent)
    print "The %sth convergent to %s is %s."%(convergent, number, x)
    print "This is approximately %s, which differs from %s by %s."%(x.n(), number, abs(number-x).n())
///
}}}

<p>One knows that the convergents $\frac{p_n}{q_n}$ of the continued fraction for $x$ satisfy</p>
<p>$$ \left| x - \frac{p_n}{q_n} \right| &lt; \frac{1}{q_n^2}. $$</p>
<p>This means that a graph of the differences $x-\frac{p_n}{q_n}$ will simply be a line along the $x$-axis; instead, we graph the differences divided by $q_n^2$.</p>

{{{id=201|
@interact
def f(number=input_box(pi), start=input_box(20), stop=input_box(100)):
    curr_prec = 100
    cf = continued_fraction(number, bits=curr_prec)
    n = len(cf)-10
    while n < stop:
        curr_prec += 50
        cf = continued_fraction(number, bits=curr_prec)
        n = len(cf)-10
    conv_ls = cf.convergents()
    ls = [(i,0) for i in range(start)]
    for i in range(start,stop):
        x = conv_ls[i]
        ls.append((i,((x.denominator()^2)*(number-x)).n(curr_prec)))
    show(list_plot(ls,plotjoined=True))
///
}}}

{{{id=193|

///
}}}

{{{id=137|

///
}}}

<h2 style="text-align: center;">Pythagorean Triples</h2>

{{{id=217|
@interact 
def __(t=('slope',1/2)): 
   t = QQ(t) 
   x = (1-t^2)/(1+t^2) 
   y = 2*t/(1+t^2) 
   r = t.numerator() 
   s = t.denominator() 
   a = s^2 - r^2; b = 2*r*s; c = s^2 + r^2 
   html('<h2 align=center>Point (x,y) = $%s$'%latex((x,y))) 
   html('Pythagorean (a,b,c) = $%s$</h2>'%latex((a,b,c))) 
   G = circle((0,0),1, rgbcolor='blue', thickness=3) 
   G += point([(-1,0), (x,y)], pointsize=150, rgbcolor='black') 
   G += line([(-1,0), (x,y)], rgbcolor='red', thickness=3) 
   G += text("(-1,0)", (-1.25,0.15), rgbcolor='black',fontsize=12) 
   try: 
      G.save('a.png',aspect_ratio=1) 
   except RuntimeError, msg: 
       print msg 
   html('<img src="cell://a.png">')
///
}}}

{{{id=216|

///
}}}

<h2 style="text-align: center;">Elliptic Curves</h2>

{{{id=225|
E = EllipticCurve([-4,4])
E
///
}}}

{{{id=224|
E.conductor()
///
}}}

{{{id=223|
E.rank()
///
}}}

{{{id=222|
L_E = E.lseries()
L_E
///
}}}

{{{id=221|
L_E.at1()
///
}}}

{{{id=234|
L_E.deriv_at1()
///
}}}

{{{id=233|
show(plot(E,xmax=10))
///
}}}

{{{id=238|
E_79 = E.change_ring(GF(79))
E_79
///
}}}

{{{id=237|
show(plot(E_79))
///
}}}

{{{id=236|
@interact
def f(E=input_box('11a'),p=input_box(97)):
    try:
        EE = EllipticCurve(E)
    except:
        print "Data (=%s) does not define an elliptic curve."%E
        return
        
    if not p.is_prime():
        print "%s is not prime."%p
        
    if p.divides(EE.conductor()): 
        print "Prime divides the conductor of E, please choose another."
        return
        
    equation = EE._repr_()
    start = equation.find('by ')
    end = equation.find(' over')
    equation = equation[start+3:end]
        
    print "Plotting elliptic curve defined by %s over F_%s ..."%(equation, p)
    print "E has %s points over F_%s."%(p+1-EE.ap(p),p)
    Ebar = EE.change_ring(GF(p))
    show(plot(Ebar))
///
}}}

{{{id=235|
@interact

def _(E=input_box([-1,0]),P=input_box([1,0]),n=input_box(10), xmax=input_box(6)):
    
    EE = EllipticCurve(E)
    try:
        PP = EE(P)
    except TypeError:
        print "P is not a point on E!"
        return
    
    print "E: ", EE
    print "P: ", PP.xy()
    print "rank(E): ", EE.rank()

    n_real = RDF(n)
    pts = [ k*PP for k in [1..n] ]
    
    img = plot(EE, xmax=xmax, rgbcolor=hue(0))
    for k in [1..n]:
        if not pts[k-1].is_zero() and pts[k-1][2] and pts[k-1].xy()[0] < xmax:
            img += plot(pts[k-1], pointsize=40, rgbcolor=hue(k/n_real))

    show(img)
///
}}}

{{{id=228|

///
}}}

<h2 style="text-align: center;">Quaternion Algebras</h2>

{{{id=73|
H = QuaternionAlgebra(-1, -1)
H
///
}}}

{{{id=70|
H.gens()
///
}}}

{{{id=69|
i,j,k = H.gens()
///
}}}

{{{id=68|
i^2
///
}}}

{{{id=67|
H.random_element()
///
}}}

{{{id=66|
x = H.random_element()
y = H.random_element()
timeit('x*y')
///
}}}

{{{id=65|
x = QQ.random_element()
y = QQ.random_element()
timeit('x*y')
///
}}}

{{{id=64|
H.discriminant()
///
}}}

{{{id=63|
H.is_division_algebra()
///
}}}

{{{id=62|
B = QuaternionAlgebra(2*3*5)
B
///
}}}

{{{id=60|
B.discriminant().factor()
///
}}}

{{{id=239|
x = B.random_element()
y = B.random_element()
timeit('x*y')
print
print "(%s) * (%s) = %s"%(x,y,x*y)
///
}}}

{{{id=240|
x = polygen(ZZ)
K.<a> = NumberField(x^4-17)
///
}}}

{{{id=241|
B = QuaternionAlgebra(K, a, a+1)
B
///
}}}

{{{id=242|
x = B.random_element()
y = B.random_element()
timeit('x*y')
print
print "(%s) * (%s) = %s"%(x,y,x*y)
///
}}}

{{{id=243|
x = K.random_element()
y = K.random_element()
timeit('x*y')
print
print "(%s) * (%s) = %s"%(x,y,x*y)
///
}}}

{{{id=248|

///
}}}

<h2 style="text-align: center;">So Much More ...</h2>
<p>There's so much more that Sage can do that we haven't even had a chance to look into. Here are just a few more bits and pieces ...</p>

{{{id=246|
ModularForms(37,2).q_expansion_basis(20)
///
}}}

{{{id=276|
eisenstein_series_qexp(4,30)
///
}}}

{{{id=274|
delta_qexp(30)
///
}}}

{{{id=275|
J = J0(33)
J
///
}}}

{{{id=278|
J.decomposition()
///
}}}

{{{id=280|
A = J[2]
A.cuspidal_subgroup()
///
}}}

{{{id=282|

///
}}}

{{{id=281|
x = polygen(ZZ)
K.<a> = NumberField(x^5-11*x^2+3)
K
///
}}}

{{{id=286|
O_K = K.ring_of_integers()
O_K.basis()
///
}}}

{{{id=284|
K.factor(7)
///
}}}

{{{id=290|
p = K.factor(7)[0][0]
p
///
}}}

{{{id=295|
p.ramification_index()
///
}}}

{{{id=294|
p.norm()
///
}}}

{{{id=293|
K.galois_group()
///
}}}

{{{id=300|

///
}}}

{{{id=299|
R = Zp(5)
R
///
}}}

{{{id=298|
R(1/3)
///
}}}

{{{id=297|
R(1/3)*3
///
}}}

{{{id=296|
R.extension(x^2-5,names='a')
///
}}}

{{{id=292|

///
}}}

{{{id=302|

///
}}}

{{{id=291|

///
}}}