<> = 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}} == Factor Trees == by William Stein {{{#!sagecell import random def ftree(rows, v, i, F): if len(v) > 0: # add a row to g at the ith level. rows.append(v) w = [] for i in range(len(v)): k, _, _ = v[i] if k is None or is_prime(k): w.append((None,None,None)) else: d = random.choice(divisors(k)[1:-1]) w.append((d,k,i)) e = k//d if e == 1: w.append((None,None)) else: w.append((e,k,i)) if len(w) > len(v): ftree(rows, w, i+1, F) def draw_ftree(rows,font): g = Graphics() for i in range(len(rows)): cur = rows[i] for j in range(len(cur)): e, f, k = cur[j] if not e is None: if is_prime(e): c = (1,0,0) else: c = (0,0,.4) g += text(str(e), (j*2-len(cur),-i), fontsize=font, rgbcolor=c) if not k is None and not f is None: g += line([(j*2-len(cur),-i), ((k*2)-len(rows[i-1]),-i+1)], alpha=0.5) return g @interact def factor_tree(n=100, font=(10, (8..20)), redraw=['Redraw']): n = Integer(n) rows = [] v = [(n,None,0)] ftree(rows, v, 0, factor(n)) show(draw_ftree(rows, font), axes=False) }}} {{attachment: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/ {{{#!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 = == Illustrating the prime number theorem == by William Stein {{{#!sagecell @interact def _(N=(100,list(range(2,2000)))): pretty_print(html(r"$\pi(x)$ and $x/(\log(x)-1)$ for $x < %s$"%N)) show(plot(prime_pi, 0, N, color='red') + plot(x/(log(x)-1), 5, N, color='blue')) }}} {{attachment:primes.png}} == Prime Spiral - Square FIXME == by David Runde {{{#!sagecell @interact def square_prime_spiral(start=1, end=100, size_limit = 10, show_lines=false, invert=false, x_cord=0, y_cord=0, n = 0): """ REFERENCES: Alpern, Dario. "Ulam's Spiral". http://www.alpertron.com.ar/ULAM.HTM Sacks, Robert. http://www.NumberSpiral.com Ventrella, Jeffery. "Prime Numbers are the Holes Behind Complex Composite Patterns". http://www.divisorplot.com Williamson, John. Number Spirals. http://www.dcs.gla.ac.uk/~jhw/spirals/index.html jhw@dcs.gla.ac.uk Weisstein, Eric W. "Prime-Generating Polynomial." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/Prime-GeneratingPolynomial.html """ #Takes an (x,y) coordinate (and the start of the spiral) and gives its corresponding n value def find_n(x,y, start): if x>0 and y>-x and y<=x: return 4*(x-1)^2 + 5*(x-1) + (start+1) + y elif x<=0 and y>=x and y<=-x: return 4*x^2 - x + (start) -y elif y>=0 and -y+1 <= x and y-1 >= x: return 4*y^2 -y + start -x elif y<0 and -x >= y and y end: print("WARNING: n is larger than the end value") #Changes the entry of a matrix by taking the old matrix and the (x,y) coordinate (in matrix coordinates) and returns the changed matrix def matrix_morph(M, x, y, set): N = copy(M) N[x-1,y] = set M = N return M #These functions return an int based on where the t is located in the spiral def SW_NE(t, x, y, start): if -y n =', find_n(x_cord, y_cord, start)) print(' ') print("SW/NE line") if -y_cord pixel on) #matrix_plot(M) }}} {{attachment:SquareSpiral.PNG}} == Prime Spiral - Polar == by David Runde 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") def f(n): return (sqrt(n)*cos(2*pi*sqrt(n)), sqrt(n)*sin(2*pi*sqrt(n))) list =[] list2=[] if show_factors == false: for i in [start..end]: if i.is_pseudoprime(): list.append(f(i-start+1)) #Primes list else: list2.append(f(i-start+1)) #Composites list P = points(list) R = points(list2, alpha = .1) #Faded Composites else: for i in [start..end]: list.append(disk((f(i-start+1)),0.05*pow(2,len(factor(i))-1), (0,2*pi))) #resizes each of the dots depending of the number of factors of each number if i.is_pseudoprime() and highlight_primes: list2.append(f(i-start+1)) P = plot(list) p_size = 5 #the orange dot size of the prime markers if not highlight_primes: list2 = [(f(n-start+1))] R=points(list2, hue = .1, pointsize = p_size) if n > 0: print('n = {}'.format(factor(n))) p = 1 #The X which marks the given n W1 = disk((f(n-start+1)), p, (pi/6, 2*pi/6)) W2 = disk((f(n-start+1)), p, (4*pi/6, 5*pi/6)) W3 = disk((f(n-start+1)), p, (7*pi/6, 8*pi/6)) W4 = disk((f(n-start+1)), p, (10*pi/6, 11*pi/6)) Q = plot(W1+W2+W3+W4, alpha = .1) n=n-start+1 #offsets the n for different start values to ensure accurate plotting if show_curves: begin_curve = 0 t = var('t') a=1 b=0 if n > (floor(sqrt(n)))^2 and n <= (floor(sqrt(n)))^2 + floor(sqrt(n)): c = -((floor(sqrt(n)))^2 - n) c2= -((floor(sqrt(n)))^2 + floor(sqrt(n)) - n) else: c = -((ceil(sqrt(n)))^2 - n) c2= -((floor(sqrt(n)))^2 + floor(sqrt(n)) - n) print('Pink Curve: n^2 +', c) print('Green Curve: n^2 + n +', c2) def g(m): return (a*m^2+b*m+c); def r(m) : return sqrt(g(m)) def theta(m) : return r(m)- m*sqrt(a) S1 = parametric_plot(((r(t))*cos(2*pi*(theta(t))),(r(t))*sin(2*pi*(theta(t)))), begin_curve, ceil(sqrt(end-start)), rgbcolor=hue(0.8), thickness = .2) #Pink Line b=1 c= c2; S2 = parametric_plot(((r(t))*cos(2*pi*(theta(t))),(r(t))*sin(2*pi*(theta(t)))), begin_curve, ceil(sqrt(end-start)), rgbcolor=hue(0.6), thickness = .2) #Green Line show(R+P+S1+S2+Q, aspect_ratio = 1, axes = false) else: show(R+P+Q, aspect_ratio = 1, axes = false) else: show(R+P, aspect_ratio = 1, axes = false) }}} {{attachment:PolarSpiral.PNG}} = 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('

Cuspidal Subgroups of Modular Jacobians J0(N)

')) @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("

Charpoly and Hecke Graph: Level %s, T_%s

"%(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 == by Emily Kirkman {{{#!sagecell from numpy import array as narray @interact def quad_res_plot(first_n_odd_primes = (20,200),display_size=[7..15]): # Compute numpy matrix of legendre symbols r = int(first_n_odd_primes) np = [nth_prime(i+2) for i in range(r)] leg = [[legendre_symbol(np[i], np[j]) for i in range(r)] for j in range(r)] na = narray(leg) for i in range(r): for j in range(r): if na[i][j] == 1 and Mod((np[i]-1)*(np[j]-1)//4,2) == 0: na[i][j] = 2 m = matrix(na) # Define plot structure MP = matrix_plot(m, cmap='Oranges') for i in range(r): if np[-1] < 100: MP += text('%d'%nth_prime(i+2),(-.75,r-i-.5), rgbcolor='black') MP += text('%d'%nth_prime(i+2), (i+.5,r+.5), rgbcolor='black') if len(np) < 75: MP += line([(0,i),(r,i)], rgbcolor='black') MP += line([(i,0),(i,r)], rgbcolor='black') if np[-1] < 100: for i in range(r): # rows for j in range(r): # cols if m[j][i] == 0: MP += text('0',(i+.5,r-j-.5),rgbcolor='black') elif m[j][i] == -1: MP += text('N',(i+.5,r-j-.5),rgbcolor='black') elif m[j][i] == 1: MP += text('A',(i+.5,r-j-.5),rgbcolor='black') elif m[j][i] == 2: MP += text('S',(i+.5,r-j-.5),rgbcolor='black') MP += line([(0,r),(r,r)], rgbcolor='black') MP += line([(r,0),(r,r)], rgbcolor='black') MP += line([(0,0),(r,0)], rgbcolor='black') MP += line([(0,0),(0,r)], rgbcolor='black') if len(np) < 75: MP += text('q',(r/2,r+2), rgbcolor='black', fontsize=15) MP += text('p',(-2.5,r/2), rgbcolor='black', fontsize=15) MP.show(axes=False, ymax=r, figsize=[display_size,display_size]) html('Symmetry of Prime Quadratic Residues mod the first %d odd primes.'%r) }}} {{attachment:quadres.png}} {{attachment:quadresbig.png}} == Cubic Residue Table FIXME == by Emily Kirkman {{{#!sagecell def power_residue_symbol(alpha, p, m): if p.divides(alpha): return 0 if not p.is_prime(): return prod(power_residue_symbol(alpha,ell,m)^e for ell, e in p.factor()) F = p.residue_field() N = p.norm() r = F(alpha)^((N-1)/m) k = p.number_field() for kr in k.roots_of_unity(): if r == F(kr): return kr def cubic_is_primary(n): g = n.gens_reduced()[0] a,b = g.polynomial().coefficients() return Mod(a,3)!=0 and Mod(b,3)==0 from numpy import array as narray @interact def cubic_sym(n=(10..35),display_size=[7..15]): # Compute numpy matrix of primary cubic residue symbols r = n m=3 D. = NumberField(x^2+x+1) it = D.primes_of_degree_one_iter() pp = [] while len(pp) < r: k = it.next() if cubic_is_primary(k): pp.append(k) n = narray([ [ power_residue_symbol(pp[i].gens_reduced()[0], pp[j], m) \ for i in range(r) ] for j in range(r) ]) # Convert to integer matrix for gradient colors for i in range(r): for j in range(r): if n[i][j] == w: n[i][j] = int(-1) elif n[i][j] == w^2: n[i][j] = int(-2) elif n[i][j] == 1: n[i][j] = int(1) m = matrix(n) # Define plot structure MP = matrix_plot(m,cmap="Blues") for i in range(r): MP += line([(0,i),(r,i)], rgbcolor='black') MP += line([(i,0),(i,r)], rgbcolor='black') for j in range(r): if m[i][j] == -2: MP += text(r'$\omega^2$',(i+.5,r-j-.5),rgbcolor='black') if m[i][j] == -1: MP += text(r'$\omega $',(i+.5,r-j-.5),rgbcolor='black') if m[i][j] == 0: MP += text('0',(i+.5,r-j-.5),rgbcolor='black') if m[i][j] == 1: MP += text('R',(i+.5,r-j-.5),rgbcolor='white') MP += line([(0,r),(r,r)], rgbcolor='black') MP += line([(r,0),(r,r)], rgbcolor='black') MP += line([(0,0),(r,0)], rgbcolor='black') MP += line([(0,0),(0,r)], rgbcolor='black') 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)) MP.show(axes=False, figsize=[display_size,display_size]) }}} {{attachment:cubres.png}} = Cyclotomic Fields = == Gauss and Jacobi Sums in Complex Plane == by Emily Kirkman {{{#!sagecell def jacobi_sum(e,f): # If they are both trivial, return p if e.is_trivial() and f.is_trivial(): return (e.parent()).order() + 1 # If they are inverses of each other, return -e(-1) g = e*f if g.is_trivial(): return -e(-1) # If both are nontrivial, apply mult. formula: elif not e.is_trivial() and not f.is_trivial(): return e.gauss_sum()*f.gauss_sum()/g.gauss_sum() # If exactly one is trivial, return 0 else: return 0 def latex2(e): return latex(list(e.values_on_gens())) def jacobi_plot(p, e_index, f_index, with_text=True): # Set values G = DirichletGroup(p) e = G[e_index] f = G[f_index] ef = e*f js = jacobi_sum(e,f) e_gs = e.gauss_sum() f_gs = f.gauss_sum() ef_gs = (e*f).gauss_sum() # Compute complex coordinates f_pt = list(f.values_on_gens()[0].complex_embedding()) e_pt = list(e.values_on_gens()[0].complex_embedding()) ef_pt = list(ef.values_on_gens()[0].complex_embedding()) f_gs_pt = list(f_gs.complex_embedding()) e_gs_pt = list(e_gs.complex_embedding()) ef_gs_pt = list(ef_gs.complex_embedding()) try: js = int(js) js_pt = [CC(js)] except: js_pt = list(js.complex_embedding()) # Define plot structure 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') if with_text: S += text('$J(%s,%s) = %s$'%(latex2(e),latex2(f),latex(js)), (3,2.5),fontsize=15, rgbcolor='black') else: html('$$J(%s,%s) = %s$$'%(latex2(e),latex2(f),latex(js))) return S @interact def single_jacobi_plot(p=prime_range(3,100), e_range=(0..100), f_range=(0..100)): e_index = floor((p-2)*e_range/100) f_index = floor((p-2)*f_range/100) S = jacobi_plot(p,e_index,f_index,with_text=False) S.show(aspect_ratio=1) }}} {{attachment:jacobising.png}} == Exhaustive Jacobi Plotter == by Emily Kirkman {{{#!sagecell def jacobi_sum(e,f): # If they are both trivial, return p if e.is_trivial() and f.is_trivial(): return (e.parent()).order() + 1 # If they are inverses of each other, return -e(-1) g = e*f if g.is_trivial(): return -e(-1) # If both are nontrivial, apply mult. formula: elif not e.is_trivial() and not f.is_trivial(): return e.gauss_sum()*f.gauss_sum()/g.gauss_sum() # If exactly one is trivial, return 0 else: return 0 def latex2(e): return latex(list(e.values_on_gens())) def jacobi_plot(p, e_index, f_index, with_text=True): # Set values G = DirichletGroup(p) e = G[e_index] f = G[f_index] ef = e*f js = jacobi_sum(e,f) e_gs = e.gauss_sum() f_gs = f.gauss_sum() ef_gs = (e*f).gauss_sum() # Compute complex coordinates f_pt = list(f.values_on_gens()[0].complex_embedding()) e_pt = list(e.values_on_gens()[0].complex_embedding()) ef_pt = list(ef.values_on_gens()[0].complex_embedding()) f_gs_pt = list(f_gs.complex_embedding()) e_gs_pt = list(e_gs.complex_embedding()) ef_gs_pt = list(ef_gs.complex_embedding()) try: js = int(js) js_pt = [CC(js)] except: js_pt = list(js.complex_embedding()) # Define plot structure 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') if with_text: S += text('$J(%s,%s) = %s$'%(latex2(e),latex2(f),latex(js)), (3,2.5),fontsize=15, rgbcolor='black') else: pretty_print(html('$$J(%s,%s) = %s$$'%(latex2(e),latex2(f),latex(js)))) return S @interact def exhaustive_jacobi_plot(p=prime_range(3,8)): ga = [jacobi_plot(p,i,j) for i in range(p-1) for j in range(p-1)[i:]] for i in range(len(ga)): ga[i].save('j%d.png'%i,figsize=4,aspect_ratio=1, xmin=-2.5,xmax=5, ymin=-2.5, ymax=2.5) # Since p is odd, will have n = p-1 even. So 1+2+...+n = (n/2)*(n+1). # We divide this by rows of 2. rows = ceil(p*(p-1)/4) s='' for i in range(rows): s+=''%(2*i) s+=''%(2*i+1) s+='
' pretty_print(html(s)) }}} {{attachment:jacobiexh.png}} = 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("""

%s-Bit Diffie-Hellman Key Exchange

  1. Alice and Bob agree to use the prime number p = %s and base g = %s.
  2. Alice chooses the secret integer a = %s, then sends Bob (ga mod p):
    %s%s mod %s = %s.
  3. Bob chooses the secret integer b=%s, then sends Alice (gb mod p):
    %s%s mod %s = %s.
  4. Alice computes (gb mod p)a mod p:
    %s%s mod %s = %s.
  5. Bob computes (ga mod p)b mod p:
    %s%s mod %s = %s.
""" % (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"

First n=%s Bernoulli numbers attached to characters with modulus m=%s

"%(n,m) s += r'' s += r'' + \ ''.join(r''%k for k in [1..n]) + '' for eps in G.list(): v = ''.join([''%latex(eps.bernoulli(k)) for k in [1..n]]) s += '%s\n'%( eps, eps.conductor(), v) s += '
$\chi$Conductor$B_{%s,\chi}$
$%s$
%s%s
' 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