Chapter Ten:

Additional Examples

10.1 Animating the Derivative

Clear[h, t]

Do[Plot[{(Sin[t+h] - Sin[t])/h, Cos[t]}, {t, -5, 5},
	PlotRange->{{-5, 5}, {-1, 1}}, Ticks->{Automatic, {-1, 1}},
	PlotStyle->{{}, Dashing[{.02, .02}]}],    {h, 1, .1, -.1}]

DerivativeApproximationMovie[
        f_, xmin_, xmax_, ymin_, ymax_, h0_:1, h1_:.1, n_:9] :=
Block[{ff, ffder},
	ff[t_] :=   f /. x -> t;
	ffder[t_] := D[f, x] /. x -> t;
	Do[Plot[{(ff[t+h] - ff[t])/h, ffder[t]}, {t, xmin, xmax},
		PlotRange->{{xmin, xmax}, {ymin, ymax}},
		PlotStyle->{{}, Dashing[{.02}]}],
	{h, h0, h1, (h1 - h0)/(n - 1)}]]

DerivativeApproximationMovie[x^10, 0, 1, 0, 30, 1, .02, 10]

f[x_] := Sin[x];
xmin = -2; xmax = 2; delta = (xmax - xmin)/10;
plotf = Plot[f[x], {x, xmin, xmax}, DisplayFunction->Identity];
plotder[t_] := Plot[f'[x], {x, xmin - delta/100, t},
								DisplayFunction->Identity];
