Differences between revisions 8 and 68 (spanning 60 versions)
Revision 8 as of 2008-11-14 13:42:08
Size: 26507
Editor: anonymous
Comment: converted to 1.6 markup
Revision 68 as of 2020-06-14 09:10:48
Size: 48543
Editor: chapoton
Comment: fixone
Deletions are marked like this. Additions are marked like this.
Line 1: Line 1:
= Sage Interactions - Number Theory =
goto [[interact|interact main page]]
Line 5: Line 2:

= Integer Factorization =

== Divisibility Poset ==
by William Stein
{{{#!sagecell
@interact
def _(n=(5..100)):
    Poset(([1..n], lambda x, y: y%x == 0) ).show()
}}}

{{attachment:divposet.png}}
Line 8: Line 18:
{{{ {{{#!sagecell
Line 41: Line 51:
                    g += line([(j*2-len(cur),-i), ((k*2)-len(rows[i-1]),-i+1)],                      g += line([(j*2-len(cur),-i), ((k*2)-len(rows[i-1]),-i+1)],
Line 55: Line 65:
== Continued Fraction Plotter ==
by William Stein
{{{
@interact
def _(number=e, ymax=selector([None,5,20,..,400],nrows=2), clr=Color('purple'), prec=[500,1000,..,5000]):
    c = list(continued_fraction(RealField(prec)(number))); print c
    show(line([(i,z) for i, z in enumerate(c)],rgbcolor=clr),ymax=ymax,figsize=[10,2])
}}}
{{attachment:contfracplot.png}}
More complicated demonstration using Mathematica: http://demonstrations.wolfram.com/FactorTrees/

== Factoring an Integer ==
by Timothy Clemans

Sage implementation of the Mathematica demonstration of the same name. http://demonstrations.wolfram.com/FactoringAnInteger/

{{{#!sagecell
@interact
def _(r=selector(range(0,10000,1000), label='range', buttons=True), n=slider(0,1000,1,2,'n',False)):
    if not r and n in (0, 1):
        n = 2
    s = '$%d = %s$' % (r + n, factor(r + n))
    s = s.replace('*', '\\times')
    pretty_print(html(s))
}}}

= Prime Numbers =
Line 67: Line 86:
{{{
@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'))
{{{#!sagecell
@interact
def _(N=(100,list(range(2,2000)))):
    pretty_print(html(r"<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, color='red') + plot(x/(log(x)-1), 5, N, color='blue'))
Line 75: Line 94:
== Computing Generalized Bernoulli Numbers ==
by William Stein (Sage-2.10.3)
{{{
@interact
def _(m=selector([1..15],nrows=2), n=(7,(3..10))):
    G = DirichletGroup(m)
    s = "<h3>First n=%s Bernoulli numbers attached to characters with modulus m=%s</h3>"%(n,m)
    s += '<table border=1>'
    s += '<tr bgcolor="#edcc9c"><td align=center>$\\chi$</td><td>Conductor</td>' + \
           ''.join('<td>$B_{%s,\chi}$</td>'%k for k in [1..n]) + '</tr>'
    for eps in G.list():
        v = ''.join(['<td align=center bgcolor="#efe5cd">$%s$</td>'%latex(eps.bernoulli(k)) for k in [1..n]])
        s += '<tr><td bgcolor="#edcc9c">%s</td><td bgcolor="#efe5cd" align=center>%s</td>%s</tr>\n'%(
             eps, eps.conductor(), v)
    s += '</table>'
    html(s)
}}}

{{attachment:bernoulli.png}}


== Fundamental Domains of SL_2(ZZ) ==
by Robert Miller
{{{
L = [[-0.5, 2.0^(x/100.0) - 1 + sqrt(3.0)/2] for x in xrange(1000, -1, -1)]
R = [[0.5, 2.0^(x/100.0) - 1 + sqrt(3.0)/2] for x in xrange(1000)]
xes = [x/1000.0 for x in xrange(-500,501,1)]
M = [[x,abs(sqrt(x^2-1))] for x in xes]
fundamental_domain = L+M+R
fundamental_domain = [[x-1,y] for x,y in fundamental_domain]
@interact
def _(gen = selector(['t+1', 't-1', '-1/t'], nrows=1)):
    global fundamental_domain
    if gen == 't+1':
        fundamental_domain = [[x+1,y] for x,y in fundamental_domain]
    elif gen == 't-1':
        fundamental_domain = [[x-1,y] for x,y in fundamental_domain]
    elif gen == '-1/t':
        new_dom = []
        for x,y in fundamental_domain:
            sq_mod = x^2 + y^2
            new_dom.append([(-1)*x/sq_mod, y/sq_mod])
        fundamental_domain = new_dom
    P = polygon(fundamental_domain)
    P.ymax(1.2); P.ymin(-0.1)
    P.show()
}}}

{{attachment:fund_domain.png}}

== Computing modular forms ==
by William Stein
{{{
j = 0
@interact
def _(N=[1..100], k=selector([2,4,..,12],nrows=1), prec=(3..40),
      group=[(Gamma0, 'Gamma0'), (Gamma1, 'Gamma1')]):
    M = CuspForms(group(N),k)
    print j; global j; j += 1
    print M; print '\n'*3
    print "Computing basis...\n\n"
    if M.dimension() == 0:
         print "Space has dimension 0"
    else:
        prec = max(prec, M.dimension()+1)
        for f in M.basis():
             view(f.q_expansion(prec))
    print "\n\n\nDone computing basis."
}}}

{{attachment:modformbasis.png}}


== Computing the cuspidal subgroup ==
by William Stein
{{{
html('<h1>Cuspidal Subgroups of Modular Jacobians J0(N)</h1>')
@interact
def _(N=selector([1..8*13], ncols=8, width=10, default=10)):
    A = J0(N)
    print A.cuspidal_subgroup()
}}}

{{attachment:cuspgroup.png}}

== A Charpoly and Hecke Operator Graph ==
by William Stein

{{{
# Note -- in Sage-2.10.3; multiedges are missing in plots; loops are missing in 3d plots
@interact
def f(N = prime_range(11,400),
      p = selector(prime_range(2,12),nrows=1),
      three_d = ("Three Dimensional", False)):
    S = SupersingularModule(N)
    T = S.hecke_matrix(p)
    G = Graph(T, multiedges=True, loops=not three_d)
    html("<h1>Charpoly and Hecke Graph: Level %s, T_%s</h1>"%(N,p))
    show(T.charpoly().factor())
    if three_d:
        show(G.plot3d(), aspect_ratio=[1,1,1])
    else:
        show(G.plot(),figsize=7)
}}}

{{attachment:heckegraph.png}}

== Demonstrating the Diffie-Hellman Key Exchange Protocol ==
by Timothy Clemans (refereed by William Stein)
{{{
@interact
def diffie_hellman(button=selector(["New example"],label='',buttons=True),
    bits=("Number of bits of prime", (8,12,..512))):
    maxp = 2^bits
    p = random_prime(maxp)
    k = GF(p)
    if bits>100:
        g = k(2)
    else:
        g = k.multiplicative_generator()
    a = ZZ.random_element(10, maxp)
    b = ZZ.random_element(10, maxp)

    print """
<html>
<style>
.gamodp {
background:yellow
}
.gbmodp {
background:orange
}
.dhsame {
color:green;
font-weight:bold
}
</style>
<h2>%s-Bit Diffie-Hellman Key Exchange</h2>
<ol style="color:#000;font:12px Arial, Helvetica, sans-serif">
<li>Alice and Bob agree to use the prime number p=%s and base g=%s.</li>
<li>Alice chooses the secret integer a=%s, then sends Bob (<span class="gamodp">g<sup>a</sup> mod p</span>):<br/>%s<sup>%s</sup> mod %s = <span class="gamodp">%s</span>.</li>
<li>Bob chooses the secret integer b=%s, then sends Alice (<span class="gbmodp">g<sup>b</sup> mod p</span>):<br/>%s<sup>%s</sup> mod %s = <span class="gbmodp">%s</span>.</li>
<li>Alice computes (<span class="gbmodp">g<sup>b</sup> mod p</span>)<sup>a</sup> mod p:<br/>%s<sup>%s</sup> mod %s = <span class="dhsame">%s</span>.</li>
<li>Bob computes (<span class="gamodp">g<sup>a</sup> mod p</span>)<sup>b</sup> mod p:<br/>%s<sup>%s</sup> mod %s = <span class="dhsame">%s</span>.</li>
</ol></html>
    """ % (bits, p, g, a, g, a, p, (g^a), b, g, b, p, (g^b), (g^b), a, p,
       (g^ b)^a, g^a, b, p, (g^a)^b)
}}}

{{attachment:dh.png}}

== Plotting an elliptic curve over a finite field ==
{{{
E = EllipticCurve('37a')
@interact
def _(p=slider(prime_range(1000), default=389)):
    show(E)
    print "p = %s"%p
    show(E.change_ring(GF(p)).plot(),xmin=0,ymin=0)
}}}

{{attachment:ellffplot.png}}

== Prime Spiral - Square ==
== Prime Spiral - Square FIXME ==
Line 240: Line 96:
{{{ {{{#!sagecell
Line 245: Line 101:
    REFERENCES:      REFERENCES:
Line 250: Line 106:
        Weisstein, Eric W. "Prime-Generating Polynomial." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Prime-GeneratingPolynomial.html          Weisstein, Eric W. "Prime-Generating Polynomial." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Prime-GeneratingPolynomial.html
Line 258: Line 114:
        elif y<0 and -x >= y and y<x: return 4*(y+1)^2 -11*(y+1) + (start+7) +x 
        else: print 'NaN'

    #Takes in an n and the start value of the spiral and gives its (x,y) coordinate 
        elif y<0 and -x >= y and y<x: return 4*(y+1)^2 -11*(y+1) + (start+7) +x
        else: print('NaN')

    #Takes in an n and the start value of the spiral and gives its (x,y) coordinate
Line 263: Line 119:
        num = num - start +1          num = num - start +1
Line 265: Line 121:
        top = ceil(sqrt(num))             top = ceil(sqrt(num))
Line 270: Line 126:
            else:              else:
Line 277: Line 133:
            else:              else:
Line 284: Line 140:
    if start < 1 or end <=start: print "invalid start or end value"
    if n > end: print "WARNING: n is larger than the end value"
    if start < 1 or end <=start: print("invalid start or end value")
    if n > end: print("WARNING: n is larger than the end value")
Line 289: Line 145:
        N = M.copy()         N = copy(M)
Line 293: Line 149:
 
    #These functions return an int based on where the t is located in the spiral 

    #These functions return an int based on where the t is located in the spiral
Line 306: Line 162:
    if n !=0: x_cord, y_cord = find_xy(n, start) #Overrides the user given x and y coordinates      if n !=0: x_cord, y_cord = find_xy(n, start) #Overrides the user given x and y coordinates
Line 314: Line 170:
    
Line 316: Line 172:
    #print x_cord, y_cord
if show_lines: 
        for t in [(-size-1)..size+1]: 
    if show_lines:
        for t in [(-size-1)..size+1]:
Line 320: Line 175:
            if m.is_pseudoprime(): main_list.add(m)              if m.is_pseudoprime(): main_list.add(m)
Line 325: Line 180:
    #This for loop changes the matrix by spiraling out from the center and changing each entry as it goes. It is faster than the find_xy function above.      #This for loop changes the matrix by spiraling out from the center and changing each entry as it goes. It is faster than the find_xy function above.
Line 327: Line 182:
        #print x, "=x y=", y, " num =", num
Line 330: Line 184:
            else: x-=1              else: x-=1
Line 332: Line 186:
        
        elif county < overcount: 

        elif county < overcount:
Line 335: Line 189:
            else: y-=1              else: y-=1
Line 337: Line 191:
        else:          else:
Line 343: Line 197:
    
        if not invert and num in main_list: 

        if not invert and num in main_list:
Line 351: Line 205:
    
    if n != 0: 
        print '(to go from x,y coords to an n, reset by setting n=0)'

    if n != 0:
        print('(to go from x,y coords to an n, reset by setting n=0)')
Line 355: Line 209:
        #print 'if n =', n, 'then (x,y) =', (x_cord, y_cord)

    print
'(x,y) =', (x_cord, y_cord), '<=> n =', find_n(x_cord, y_cord, start)
    print ' '
    print
"SW/NE line"
    if -y_cord<x_cord: print '4*t^2 + 2*t +', -x_cord+y_cord+start
    else: print '4*t^2 + 2*t +', +x_cord-y_cord+start

    print
"NW/SE line"
    if x_cord<y_cord: print '4*t^2 +', -x_cord-y_cord+start
    else: print '4*t^2 + 4*t +', +x_cord+y_cord+start

    print(
'(x,y) =', (x_cord, y_cord), '<=> n =', find_n(x_cord, y_cord, start))
    print(' ')
    print(
"SW/NE line")
    if -y_cord<x_cord: print('4*t^2 + 2*t +', -x_cord+y_cord+start)
    else: print('4*t^2 + 2*t +', +x_cord-y_cord+start)

    print(
"NW/SE line")
    if x_cord<y_cord: print('4*t^2 +', -x_cord-y_cord+start)
    else: print('4*t^2 + 4*t +', +x_cord+y_cord+start)
Line 368: Line 221:
    else:      else:
Line 377: Line 230:
{{{
@interact
def polar_prime_spiral(start=1, end=2000, show_factors = false, highlight_primes = false, show_curves=true, n = 0): 

    #For more information about the factors in the spiral, visit http://www.dcs.gla.ac.uk/~jhw/spirals/index.html by John Williamson. 

    if start < 1 or end <=start: print "invalid start or end value"
    if n > end: print "WARNING: n is greater than end value"

Needs fix for show_factors
{{{#!sagecell
@interact
def polar_prime_spiral(start=1, end=2000, show_factors = false, highlight_primes = false, show_curves=true, n = 0):

    #For more information about the factors in the spiral, visit http://www.dcs.gla.ac.uk/~jhw/spirals/index.html by John Williamson.

    if start < 1 or end <=start: print("invalid start or end value")
    if n > end: print("WARNING: n is greater than end value")
Line 387: Line 242:
    
Line 395: Line 250:
        R = points(list2, alpha = .1) #Faded Composites 
    else: 
        R = points(list2, alpha = .1) #Faded Composites
    else:
Line 403: Line 258:
        R=points(list2, hue = .1, pointsize = p_size) 
    
        R=points(list2, hue = .1, pointsize = p_size)
Line 406: Line 261:
        print 'n =', factor(n)
        
        print('n = {}'.format(factor(n)))
Line 414: Line 269:
        Q = plot(W1+W2+W3+W4, alpha = .1)                   Q = plot(W1+W2+W3+W4, alpha = .1)
Line 417: Line 272:
        if show_curves:          if show_curves:
Line 422: Line 277:
            if n > (floor(sqrt(n)))^2 and n <= (floor(sqrt(n)))^2 + floor(sqrt(n)):              if n > (floor(sqrt(n)))^2 and n <= (floor(sqrt(n)))^2 + floor(sqrt(n)):
Line 425: Line 280:
            else:              else:
Line 428: Line 283:
            print 'Pink Curve: n^2 +', c
            print 'Green Curve: n^2 + n +', c2
            def g(m): return (a*m^2+b*m+c); 
            print('Pink Curve: n^2 +', c)
            print('Green Curve: n^2 + n +', c2)
            def g(m): return (a*m^2+b*m+c);
Line 436: Line 291:
            c= c2;              c= c2;
Line 446: Line 301:
== Quadratic Residue Table ==
= Modular Forms =

== Computing modular forms ==
by William Stein
{{{#!sagecell
@interact
def _(N=[1..100], k=selector([2,4,..,12],nrows=1), prec=(3..40),
      group=[(Gamma0, 'Gamma0'), (Gamma1, 'Gamma1')]):
    M = CuspForms(group(N),k)
    print(M)
    print('\n' * 3)
    print("Computing basis...\n\n")
    if M.dimension() == 0:
         print("Space has dimension 0")
    else:
        prec = max(prec, M.dimension() + 1)
        for f in M.basis():
             view(f.q_expansion(prec))
    print("\n\n\nDone computing basis.")
}}}

{{attachment:modformbasis.png}}


== Computing the cuspidal subgroup ==
by William Stein

ncols not working
{{{#!sagecell
pretty_print(html('<h1>Cuspidal Subgroups of Modular Jacobians J0(N)</h1>'))
@interact
def _(N=selector([1..8*13], ncols=8, width=10, default=10)):
    A = J0(N)
    print(A.cuspidal_subgroup())
}}}

{{attachment:cuspgroup.png}}

== A Charpoly and Hecke Operator Graph ==
by William Stein

{{{#!sagecell
# Note -- in Sage-2.10.3; multiedges are missing in plots; loops are missing in 3d plots
@interact
def f(N = prime_range(11,400),
      p = selector(prime_range(2,12),nrows=1),
      three_d = ("Three Dimensional", False)):
    S = SupersingularModule(N)
    T = S.hecke_matrix(p)
    G = DiGraph(T, multiedges=not three_d)
    if three_d:
        G.remove_loops()
    html("<h1>Charpoly and Hecke Graph: Level %s, T_%s</h1>"%(N,p))
    show(T.charpoly().factor())
    if three_d:
        show(G.plot3d(), aspect_ratio=[1,1,1])
    else:
        show(G.plot(),figsize=7)
}}}

{{attachment:heckegraph.png}}

= Modular Arithmetic =

== Quadratic Residue Table FIXME ==
Line 448: Line 368:
{{{ {{{#!sagecell
Line 499: Line 419:
== Cubic Residue Table == == Cubic Residue Table FIXME ==
Line 501: Line 421:
{{{ {{{#!sagecell
Line 519: Line 439:
    if Mod(a,3)!=0 and Mod(b,3)==0:
        return True
    else:
        return False
    return Mod(a,3)!=0 and Mod(b,3)==0
Line 557: Line 474:
        MP += line([(i,0),(i,r)], rgbcolor='black')          MP += line([(i,0),(i,r)], rgbcolor='black')
Line 560: Line 477:
                MP += text('$\omega^2$',(i+.5,r-j-.5),rgbcolor='black')                 MP += text(r'$\omega^2$',(i+.5,r-j-.5),rgbcolor='black')
Line 562: Line 479:
                MP += text('$\omega $',(i+.5,r-j-.5),rgbcolor='black')                 MP += text(r'$\omega $',(i+.5,r-j-.5),rgbcolor='black')
Line 571: Line 488:
    MP += text('$ \pi_1$',(r/2,r+2), rgbcolor='black', fontsize=25)
    MP += text('$ \pi_2$',(-2.5,r/2), rgbcolor='black', fontsize=25)

    html('Symmetry of Primary Cubic Residues mod ' \
          + '%d primary primes in $ \mathbf Z[\omega]$.'%r)
    MP += text(r'$ \pi_1$',(r/2,r+2), rgbcolor='black', fontsize=25)
    MP += text(r'$ \pi_2$',(-2.5,r/2), rgbcolor='black', fontsize=25)

    pretty_print(html('Symmetry of Primary Cubic Residues mod ' \
          + r'%d primary primes in $ \mathbf Z[\omega]$.'%r))
Line 580: Line 497:

= Cyclotomic Fields =
Line 583: Line 502:
{{{ {{{#!sagecell
Line 632: Line 551:
    S = circle((0,0),1,rgbcolor='yellow')  \
    +
line([e_pt,e_gs_pt], rgbcolor='red', thickness=4) \
    +
line([f_pt,f_gs_pt], rgbcolor='blue', thickness=3) \
    +
line([ef_pt,ef_gs_pt], rgbcolor='purple',thickness=2) \
    +
point(e_pt,pointsize=50, rgbcolor='red')  \
    +
point(f_pt,pointsize=50, rgbcolor='blue') \
    +
point(ef_pt,pointsize=50,rgbcolor='purple') \
    +
point(f_gs_pt,pointsize=75, rgbcolor='black') \
    +
point(e_gs_pt,pointsize=75, rgbcolor='black') \
    +
point(ef_gs_pt,pointsize=75, rgbcolor='black') \
    +
point(js_pt,pointsize=100,rgbcolor='green')
    S = circle((0,0),1,rgbcolor='yellow')
    S +=
line([e_pt,e_gs_pt], rgbcolor='red', thickness=4)
    S +=
line([f_pt,f_gs_pt], rgbcolor='blue', thickness=3)
    S +=
line([ef_pt,ef_gs_pt], rgbcolor='purple',thickness=2)
    S +=
point(e_pt,pointsize=50, rgbcolor='red')
    S +=
point(f_pt,pointsize=50, rgbcolor='blue')
    S +=
point(ef_pt,pointsize=50,rgbcolor='purple')
    S +=
point(f_gs_pt,pointsize=75, rgbcolor='black')
    S +=
point(e_gs_pt,pointsize=75, rgbcolor='black')
    S +=
point(ef_gs_pt,pointsize=75, rgbcolor='black')
    S +=
point(js_pt,pointsize=100,rgbcolor='green')
Line 644: Line 563:
        S += text('$J(%s,%s) = %s$'%(latex2(e),latex2(f),latex(js)), \         S += text('$J(%s,%s) = %s$'%(latex2(e),latex2(f),latex(js)),
Line 663: Line 582:
{{{ {{{#!sagecell
Line 712: Line 631:
    S = circle((0,0),1,rgbcolor='yellow')  \
    +
line([e_pt,e_gs_pt], rgbcolor='red', thickness=4) \
    +
line([f_pt,f_gs_pt], rgbcolor='blue', thickness=3) \
    +
line([ef_pt,ef_gs_pt], rgbcolor='purple',thickness=2) \
    +
point(e_pt,pointsize=50, rgbcolor='red')  \
    +
point(f_pt,pointsize=50, rgbcolor='blue') \
    +
point(ef_pt,pointsize=50,rgbcolor='purple') \
    +
point(f_gs_pt,pointsize=75, rgbcolor='black') \
    +
point(e_gs_pt,pointsize=75, rgbcolor='black') \
    +
point(ef_gs_pt,pointsize=75, rgbcolor='black') \
    +
point(js_pt,pointsize=100,rgbcolor='green')
    S = circle((0,0),1,rgbcolor='yellow')
    S +=
line([e_pt,e_gs_pt], rgbcolor='red', thickness=4)
    S +=
line([f_pt,f_gs_pt], rgbcolor='blue', thickness=3)
    S +=
line([ef_pt,ef_gs_pt], rgbcolor='purple',thickness=2)
    S +=
point(e_pt,pointsize=50, rgbcolor='red')
    S +=
point(f_pt,pointsize=50, rgbcolor='blue')
    S +=
point(ef_pt,pointsize=50,rgbcolor='purple')
    S +=
point(f_gs_pt,pointsize=75, rgbcolor='black')
    S +=
point(e_gs_pt,pointsize=75, rgbcolor='black')
    S +=
point(ef_gs_pt,pointsize=75, rgbcolor='black')
    S +=
point(js_pt,pointsize=100,rgbcolor='green')
Line 724: Line 643:
        S += text('$J(%s,%s) = %s$'%(latex2(e),latex2(f),latex(js)), \         S += text('$J(%s,%s) = %s$'%(latex2(e),latex2(f),latex(js)),
Line 727: Line 646:
        html('$$J(%s,%s) = %s$$'%(latex2(e),latex2(f),latex(js)))         pretty_print(html('$$J(%s,%s) = %s$$'%(latex2(e),latex2(f),latex(js))))
Line 736: Line 655:
        ga[i].save('j%d.PNG'%i,figsize=4,aspect_ratio=1, \         ga[i].save('j%d.png'%i,figsize=4,aspect_ratio=1,
Line 742: Line 661:
    html('<table bgcolor=lightgrey cellpadding=2>')     s='<table bgcolor=lightgrey cellpadding=2>'
Line 744: Line 663:
        html('<tr><td align="center"><img src="cell://j%d.PNG"></td>'%(2*i))
        html('<td align="center"><img src="cell://j%d.PNG"></td></tr>'%(2*i+1))
    html('</table>')
        s+='<tr><td align="center"><img src="cell://j%d.png"></td>'%(2*i)
        s+='<td align="center"><img src="cell://j%d.png"></td></tr>'%(2*i+1)
    s+='</table>'
    pretty_print(html(s
))
Line 750: Line 670:

= Elliptic Curves =

== Adding points on an elliptic curve ==
by David Møller Hansen
{{{#!sagecell
def point_txt(P,name,rgbcolor):
    if (P.xy()[1]) < 0:
        r = text(name,[float(P.xy()[0]),float(P.xy()[1])-1],rgbcolor=rgbcolor)
    elif P.xy()[1] == 0:
        r = text(name,[float(P.xy()[0]),float(P.xy()[1])+1],rgbcolor=rgbcolor)
    else:
        r = text(name,[float(P.xy()[0]),float(P.xy()[1])+1],rgbcolor=rgbcolor)
    return r

E = EllipticCurve('37a')
list_of_points = E.integral_points()
html("Graphical addition of two points $P$ and $Q$ on the curve $ E: %s $"%latex(E))

def line_from_curve_points(E,P,Q,style='-',rgb=(1,0,0),length=25):
 """
 P,Q two points on an elliptic curve.
 Output is a graphic representation of the straight line intersecting with P,Q.
 """
 # The function tangent to P=Q on E
 if P == Q:
  if P[2]==0:
   return line([(1,-length),(1,length)],linestyle=style,rgbcolor=rgb)
  else:
   # Compute slope of the curve E in P
   l=-(3*P[0]^2 + 2*E.a2()*P[0] + E.a4() - E.a1()*P[1])/((-2)*P[1] - E.a1()*P[0] - E.a3())
   f(x) = l * (x - P[0]) + P[1]
   return plot(f(x),-length,length,linestyle=style,rgbcolor=rgb)
 # Trivial case of P != R where P=O or R=O then we get the vertical line from the other point
 elif P[2] == 0:
  return line([(Q[0],-length),(Q[0],length)],linestyle=style,rgbcolor=rgb)
 elif Q[2] == 0:
  return line([(P[0],-length),(P[0],length)],linestyle=style,rgbcolor=rgb)
 # Non trivial case where P != R
 else:
  # Case where x_1 = x_2 return vertical line evaluated in Q
  if P[0] == Q[0]:
   return line([(P[0],-length),(P[0],length)],linestyle=style,rgbcolor=rgb)

  #Case where x_1 != x_2 return line trough P,R evaluated in Q"
  l=(Q[1]-P[1])/(Q[0]-P[0])
  f(x) = l * (x - P[0]) + P[1]
  return plot(f(x),-length,length,linestyle=style,rgbcolor=rgb)

@interact
def _(P=selector(list_of_points,label='Point P'),Q=selector(list_of_points,label='Point Q'), marked_points = checkbox(default=True,label = 'Points'), Lines = selector([0..2],nrows=1), Axes=True):
 curve = E.plot(rgbcolor = (0,0,1),xmin=-5,xmax=5,plot_points=300)
 R = P + Q
 Rneg = -R
 l1 = line_from_curve_points(E,P,Q)
 l2 = line_from_curve_points(E,R,Rneg,style='--')
 p1 = plot(P,rgbcolor=(1,0,0),pointsize=40)
 p2 = plot(Q,rgbcolor=(1,0,0),pointsize=40)
 p3 = plot(R,rgbcolor=(1,0,0),pointsize=40)
 p4 = plot(Rneg,rgbcolor=(1,0,0),pointsize=40)
 textp1 = point_txt(P,"$P$",rgbcolor=(0,0,0))
 textp2 = point_txt(Q,"$Q$",rgbcolor=(0,0,0))
 textp3 = point_txt(R,"$P+Q$",rgbcolor=(0,0,0))
 if Lines==0:
  g=curve
 elif Lines ==1:
  g=curve+l1
 elif Lines == 2:
  g=curve+l1+l2
 if marked_points:
  g=g+p1+p2+p3+p4
 if P != Q:
  g=g+textp1+textp2+textp3
 else:
  g=g+textp1+textp3
 g.axes_range(xmin=-5,xmax=5,ymin=-13,ymax=13)
 show(g,axes = Axes)
}}}
{{attachment:PointAddEllipticCurve.png}}


== Plotting an elliptic curve over a finite field ==
{{{#!sagecell
E = EllipticCurve('37a')
@interact
def _(p=slider(prime_range(1000), default=389)):
    show(E)
    print("p = %s" % p)
    show(E.change_ring(GF(p)).plot(), xmin=0, ymin=0)
}}}

{{attachment:ellffplot.png}}

= Cryptography =

== The Diffie-Hellman Key Exchange Protocol ==
by Timothy Clemans and William Stein
{{{#!sagecell
@interact
def diffie_hellman(bits=slider(8, 513, 4, 8, 'Number of bits', False),
    button=selector(["Show new example"],label='',buttons=True)):
    maxp = 2 ^ bits
    p = random_prime(maxp)
    k = GF(p)
    if bits > 100:
        g = k(2)
    else:
        g = k.multiplicative_generator()
    a = ZZ.random_element(10, maxp)
    b = ZZ.random_element(10, maxp)

    pretty_print(html("""
<style>
.gamodp, .gbmodp {
color:#000;
padding:5px
}
.gamodp {
background:#846FD8
}
.gbmodp {
background:#FFFC73
}
.dhsame {
color:#000;
font-weight:bold
}
</style>
<h2 style="color:#000;font-family:Arial, Helvetica, sans-serif">%s-Bit Diffie-Hellman Key Exchange</h2>
<ol style="color:#000;font-family:Arial, Helvetica, sans-serif">
<li>Alice and Bob agree to use the prime number p = %s and base g = %s.</li>
<li>Alice chooses the secret integer a = %s, then sends Bob (<span class="gamodp">g<sup>a</sup> mod p</span>):<br/>%s<sup>%s</sup> mod %s = <span class="gamodp">%s</span>.</li>
<li>Bob chooses the secret integer b=%s, then sends Alice (<span class="gbmodp">g<sup>b</sup> mod p</span>):<br/>%s<sup>%s</sup> mod %s = <span class="gbmodp">%s</span>.</li>
<li>Alice computes (<span class="gbmodp">g<sup>b</sup> mod p</span>)<sup>a</sup> mod p:<br/>%s<sup>%s</sup> mod %s = <span class="dhsame">%s</span>.</li>
<li>Bob computes (<span class="gamodp">g<sup>a</sup> mod p</span>)<sup>b</sup> mod p:<br/>%s<sup>%s</sup> mod %s = <span class="dhsame">%s</span>.</li>
</ol>
    """ % (bits, p, g, a, g, a, p, (g^a), b, g, b, p, (g^b), (g^b), a, p,
       (g^ b)^a, g^a, b, p, (g^a)^b)))
}}}


{{attachment:dh.png}}

= Other =

== Continued Fraction Plotter ==
by William Stein

crows not working
{{{#!sagecell
@interact
def _(number=e, ymax=selector([5,20,..,400],nrows=2), clr=Color('purple'), prec=[500,1000,..,5000]):
    c = list(continued_fraction(RealField(prec)(number))); print(c)
    show(line([(i,z) for i, z in enumerate(c)],rgbcolor=clr),ymax=ymax,figsize=[10,2])
}}}
{{attachment:contfracplot.png}}

== Computing Generalized Bernoulli Numbers ==
by William Stein (Sage-2.10.3)
{{{#!sagecell
@interact
def _(m=selector([1..15],nrows=2), n=(7,[3..10])):
    G = DirichletGroup(m)
    s = r"<h3>First n=%s Bernoulli numbers attached to characters with modulus m=%s</h3>"%(n,m)
    s += r'<table border=1>'
    s += r'<tr bgcolor="#edcc9c"><td align=center>$\chi$</td><td>Conductor</td>' + \
           ''.join(r'<td>$B_{%s,\chi}$</td>'%k for k in [1..n]) + '</tr>'
    for eps in G.list():
        v = ''.join(['<td align=center bgcolor="#efe5cd">$%s$</td>'%latex(eps.bernoulli(k)) for k in [1..n]])
        s += '<tr><td bgcolor="#edcc9c">%s</td><td bgcolor="#efe5cd" align=center>%s</td>%s</tr>\n'%(
             eps, eps.conductor(), v)
    s += '</table>'
    pretty_print(html(s))
}}}

{{attachment:bernoulli.png}}


== Fundamental Domains of SL_2(ZZ) ==
by Robert Miller
{{{#!sagecell
L = [[-0.5, 2.0^(x/100.0) - 1 + sqrt(3.0)/2] for x in range(1000, -1, -1)]
R = [[0.5, 2.0^(x/100.0) - 1 + sqrt(3.0)/2] for x in range(1000)]
xes = [x/1000.0 for x in range(-500,501,1)]
M = [[x,abs(sqrt(x^2-1))] for x in xes]
fundamental_domain = L+M+R
fundamental_domain = [[x-1,y] for x,y in fundamental_domain]
@interact
def _(gen = selector(['t+1', 't-1', '-1/t'], buttons=True,nrows=1)):
    global fundamental_domain
    if gen == 't+1':
        fundamental_domain = [[x+1,y] for x,y in fundamental_domain]
    elif gen == 't-1':
        fundamental_domain = [[x-1,y] for x,y in fundamental_domain]
    elif gen == '-1/t':
        new_dom = []
        for x,y in fundamental_domain:
            sq_mod = x^2 + y^2
            new_dom.append([(-1)*x/sq_mod, y/sq_mod])
        fundamental_domain = new_dom
    P = polygon(fundamental_domain)
    P.ymax(1.2); P.ymin(-0.1)
    P.show()
}}}

{{attachment:fund_domain.png}}

= Multiple Zeta Values =
by Akhilesh P.
== Computing Multiple Zeta values ==
=== Word Input ===
{{{#!sagecell
R=RealField(10)
@interact
def _( weight=(5,(2..100))):
 n=weight
 a=[0 for i in range(n-1)]
 a.append(1)
 @interact
 def _(v=('word', input_grid(1, n, default=[a], to_value=lambda x: vector(flatten(x)))), accuracy=(100..100000)):
  D=accuracy
  a=[v[i] for i in range(len(v))]
  DD=int(3.321928*D)+int(R(log(3.321928*D))/R(log(10)))+4
  RIF=RealIntervalField(DD)
  def Li(word):
        n=int(DD*log(10)/log(2))+1
        B=[]
        L=[]
        S=[]
        count=-1
        k=len(word)
        for i in range(k):
                B.append(RIF('0'))
                L.append(RIF('0'))
                if(word[i]==1 and i<k-1):
                        S.append(RIF('0'))
                        count=count+1
        T=RIF('1')
        for m in range(n):
                T=T/2
                B[k-1]=RIF('1')/(m+1)
                j=count
                for i in range(k-2,-1,-1):
                        if(word[i]==0):
                                B[i]=B[i+1]/(m+1)
                        elif(word[i]==1):
                                B[i]=S[j]/(m+1)
                                S[j]=S[j]+B[i+1]
                                j=j-1
                        L[i]=T*B[i]+L[i]
                L[k-1]=T*B[k-1]+L[k-1]
        return(L)
  def dual(a):
        b=list()
        b=a
        b=b[::-1]
        for i in range(len(b)):
                b[i]=1-b[i]
        return(b)
  def zeta(a):
        b=dual(a)
        l1=Li(a)+[1]
        l2=Li(b)+[1]
        Z=RIF('0')
        for i in range(len(l1)):
                Z=Z+l1[i]*l2[len(a)-i]
        return(Z)
  u=zeta(a)
  RIF=RealIntervalField(int(3.321928*D))
  u=u/1
  print(u)
}}}
{{attachment:akhi1.png}}
=== Composition Input ===
{{{#!sagecell
R=RealField(10)
@interact
def _( Depth=(5,(2..100))):
 n=Depth
 a=[2]
 a=a+[1 for i in range(n-1)]
 @interact
 def _(v=('Composition', input_grid(1, n, default=[a], to_value=lambda x: vector(flatten(x)))), accuracy=(100..100000)):
  D=accuracy
  a=[v[i] for i in range(len(v))]
  def comptobin(a):
        word=[]
        for i in range(len(a)):
                word=word+[0]*(a[i]-1)+[1]
        return(word)
  a=comptobin(a)
  DD=int(3.321928*D)+int(R(log(3.321928*D))/R(log(10)))+4
  RIF=RealIntervalField(DD)
  def Li(word):
        n=int(DD*log(10)/log(2))+1
        B=[]
        L=[]
        S=[]
        count=-1
        k=len(word)
        for i in range(k):
                B.append(RIF('0'))
                L.append(RIF('0'))
                if(word[i]==1 and i<k-1):
                        S.append(RIF('0'))
                        count=count+1
        T=RIF('1')
        for m in range(n):
                T=T/2
                B[k-1]=RIF('1')/(m+1)
                j=count
                for i in range(k-2,-1,-1):
                        if(word[i]==0):
                                B[i]=B[i+1]/(m+1)
                        elif(word[i]==1):
                                B[i]=S[j]/(m+1)
                                S[j]=S[j]+B[i+1]
                                j=j-1
                        L[i]=T*B[i]+L[i]
                L[k-1]=T*B[k-1]+L[k-1]
        return(L)
  def dual(a):
        b=list()
        b=a
        b=b[::-1]
        for i in range(len(b)):
                b[i]=1-b[i]
        return(b)
  def zeta(a):
        b=dual(a)
        l1=Li(a)+[1]
        l2=Li(b)+[1]
        Z=RIF('0')
        for i in range(len(l1)):
                Z=Z+l1[i]*l2[len(a)-i]
        return(Z)
  u=zeta(a)
  RIF=RealIntervalField(int(3.321928*D))
  u=u/1
  print(u)
}}}
{{attachment:akhi5.png}}
== Program to Compute Integer Relation between Multiple Zeta Values ==
{{{#!sagecell
from mpmath import *
print("Enter the number of composition")
@interact
def _( n=(5,(2..100))):
 a=[]
 for i in range(n):
        a.append([i+2,1])
 print("In each box Enter composition as an array")
 @interact
 def _(v=('Compositions', input_box( default=a, to_value=lambda x: vector(flatten(x)))), accuracy=(100..100000)):
  D=accuracy
  R=RealField(10)
  a=v
  def comptobin(a):
        word=[]
        for i in range(len(a)):
                word=word+[0]*(a[i]-1)+[1]
        return(word)
  DD=int(D)+int(R(log(3.321928*D))/R(log(10)))+4
  RIF=RealIntervalField(DD)
  mp.dps=DD
  def Li(word):
        n=int(DD*log(10)/log(2))+1
        B=[]
        L=[]
        S=[]
        count=-1
        k=len(word)
        for i in range(k):
                B.append(mpf('0'))
                L.append(mpf('0'))
                if(word[i]==1 and i<k-1):
                        S.append(mpf('0'))
                        count=count+1
        T=mpf('1')
        for m in range(n):
                T=T/2
                B[k-1]=mpf('1')/(m+1)
                j=count
                for i in range(k-2,-1,-1):
                        if(word[i]==0):
                                B[i]=B[i+1]/(m+1)
                        elif(word[i]==1):
                                B[i]=S[j]/(m+1)
                                S[j]=S[j]+B[i+1]
                                j=j-1
                        L[i]=T*B[i]+L[i]
                L[k-1]=T*B[k-1]+L[k-1]
        return(L)
  def dual(a):
        b=list()
        b=a
        b=b[::-1]
        for i in range(len(b)):
                b[i]=1-b[i]
        return(b)
  def zeta(a):
        b=dual(a)
        l1=Li(a)+[1]
        l2=Li(b)+[1]
        Z=mpf('0')
        for i in range(len(l1)):
                Z=Z+l1[i]*l2[len(a)-i]
        return(Z)
  zet=[]
  for i in range(n):
        zet.append((zeta(comptobin(a[i]))))
  mp.dps=D
  for i in range(n):
        zet[i]=zet[i]/1
        print("zeta(", a[i], ")=", zet[i])
  u=pslq(zet,tol=10**-D,maxcoeff=100,maxsteps=10000)
  print("the Intger Relation between the above zeta values given by the vector")
  print(u)
}}}
{{attachment:akhi10.png}}
== Word to composition ==
{{{#!sagecell
@interact
def _( weight=(7,(2..100))):
 n=weight
 a=[0 for i in range(n-1)]
 a.append(1)
 @interact
 def _(v=('word', input_grid(1, n, default=[a], to_value=lambda x: vector(flatten(x))))):
  a=[v[i] for i in range(len(v))]
  def bintocomp(a):
 b=[]
 count=1
 for j in range(len(a)):
  if(a[j]==0):
   count=count+1
  else:
   b.append(count)
   count=1
 return(b)
  print("Composition is {}".format(bintocomp(a)))
}}}

{{attachment:akhi2.png}}
== Composition to Word ==
{{{#!sagecell
@interact
def _( Depth=(7,(1..100))):
 n=Depth
 a=[]
 a.append(2)
 a=a+[1 for i in range(1,n)]
 @interact
 def _(v=('composition', input_grid(1, n, default=[a], to_value=lambda x: vector(flatten(x))))):
  a=[v[i] for i in range(len(v))]
  def comptobin(a):
 word=[]
 for i in range(len(a)):
  word=word+[0]*(a[i]-1)+[1]
 return(word)

  print("Word is {}".format(comptobin(a)))
}}}

{{attachment:akhi3.png}}
== Dual of a Word ==
{{{#!sagecell
@interact
def _( weight=(7,(2..100))):
 n=weight
 a=[0 for i in range(n-1)]
 a.append(1)
 @interact
 def _(v=('word', input_grid(1, n, default=[a], to_value=lambda x: vector(flatten(x))))):
  a=[v[i] for i in range(len(v))]
  def dual(a):
 b=list()
 b=a
 b=b[::-1]
 for i in range(len(b)):
  b[i]=1-b[i]
 return(b)

  print("Dual word is {}"?format(dual(a)))
}}}

{{attachment:akhi4.png}}


== Shuffle product of two Words ==
{{{#!sagecell
@interact
def _( w1=(2,(2..100)), w2=(2,(2..100))):
 a=[0]
 b=[0 for i in range(w2-1)]
 a=a+[1 for i in range(1,w1)]
 b=b+[1]
 import itertools
 #this program gives the list of all binary words of weight n and depth k
 @interact
 def _(v1=('word1', input_grid(1, w1, default=[a], to_value=lambda x: vector(flatten(x)))), v2=('word2', input_grid(1, w2, default=[b], to_value=lambda x: vector(flatten(x))))):
  a=[v1[i] for i in range(len(v1))]
  b=[v2[i] for i in range(len(v2))]
  def kbits(n, k):
    result = []
    for bits in itertools.combinations(range(n), k):
        s = ['0'] * n
        for bit in bits:
            s[bit] = '1'
        result.append(''.join(s))
    return result
  def sort(a,l,m):
        b=[]
        n=len(a)
        for i in range(n):
                b.append(a[i])
        for j in range(l-1,-1,-1):
                k=0
                for t in range(m+1):
                        for i in range(n):
                                if(a[i][j]== t):
                                        b[k]=a[i]
                                        k=k+1
                for i in range(n):
                        a[i]=b[i]
        return(a)
  def count(a):
        n=len(a)
        b=[]
        b.append(a[0])
        m=[]
        m.append(1)
        c=0
        for i in range(1,n):
                if(a[i]==a[i-1]):
                        m[c]=m[c]+1
                else:
                        b.append(a[i])
                        m.append(1)
                        c=c+1
        return(b,m)
  def shuffle(a,b):
        r=len(a)
        s=len(b)
        # Generating an array of strings containing all combinations of weight r+s and depth s
        M=kbits(r+s,s)
        n=len(M)
        a1= []
        for i in range(n):
                a1.append(list(M[i]))
        # The zeroes are replaced by the entries of a and the ones by the entries of b
        a2= []
        for i in range(n):
                a2.append([])
                count0=0
                count1=0
                for j in range(s+r):
                        if(a1[i][j]=='0'):
                                a2[i].append(a[count0])
                                count0=count0+1
                        if(a1[i][j]=='1'):
                                a2[i].append(b[count1])
                                count1=count1+1
        # Reordering in lexicographic order the entries of a2: this is done by first reordering them according to the last digit, then the next to last digit, etc
        a3=sort(a2,r+s,max(a+b+[0]))
        # Getting the same list without repetitions and with multiplicities
        a4=count(a3)
        return(a4)
  c=shuffle(a,b)
  for i in range(len(c[0])-1):
    print(c[1][i],"*",c[0][i] ,"+ ")
  print(c[1][len(c[0])-1],"*",c[0][len(c[0])-1])


}}}
{{attachment:akhi6.png}}
== Shuffle Regularization at 0 ==
{{{#!sagecell
@interact
def _( w=(2,(2..100))):
 a=[0]
 a=a+[1 for i in range(1,w)]
 import itertools
 #this program gives the list of all binary words of weight n and depth k
 @interact
 def _(v=('word', input_grid(1, w, default=[a], to_value=lambda x: vector(flatten(x))))):
  a=[v[i] for i in range(len(v))]
  def kbits(n, k):
    result = []
    for bits in itertools.combinations(range(n), k):
        s = ['0'] * n
        for bit in bits:
            s[bit] = '1'
        result.append(''.join(s))
    return result
  def sort(a,l,m):
 b=[]
 n=len(a)
 for i in range(n):
  b.append(a[i])
 for j in range(l-1,-1,-1):
  k=0
  for t in range(m+1):
   for i in range(n):
    if(a[i][j]== t):
     b[k]=a[i]
     k=k+1
  for i in range(n):
   a[i]=b[i]
 return(a)

  def sort1(a,l,m):
 b=[]
 b.append([])
 b.append([])
 n=len(a[0])
 for i in range(n):
  b[0].append(a[0][i])
  b[1].append(a[1][i])
 for j in range(l-1,-1,-1):
  k=0
  for t in range(m+1):
   for i in range(n):
    if(a[0][i][j]== t):
     b[0][k]=a[0][i]
     b[1][k]=a[1][i]
     k=k+1
  for i in range(n):
   a[0][i]=b[0][i]
   a[1][i]=b[1][i]
 return(a)

  def count(a):
 n=len(a)
 b=[]
 b.append(a[0])
 m=[]
 m.append(1)
 c=0
 for i in range(1,n):
  if(a[i]==a[i-1]):
   m[c]=m[c]+1
  else:
   b.append(a[i])
   m.append(1)
   c=c+1
 return(b,m)


  def count1(a):
 n=len(a[0])
 b=[]
 b.append([])
 b.append([])
 b[0].append(a[0][0])
 b[1].append(a[1][0])
 c=0
 for i in range(1,n):
  if(a[0][i]==a[0][i-1]):
   b[1][c]=b[1][c]+a[1][i]
  else:
   b[0].append(a[0][i])
   b[1].append(a[1][i])
   c=c+1

 return(b)
  def shuffle(a,b):
        r=len(a)
        s=len(b)
        # Generating an array of strings containing all combinations of weight r+s and depth s
        M=kbits(r+s,s)
        n=len(M)
        a1= []
        for i in range(n):
                a1.append(list(M[i]))
        # The zeroes are replaced by the entries of a and the ones by the entries of b
        a2= []
        for i in range(n):
                a2.append([])
                count0=0
                count1=0
                for j in range(s+r):
                        if(a1[i][j]=='0'):
                                a2[i].append(a[count0])
                                count0=count0+1
                        if(a1[i][j]=='1'):
                                a2[i].append(b[count1])
                                count1=count1+1
        # Reordering in lexicographic order the entries of a2: this is done by first reordering them according to the last digit, then the next to last digit, etc
        a3=sort(a2,r+s,max(a+b+[0]))
        # Getting the same list without repetitions and with multiplicities
        a4=count(a3)
        return(a4)
  def Regshuf0(a):
        r=[]
        r.append([])
        r.append([])
        t=0
        c=1
        for i in range(len(a)+1):
                if(t==0):
                        b=shuffle(a[:len(a)-i],a[len(a)-i:])
                        for j in range(len(b[0])):
                                r[0].append(b[0][j])
                                r[1].append(b[1][j]*c)
                        c=-c
                        if(i<len(a)):
                                if(a[len(a)-1-i]==1):
                                        t=1
        r=sort1(r,len(a),max(a+[0]))
        r=count1(r)
        rg=[]
        rg.append([])
        rg.append([])
        for i in range(len(r[0])):
                if(r[1][i] is not 0):
                        rg[0].append(r[0][i])
                        rg[1].append(r[1][i])
        return(rg)
  c = Regshuf0(a)
  for i in range(len(c[0])-1):
    if(c[1][i] != 0):
      print(c[1][i],"*",c[0][i] ,"+ ")
  if(c[1][len(c[0])-1] != 0):
    print(c[1][len(c[0])-1],"*",c[0][len(c[0])-1])


}}}
{{attachment:akhi7.png}}
== Shuffle Regularization at 1 ==
{{{#!sagecell
@interact
def _( w=(2,(2..20))):
 a=[0]
 a=a+[1 for i in range(1,w)]
 import itertools
 #this program gives the list of all binary words of weight n and depth k
 @interact
 def _(v=('word', input_grid(1, w, default=[a], to_value=lambda x: vector(flatten(x))))):
  a=[v[i] for i in range(len(v))]
  def kbits(n, k):
    result = []
    for bits in itertools.combinations(range(n), k):
        s = ['0'] * n
        for bit in bits:
            s[bit] = '1'
        result.append(''.join(s))
    return result
  def sort(a,l,m):
 b=[]
 n=len(a)
 for i in range(n):
  b.append(a[i])
 for j in range(l-1,-1,-1):
  k=0
  for t in range(m+1):
   for i in range(n):
    if(a[i][j]== t):
     b[k]=a[i]
     k=k+1
  for i in range(n):
   a[i]=b[i]
 return(a)

  def sort1(a,l,m):
 b=[]
 b.append([])
 b.append([])
 n=len(a[0])
 for i in range(n):
  b[0].append(a[0][i])
  b[1].append(a[1][i])
 for j in range(l-1,-1,-1):
  k=0
  for t in range(m+1):
   for i in range(n):
    if(a[0][i][j]== t):
     b[0][k]=a[0][i]
     b[1][k]=a[1][i]
     k=k+1
  for i in range(n):
   a[0][i]=b[0][i]
   a[1][i]=b[1][i]
 return(a)

  def count(a):
 n=len(a)
 b=[]
 b.append(a[0])
 m=[]
 m.append(1)
 c=0
 for i in range(1,n):
  if(a[i]==a[i-1]):
   m[c]=m[c]+1
  else:
   b.append(a[i])
   m.append(1)
   c=c+1
 return(b,m)


  def count1(a):
 n=len(a[0])
 b=[]
 b.append([])
 b.append([])
 b[0].append(a[0][0])
 b[1].append(a[1][0])
 c=0
 for i in range(1,n):
  if(a[0][i]==a[0][i-1]):
   b[1][c]=b[1][c]+a[1][i]
  else:
   b[0].append(a[0][i])
   b[1].append(a[1][i])
   c=c+1

 return(b)
  def shuffle(a,b):
        r=len(a)
        s=len(b)
        # Generating an array of strings containing all combinations of weight r+s and depth s
        M=kbits(r+s,s)
        n=len(M)
        a1= []
        for i in range(n):
                a1.append(list(M[i]))
        # The zeroes are replaced by the entries of a and the ones by the entries of b
        a2= []
        for i in range(n):
                a2.append([])
                count0=0
                count1=0
                for j in range(s+r):
                        if(a1[i][j]=='0'):
                                a2[i].append(a[count0])
                                count0=count0+1
                        if(a1[i][j]=='1'):
                                a2[i].append(b[count1])
                                count1=count1+1
        # Reordering in lexicographic order the entries of a2: this is done by first reordering them according to the last digit, then the next to last digit, etc
        a3=sort(a2,r+s,max(a+b+[0]))
        # Getting the same list without repetitions and with multiplicities
        a4=count(a3)
        return(a4)
  def Regshuf1(a):
 r=[]
 r.append([])
 r.append([])
 t=0
 c=1
 for i in range(len(a)+1):
  if(t==0):
   b=shuffle(a[:i],a[i:])
   for j in range(len(b[0])):
    r[0].append(b[0][j])
    r[1].append(b[1][j]*c)
   c=-c
   if(i<len(a)):
    if(a[i]==0):
     t=1
 r=sort1(r,len(a),max(a+[0]))
 r=count1(r)
 rg=[]
 rg.append([])
 rg.append([])
 for i in range(len(r[0])):
  if(r[1][i] is not 0):
   rg[0].append(r[0][i])
   rg[1].append(r[1][i])
 return(rg)
  c = Regshuf1(a)
  for i in range(len(c[0])-1):
    if(c[1][i] != 0):
      print(c[1][i],"*",c[0][i] ,"+ ")
  if(c[1][len(c[0])-1] != 0):
    print(c[1][len(c[0])-1],"*",c[0][len(c[0])-1])


}}}
{{attachment:akhi8.png}}

Integer Factorization

Divisibility Poset

by William Stein

divposet.png

Factor Trees

by William Stein

factortree.png

More complicated demonstration using Mathematica: http://demonstrations.wolfram.com/FactorTrees/

Factoring an Integer

by Timothy Clemans

Sage implementation of the Mathematica demonstration of the same name. http://demonstrations.wolfram.com/FactoringAnInteger/

Prime Numbers

Illustrating the prime number theorem

by William Stein

primes.png

Prime Spiral - Square FIXME

by David Runde

SquareSpiral.PNG

Prime Spiral - Polar

by David Runde

Needs fix for show_factors

PolarSpiral.PNG

Modular Forms

Computing modular forms

by William Stein

modformbasis.png

Computing the cuspidal subgroup

by William Stein

ncols not working

cuspgroup.png

A Charpoly and Hecke Operator Graph

by William Stein

heckegraph.png

Modular Arithmetic

Quadratic Residue Table FIXME

by Emily Kirkman

quadres.png

quadresbig.png

Cubic Residue Table FIXME

by Emily Kirkman

cubres.png

Cyclotomic Fields

Gauss and Jacobi Sums in Complex Plane

by Emily Kirkman

jacobising.png

Exhaustive Jacobi Plotter

by Emily Kirkman

jacobiexh.png

Elliptic Curves

Adding points on an elliptic curve

by David Møller Hansen

PointAddEllipticCurve.png

Plotting an elliptic curve over a finite field

ellffplot.png

Cryptography

The Diffie-Hellman Key Exchange Protocol

by Timothy Clemans and William Stein

dh.png

Other

Continued Fraction Plotter

by William Stein

crows not working

contfracplot.png

Computing Generalized Bernoulli Numbers

by William Stein (Sage-2.10.3)

bernoulli.png

Fundamental Domains of SL_2(ZZ)

by Robert Miller

fund_domain.png

Multiple Zeta Values

by Akhilesh P.

Computing Multiple Zeta values

Word Input

akhi1.png

Composition Input

akhi5.png

Program to Compute Integer Relation between Multiple Zeta Values

akhi10.png

Word to composition

akhi2.png

Composition to Word

akhi3.png

Dual of a Word

akhi4.png

Shuffle product of two Words

akhi6.png

Shuffle Regularization at 0

akhi7.png

Shuffle Regularization at 1

akhi8.png

interact/number_theory (last edited 2020-06-14 09:10:48 by chapoton)