yrange = PlotRange[Plot[{f[x], f'[x]}, {x, xmin, xmax},
		PlotRange->All, DisplayFunction->Identity]][[2]];

Do[Show[plotf, plotder[t],
		Graphics[{Line[{{t, f[t]} - delta/3 {1, f'[t]},
                        {t, f[t]} + delta/3 {1, f'[t]}}],
  		Dashing[{.008}], Thickness[.003],
  		Line[{{t, f[t]}, {t, f'[t]}}]}],
  	PlotRange->{{xmin - delta/3, xmax + delta/3}, yrange},
  	DisplayFunction->$DisplayFunction],
{t, xmin, xmax, delta}]

10.2 Billiard Paths on Elliptical Tables

Clear[f]; 
{a, b}  = {2, 1};
pi = N[Pi];
f[t_]  := {a Cos[t], b Sin[t]}; Attributes[f] = Listable;
ellipse = ParametricPlot[f[t], {t, 0, 2 Pi}];
foci    = {{Sqrt[a^2-b^2], 0}, {-Sqrt[a^2-b^2], 0}};

slope1[t0_, t1_] := Divide @@ Reverse[f[t1] - f[t0]] /;
									   	    t0+t1 != 2 pi;
slope2[t1_] := -b Cot[t1] / a           /; t1 != 0 && t1 != pi;

slope3[t0_, t1_] :=
	Tan[2 ArcTan[slope2[t1]] - ArcTan[slope1[t0,t1]]] /; 
 					   t1 != 0 && t1 != pi && t0+t1 != 2 pi;
slope3[t0_, t1_] := -slope1[t0, t1]     /; t1 == 0 || t1 == pi;
slope3[t0_, t1_] := Tan[2 ArcTan[slope2[t1]] - pi/2] /;
										    t0+t1 == 2 pi;

newpoint[{t0_, t1_}] := {t1,
 First[Select[ArcTan[b x, a y] /. 
  N[Solve[(x/a)^2  + (y/b)^2  == 1 && 
	      y - f[t1][[2]] == (x - f[t1][[1]]) slope3[t0, t1]]], 
	  Abs[# - t1] > 10^-15 & ]]}

BilliardPath[start_, n_] := Show[ellipse, Graphics[{
	Thickness[.0001], PointSize[.02],
	Point /@ f[start], Point /@ foci,
     Line[f[First /@ NestList[newpoint, N[start], n]]]}]]

BilliardPath[{Pi - .2, Pi/2}, 35];

BilliardPath[{Pi + .5, Pi/2 - .6}, 35];

BilliardPath[{Pi+ArcCos[Sqrt[3]/2], ArcCos[-Sqrt[3]/2]}, 4];

BilliardPath[{Pi, Pi/2}, 8];

10.3 The Art Gallery Theorem

SignedArea[{{x1_,y1_}, {x2_,y2_}, {x3_,y3_}}] :=
	Chop[Det[{{x1,y1,1}, {x2,y2,1}, {x3,y3,1}}]/2.]

Orientation[{{x1_,y1_}, {x2_,y2_}, {x3_,y3_}}] :=
	Sign[Chop[Det[{{x1,y1,1}, {x2,y2,1}, {x3,y3,1}}]]]

Leftmost[polygon_] :=
		First@First@Position[polygon, First[Sort[polygon]]]
   		
Orientation[polygon_] := Block[{n = Length[polygon]}, 	Orientation[polygon[[Leftmost[polygon] + {-1,0,1} /.
		{0 -> n, n+1 -> 1}]]]] 		/; Length[polygon] > 3

SignedArea[{{0,0}, {0,1}, {1,0}}]

Orientation[{{0,0}, {5,0}, {0,5}}]

LeftOf[p_, q_, r_] := Orientation[{p, q, r}] == 1

norm[p_] := Sqrt[p.p]

Between[p_, q_, r_] := Orientation[{p, q, r}] == 0 &&
	       (r-p) . (q-p) >= 0 && norm[r-p] <= norm[q-p]

Intersect[{a_, b_}, {c_, d_}] :=  
	 (LeftOf[a,b,c] ~Xor~ LeftOf[a,b,d]) &&
	 (LeftOf[c,d,a] ~Xor~ LeftOf[c,d,b])     ||
	 
	 Between[a,b,c] || Between[a,b,d]  ||
	 Between[c,d,a] || Between[c,d,b]

randompoint := Random[]*
		{Cos[angle = Random[Real, 2 Pi//N]], Sin[angle]};

good[point_] := And @@ Table[
	!Intersect[walk[[{i, i+1}]], {Last@walk, point}],
		{i, Length@walk - 2}]
		
goodlaststep[walk_] := And @@ Table[
	!Intersect[walk[[{i, i+1}]], {Last@walk, First@walk}],
		{i, 2, Length@walk - 2}]

RandomPolygon[n_, maxpts_:40] := (
	walk = {{0., 0.}, randompoint}; j = 2;
	Do[While[
		 j++ < maxpts && !good[new = Last[walk]+randompoint]];
  	   AppendTo[walk, new], {n - 3}];
	While[j++ < maxpts &&
     (!good[new = Last[walk]+randompoint] ||
      !goodlaststep[Append[walk, new]])];
	If[j > maxpts, 0, Join[walk, {new, First@walk}]])

ShowRandomPolygon[n_, maxpts_:40] := 
	If[!NumberQ[temp = RandomPolygon[n, maxpts]],
	Show[Graphics[Line[temp]], PlotRange->All]]

ShowRandomPolygon[12, 50]

MakeCounterclockwise[polygon_] := 
	If[Orientation[polygon] == -1, Reverse[polygon], polygon]

(* ConvexVertex[polygon_, i_] :=
				Orientation[polygon[[{i-1, i, i+1}]]] == 1 
				
Preceding, from book, would have been better as follows. Book
version doesn't work on first and last vertices *)

ConvexVertex[polygon_, i_] := Block[{n = Length @ polygon},
	Orientation[polygon[[
	(i + Range[-1,1]) /. {0 -> n, n+1 -> 1}]]] == 1]


ConvexVertex[polygon_] := Block[{i = 2},
				While[!ConvexVertex[polygon, i], i++];  i]

Inside[{p_, q_, r_}, point_] := 
	!LeftOf[q, p, point] && !LeftOf[r, q, point] &&
	!LeftOf[p, r, point] 

PolygonDiagonal[polygon_] :=
  Block[{i = ConvexVertex[polygon], invertices},
  invertices = Select[ Drop[polygon, {i-1, i+1}],
					 Inside[polygon[[{i-1, i, i+1}]], #] &];
  If[invertices == {}, {i-1, i+1},
    Prepend[
     Flatten@Position[polygon, Last@Sort[invertices, 
	 SignedArea[{polygon[[i-1]], #1, polygon[[i+1]]}] <=
	 SignedArea[{polygon[[i-1]], #2, polygon[[i+1]]}] &]],
      i]]]

Triangulate[polygon_] := {}  /; Length[polygon] == 3

Triangulate[polygon_] := Block[{diag, poly1, poly2},
	diag = PolygonDiagonal[polygon];
	poly1 = Take[polygon, diag];  
	poly2 = Drop[polygon, {1+diag[[1]], diag[[2]]-1}];
	Join[Triangulate[poly1], Triangulate[poly2],
 				{polygon[[diag]]}]]

ShowTriangulation[polygon_]:= 
	Block[{ccwpolygon = MakeCounterclockwise[polygon]},
	Show[Graphics[
	{Line[Append[polygon, polygon[[1]]]],
	 Map[Line, Triangulate[ccwpolygon]]}],
						    AspectRatio->1, PlotRange->All]]

test = {{1.2, 1}, {0, 0}, {1.2, 0}, {2, .1}, {1, .2}, {2, .3},
      {1, .4}, {2, .5}, {1.1, .55}, {2, .7}, {1, .7}, {2, .9}};

ConvexVertex[test]

PolygonDiagonal[test]

ShowTriangulation[test];

Show[Graphics[Line[{{0,0}, {0,1}, {1, 1}, {1, 0}, {0,0}}],
			AspectRatio->1]];

posn[polygon_, diag_] :=
	         Flatten[Map[Position[polygon, #] &, diag]]

convert[polygon_, diags_]:=
	{Range[Length@polygon], Map[posn[polygon, #] &, diags]} 
				 
Color[verts_, diags_] := {verts, {1,2,3}} /; Length[verts] == 3

Color[verts_, diags_]:= Block[{i = 1, oldcolors},

     While[!FreeQ[diags, verts[[++i]]]];
	
	(* preceding finds the ear i-1, i, i+1, since that is
	   where a vertex is not on any diagonal *)

	oldcolors = Color[
		Complement[verts, {verts[[i]]}],
		Complement[diags, {{verts[[i-1]], verts[[i+1]]}}]];

	(* preceding deletes vertex i and colors by recursion;
	   following colors vertex i appropriately *)

	{Append[oldcolors[[1]], verts[[i]]],
	Append[oldcolors[[2]], 6 -
	  oldcolors[[2,
 	    First@First@Position[oldcolors[[1]], verts[[i-1]]]]] -
 	  oldcolors[[2,
 	   First@First@Position[oldcolors[[1]], verts[[i+1]]]]]]}]

ThreeColor[polygon_] :=
    Block[{ccwpoly = MakeCounterclockwise[polygon]},
    {Color @@ convert[ccwpoly, Triangulate[ccwpoly]], ccwpoly}]
	
Guards[polygon_] := Block[{coloring, ccw, colorclass, i},
  {coloring, ccw} = ThreeColor[polygon];
  colorclass = First@Sort[Table[
  (First /@ Select[Transpose[coloring], #[[2]]==i &]), {i, 3}],
                        Length[#1] < Length[#2]&];
 Show[Graphics[{Line[Append[polygon, First[polygon]]],
 		       PointSize[.04], Point /@ ccw[[colorclass]]}],
						    PlotRange->All, AspectRatio->1]]

Guards[test]

GuardsAndColors, a routine to show all 3 colors

GuardsAndColors[polygon_]:=Block[{coloring, ccw, i, colorclasses},
 {coloring, ccw} = ThreeColor[polygon];

 colorclasses = Sort[
Table[(First /@ Select[Transpose[coloring], #[[2]]==i &]), {i, 3}],
                  Length[#1] < Length[#2]&];
 Show[Graphics[Line[Append[polygon,polygon[[1]]]]],
	Graphics@Prepend[Point/@ccw[[colorclasses[[1]]]],
									PointSize@.06],
	Graphics@Prepend[Point/@ccw[[colorclasses[[2]]]],
									PointSize@.04],
	Graphics[{GrayLevel@.4,
			Prepend[Point/@ccw[[colorclasses[[3]]]],
			PointSize@.04]}],
						PlotRange->All,AspectRatio->1]]

GuardsAndColors[test]

10.4 Rational Enumeration

(* 2.0 note: Digits has been replaced by IntegerDigits 
   and RealDigits *)

Index[q_Rational] := (
t = Join@@{Digits@Numerator[q], {10}, Digits@Denominator[q]}) .
    (11^Range[Length[t] - 1, 0, -1])  		    /; Positive[q]
  
Index[n_Integer] := (t =
Join@@{Digits@n, {10,1}}) . (11^Range[Length[t]-1, 0, -1]) /;
 											 !Negative[n]

Off[Power::infy];

InverseIndex[n_] := Block[{temp, t}, 
 	temp = Take[t = Characters[ToString@BaseForm[n,11]],
	First @ First @ Position[t, "\n"] - 1] /. "a"-> "/";
   If[First[temp] == "/", Return["Failure"]];
   temp = ToExpression[StringJoin@@temp];
   If[!NumberQ[temp] || Index[temp] != n, "Failure", temp]] /; 											 !Negative[n]

(* preceding code changed from first printing of book so that
it works in both versions 1.2 and 2.0 *) 

InverseIndex[111] = 0;

SetAttributes[{Index, InverseIndex}, Listable]

Index[345/1001]

InverseIndex[66500754]

Index[Range[0,20]/20]

InverseIndex[%]

Select[InverseIndex[Range[1000]], NumberQ]

10.5 Algebraic Numbers

x /. Solve[x^3 - 9 x^2 + 27 x - 31 == 0]

Clear[n, t]

version 2.0 notes

Because 2.0 treats powers differently, some of this code may be subject to shortening. But it does work.

polynomial::usage = "polynomial[alg] returns a polynomial in  x with integer coefficients for which the algebraic integer alg is a root.";

polynomial[n_Integer] := x - n

polynomial[z_Complex] := x^2 - 2 Re[z] x + Abs[z]^2

polynomial[e_ + f_] := polysum[{polynomial[e], polynomial[f]}]

polynomial[e_ ^ x_ * f_ ^ x_] := polynomial[(e f)^x]

polynomial[e_ * f_] := polyprod[{polynomial[e], polynomial[f]}]

polysum[{p_, q_}] := Resultant[q /. x -> t, p /. x -> x-t, t]

polysum[list_] := polysum[{First@list, polysum[Rest@list]}] /;
									  Length[list] > 2

degree[p_] := Length[CoefficientList[p, x]] - 1

polyprod[{p_, q_}] := 
 Expand[(p /. x -> 0) ^ degree[q] *
	 Resultant[t^degree[p] (p /. x -> x/t), q /. x -> t, t]]

polyprod[list_] := 

polyprod[{First@list, polyprod[Rest@list]}] /; Length[list] > 2

polynomial[e_ ^ n_Integer?Positive] := polynomial[Expand[e^n]]

polynomial[e_ ^ q_Rational] :=
  polynomial[Expand[e ^ Numerator[q]]] /. x -> x^Denominator[q]

polynomial[e_ ^ x_ * f_ ^ x_] := polynomial[(e f)^x]

polynomial[q_Rational] := x - q

polynomial[e_ ^ n_Integer?Negative] := Block[{
	temp = polynomial[e ^ -n]},
       Expand[(x ^ degree[temp]) * (temp /. x -> 1/x)]]

polynomial[Sqrt[2]+Sqrt[3]]

polynomial[Sqrt[2] * Sqrt[3]]

polynomial[ (Sqrt[3] + Sqrt[2]) ^ (1/5) ]

polynomial[ (1 + 19 Sqrt[3]) ^ (1/3) * (Sqrt[2] + 5) ]

polynomial[ Sqrt[Sqrt[3] + Sqrt[2]]  *  (5 + I) ^ (1/3) ]

r = (1 + 1/Sqrt[2]) / (Sqrt[3] + 5 ^ (1/3)); polynomial[r]

% /. x -> N[r, 50]

10.6 The Riemann Zeta Function 
			(some cells take lots of time)

Zeta[Range[2, 16, 2]]

Plot[Zeta[n], {n,-10,10}]

Show[%, PlotRange->{{-10, 0}, {-.1, .1}}]

Plot3D[Abs[Zeta[x+I y]], {x, -4, 4}, {y, -10, 40},
	    PlotPoints->{70, 110}, ViewPoint->{8, 1, 3}, 
	    PlotRange->{0, 5}, Lighting->False,
	    Shading->False, Boxed->False,
	    BoxRatios->{5, 10, 2}, AxesLabel->{"x", "y", None},
	    AxesEdge->{{+1, -1}, Automatic, Automatic},
	    Ticks->{Automatic, Range[0, 30, 10], Range[0, 4]}]

Plot3D[1/Abs[Zeta[x+I y]], {x, -4, 4}, {y, -10, 40},
	    PlotPoints->{70, 110}, ViewPoint->{7, 2, 3}, 
	    PlotRange->{0, 5}, Lighting->False,
	    Shading->False, Boxed->False,
	    AxesLabel->{"x", "y", None},
	    AxesEdge->{{+1, -1}, Automatic, Automatic},
	    BoxRatios->{5, 10, 2},
	    Ticks->{Automatic, Range[0, 30, 10], Range[0, 4]}]

zeta[t_] := {Re[temp = Zeta[N[t]]], Im[temp]}
Attributes[zeta] = Listable

zeta[1/2 + I {0,1,2,3,4,5}]

ParametricPlot[zeta[.5 + I t], {t, 0, 26}]
ParametricPlot[zeta[.5 + I t], {t, 7004.1, 7005.32},
										MaxBend->20]

version 2.0 notes

The Z function will be built-in as RiemannSiegelZ[]

Z[t_] := Exp[I (Im[Log[Gamma[1/4 + I t/2]]] - t/2 Log[Pi])]*
		 Zeta[1/2 + I t]//N

Plot[Z[t]], {t, 0, 100}]

Z[100]

Z[10000]

Plot[Chop[Z[t], 10.^-6], {t, 7004.1, 7005.3},
	Ticks->{{7004.1, 7004.7,7005.3}, {-1.4, -1, -.6}},
	Axes->{7005, 0}]

Show[%, PlotRange->{{7005, 7005.15}, {-.1, .03}},
   	   Ticks->{{7005.03, 7005.08, 7005.13}, {-.1, -.02, .02}}]

approximatezeros = First[Transpose[ 
	{{14.1745, 0.00813295}, {20.9543, 0.00813295},
	 {25.0456, 0.00813295}, {30.4812, -0.00619712},
	 {32.936, -0.00619712}, {37.6117, 0.00813295}}]]

Chop[t] /. Map[FindRoot[Zeta[1/2 + I t], {t, # + {-.3, +.3}},
		AccuracyGoal->20, WorkingPrecision->30] &, 
          approximatezeros]

Z /@ {14.134725141734693790457251983500000000,
      14.134725141734693790457251983600000000,
      14.134725141734693790457251983700000000}

10.7 The Influence of the Complex Zeros of z on the Distribution of Primes  (some cells take lots of time)

N[LogIntegral[10^100], 50]

zeros = N[Re[t],20] /. 
	Map[FindRoot[Zeta[1/2 + I t], {t, # + {-.3, +.3}},
    	AccuracyGoal->20, WorkingPrecision->30] &, 
        		approximatezeros]

rho = 1/2 + I zeros;

list = Table[N[MoebiusMu[n]/n], {n, 154}];
MoebiusIndices = Select[Range[154], list[[#]] != 0 &]
MoebiusData    = -2 list[[MoebiusIndices]]

T[x_, k_] := MoebiusData .
  Re[  ExpIntegralEi[rho[[k]] / MoebiusIndices * Log[N[x]]]  ]

domain = Range[12., 100, 88./150]
T1 = T[domain, 1];

ListPlot[Transpose[{domain, T1}], PlotJoined->True]

start = 1
Do[T[domain, i] >>> "T.Data", {i, start, start + 3}]

(* Needs RiemannR from Appendix *)

RiemannData = RiemannR[domain] + N[ArcTan[Pi/Log[Domain]]/Pi]

following assumes data has been built as indicated in the text. This entire computation took me two weeks of computing two or three tables of T values each night!

Scan[ListPlot[Transpose[{domain, #}], PlotJoined->True,
		PlotRange->{{10, 100}, {5, 27}}, Axes->{10, 5},
		Ticks->{Range[20, 100, 10], Automatic}] &, 
	Accumulate[Plus, Prepend[data, RiemannData]]];