Appendix:

Supplementary Programs

Chapter 1

Wilson's Theorem Primality Test

WilsonPrimeTest[n_] :=
 Block[{i = 1}, Nest[Mod[# i++, n]&, 1, n-1] == n-1]

primeQ

BeginPackage["primeQ`"];

primeQ::usage = "primeQ[n] modifies PrimeQ by first testing for divisibility by any of the first 1,000 primes; it should be used instead of PrimeQ when numbers being tested are large";

Begin["`Private`"];

smallprimes = Prime[Range[1000]];
testforsmalldivisors[x_] :=
 (n = 1; 
  While[n < 1000 && Mod[x, smallprimes[[n]]] != 0, n++]; 
  n == 1000)

primeQ[n_] := (n <= 7919 || testforsmalldivisors[n]) && PrimeQ[n]

End[];   EndPackage[];

encode & decode (for version 1.2)

encode[message_] := ToExpression[
		Apply[StringJoin,
			Map[ToString[ToASCII[#] - 22]&,
				Characters[message] /. "z" -> "["]]]

decode[codedmessage_Integer] := Block[{out},
	out = Characters[codedmessage];
	out = Table[out[[{2i - 1, 2i}]], {i, Length[out]/2}];
	out = Map[Apply[StringJoin, #]&, out];
	out = Map[FromASCII[ToExpression[#] + 22]&, out];
	Apply[StringJoin, (out /. "[" -> "z" )]]

Version 2.0 versions of encode and decode

encode[message_] := ToExpression[ StringJoin @@ 
	ToString /@ (ToCharacterCode[message]-22 /. 100 -> 71) ]
	
decode[codedmessage_] := StringJoin @@ FromCharacterCode /@ 
	(22 + (ToExpression[StringJoin @@ #] & /@ 
	Partition[Characters[ToString @ codedmessage], 2, 2] /. 71 -> 100))

PrimePi (due to Ilan Vardi)

PrimePi::usage = "PrimePi[x] returns the number of primes 
less than or equal to x";

PrimePi[x_] := PrimePi[Floor[x]] /; Not[IntegerQ[x]]

PrimePi[n_Integer] := 0 /; n <= 1

PrimePi[4] = 2;

PrimePi[n_Integer] := Block[{t, diff, m},
	t = Round[N[Log[n]]];
	m = Round[N[LogIntegral[n]]];
	diff = n - Prime[m];
	While[ Abs[diff] > t^2 / (1.16)^2,
		m += Round[ diff / t];
		diff = n - Prime[m];] 
	If[diff < 0,
		While[n < Prime[m], m--;],
		While[n >= Prime[m], m++;]; m--];
	Return[m]]  /; n > 1 && n != 4

Attributes[PrimePi] = Listable;

RiemannR

RiemannR::usage = "Riemann[x] returns the value of Riemann's
prime-approximating function R(x).";

rlist = Map[1/(#! # N[Zeta[# + 1]]) &, Range[100]]
RiemannR[x_] := Round[1 + (Floor[N[Log[x]] ^ Range[100]]) . rlist]
RiemannR[0] = 0;
Attributes[RiemannR] = {Listable};

PrimePiModn

PrimePiModn[x_, n_] := Block[
	{i, temp, modlist = Mod[Prime[Range[x]], n]},
    Table[(temp = Flatten[Position[modlist, i]];
           Transpose[{Prime[temp], Range[Length[temp]]}]),
       {i, n - 1}]]

PrimePiMod4 exercise

ListPlot[Accumulate[Plus, Mod[Prime[Range[3000]], 4] - 2], 
	PlotJoined->True];

Chapter 2

Trochoid Animation Generator

r = Input["This program traces the path of a point on a spoke of a rolling wheel. The radius of the rolling wheel is 1. Enter the radius of the spoke."]

rr = Max[r - 1, 0]; pi = N[Pi];
cycloid[t_] := {t - r Sin[t], 1 - r Cos[t]};
rotate = {{Cos[theta],Sin[theta]}, {-Sin[theta],Cos[theta]}}
spokepoints :=
    Table[{t,1} + rotate . (cycloid[t] - {t,1}), {theta, 0, 2 pi, pi/4}]
spokes := {Thickness[.004], Map[Line[{{t,1}, #}] &, spokepoints]}
baseline = Line[{{-2,0}, {10.5,0}}];
dots = {PointSize[.006]};		
Do[  AppendTo[dots, Point[cycloid[t]]];
     Show[Graphics[{{GrayLevel[.5], Disk[{t,1}, 1]},
				Line[{{t,1}, cycloid[t]}],	(* spoke *)
				spokes, 
				baseline,
				Circle[{t,1}, 1],
				dots}],
		AspectRatio->(3+2 rr)/(12+2 rr),
		PlotRange->{{-1.5-rr, 10.5+rr}, {-.5-rr, 2.5+rr}}],
  {t, 0., 2 pi 59/40, 2 pi/40}]

Epicycloid and Hypocycloid Animation Generator

Clear[rotate]
r = Input["This program traces the path of a point on the circumference of a circle rolling around (positive r) or inside (negative r) a circle of radius 1. Enter the radius of the rolling circle."]
rr = Max[1, 1 + 2 r]; 
pi = N[Pi]
rotate[point_, theta_] :=
       {{Cos[theta], Sin[theta]}, {-Sin[theta], Cos[theta]}} . point
fixedcircle = {Thickness[.006], Circle[{0, 0}, 1]}
rollingdisk := {GrayLevel@.5, Disk[center, Abs[r]]}
rollingcircle := Circle[center, Abs[r]]
dots = {PointSize[.006]}

Do[  center = rotate[{-1-r, 0}, theta];
     spoke1 = rotate[{1,0}, theta (1+1/r)];
     spokeline1 = {center + r spoke1, center - r spoke1};
     spoke2 = rotate[{0,1}, theta (1+1/r)];
     spokeline2 = {center + r spoke2, center - r spoke2};
     spokes = {Thickness[.007], Line[spokeline1], GrayLevel[1],
 	                                         Line[spokeline2]};
     AppendTo[dots, Point[center + r spoke1]];
     Show[Graphics[
	{rollingdisk, spokes, fixedcircle, rollingcircle, dots}],
        AspectRatio->1,
        PlotRange->{{-rr - .03, rr + .03}, {-rr - .03, rr + .03}}],
  {theta, 0., 2 pi , 2 pi/60}]

Rolling Polygons

n = Input["How many sides does the rolling polygon have?"];
pi = Pi//N;
cr = Csc[pi/n];	(* circumradius of polygon *)
ir = Cot[pi/n];	(* inradius of polygon *)
a = ir ArcSinh[1/ir];
dots = {PointSize[.0018]};
rotate[p_, t_, center_] :=
           center + {{Cos[t], Sin[t]}, {-Sin[t], Cos[t]}} . (p - center) 
plot = Table[Plot[cr - ir Cosh[(x - 2 a i)/ir], {x, -a + 2 a i, a + 2 a i},
     DisplayFunction->Identity,
      PlotStyle->Thickness[.001]], {i, 0, 5}];

startpolygon = Map[# + {0, cr} &, Map[{Re[#], Im[#]} &,
	       cr N[Exp[I Range[-Pi/2 + Pi/n, 3 Pi/2 + Pi/n, 2 Pi/n]]]]];
	      

polygon[t_] := Block[{loop = N[Floor[(t + a)/(2a)]]},
					(* loop tells which loop the square is on *) 
Map[{t,0} + rotate[#, loop 2pi/n + ArcTan[Sinh[(t-loop 2 a)/ir]], {0,cr}]&,
   	startpolygon]]
   	
Do[	newpolygon = polygon[t]; 
	AppendTo[dots, Point[newpolygon[[n]]]];
	Show[
		Graphics[{PointSize[.008], Thickness[.001],
       			{GrayLevel[.75], Polygon[newpolygon]},
       			Line[newpolygon],
        			Map[Line[{{t, cr},#}] &, newpolygon],
        			Line[{{-cr, cr}, {11.5 a,cr}}],
        			Point[{t, cr}],
        			Point[newpolygon[[n]]],
        			dots}], plot,
    PlotRange->{{-cr, 11 a + cr}, {-.2, 2 cr}}, 
    AspectRatio->(2 cr +.2) / (11 a + 2 cr), 
    DisplayFunction->$DisplayFunction,
    Ticks->None, Axes->None],		{t, 0., 11 a, a/10}];

Reuleaux Triangle Rotating inside a Square

Clear[t, rotate]
centroid = {1 - 2 Sqrt[3.]/3, 0}
choose[t_] := Floor[.5 (Floor[t / (Pi/6)] + 3)] 
	
rotate[point_, t_, center_] :=
       center + {{Cos[t], -Sin[t]}, {Sin[t], Cos[t]}} . (point - center)

arc[{center_, startangle_}] :=
	Block[{t}, Table[center + 2 {Cos[t], Sin[t]},
					 {t, startangle, N[startangle + Pi/3], 0.08}]]

verts = {{1 - Sqrt[3.], -1}, {1, 0}, {1 - Sqrt[3.], 1}}
angles = {Pi/6, 5 Pi/6, 9 Pi/6};
dots = {};
Do[tempverts = Map[rotate[#, t, centroid] &, verts];	
   translate = 2 Mod[{choose[t], choose[t-Pi/2]}, 2] - 1 - 
				{tempverts[[Mod[choose[t], 3] + 1, 1]], 
				 tempverts[[Mod[choose[t - Pi/2], 3] + 1, 2]]};
   newverts = Map[ # + translate &, tempverts];
   reuleaux = Flatten[Map[arc, Transpose@{newverts, t + angles}], 1];
   AppendTo[reuleaux, First@reuleaux];
   dots = Join[dots, {newverts[[3]], centroid + translate,
                      rotate[{0,0}, t, centroid] + translate}];
   Show[Graphics[{{GrayLevel[.75], Polygon[reuleaux]},
                  Line[reuleaux],
                  Line[Append[newverts, First[newverts]]], 
                  PointSize[.02], Point /@ {
                      newverts[[3]],
                      centroid + translate, 
                      rotate[{0, 0}, t, centroid] + translate},
                  PointSize[.007], Point /@ dots,
                  {Thickness[.001], 
					Line[{{1,1},{1,-1},{-1,-1},{-1,1},{1, 1}}]},
                  Map[Line[{centroid + translate, #}] &, newverts]}],
     PlotRange->{{-1.05, 1.05}, {-1.05, 1.05}}, AspectRatio->1], 
 									{t, 0, 1.99 Pi, Pi/12}]

Cycloid as tautochrone

pi = Pi//N;
box = {Thickness[.003], Line[{{2 pi + .2, .2}, {2 pi + .2, -2.2},
                          {-.2, -2.2}, {-.2, .2}, {2 pi + .2, .2}}]};
cycloid[t_] := {t - Sin[t], -1 + Cos[t]};

cycloidplot = ParametricPlot[cycloid[s], {s, 0, 2 pi}, 
	Axes->None, PlotStyle->Thickness[.0005],
	DisplayFunction->Identity];

options = {PlotRange->{{-.2, 2 pi + .2}, {-2.2, .2}},
           AspectRatio->2.4/(2 pi + .4),
           DisplayFunction->$DisplayFunction};
stopline = {Thickness[.002], Line[{{pi, -2.2}, {pi, -1.8}}]};

showframe := Show[Graphics[{Thickness[.0004], 
	  Disk[{2 pi-posn[[1]], posn[[2]]},radius],
	  Disk[posn1,  radius],
	  stopline}], cycloidplot, options];

showstill := Show[cycloidplot,
   Graphics[{box, Thickness[.0004], list, stopline}], options];

g = 980 2 pi /15;

initial = Input["Enter initial t-value for start time on auxiliary cycloid: 0 < t < 3"] //N;
choice = InputString["Do you want the whole movie [movie] or just a single image [still]?"]
inc = Input["Enter the time increment, which controls the number of frames generated. Values between 1/100 and 1/400 are best."] //N

t = time = 0.; radius = .07; list = {};

(* MAIN LOOP BEGINS HERE. The variables posn and t refer to the case of the bead starting at the top of the cycloid; posn1 and t1 refer to the useràs starting position. *)

While[t < pi - radius/2,
	    t =  Sqrt[g]*time;
	    t1 = 2 ArcCos[Cos[initial/2.] * Cos[Sqrt[g] time / 2]];
	    posn =  cycloid[Min[t, pi - radius/2]];
	    posn1 = cycloid[Min[t1, pi - radius/2]];
	    time += inc;
	    If[choice == "movie", showframe,
			list = Join[list,
				{Disk[{2 pi - posn[[1]], posn[[2]]}, radius],
				 Disk[posn1, radius]}]]];
 If[choice != "movie", showstill]

Cycloid as brachistochrone

Clear[i, t, tt]
pi = Pi//N;
norm[p_] := Sqrt[p.p];
cycloid[t_] := {t - Sin[t], -1 + Cos[t]};
cycloidplot = ParametricPlot[cycloid[tt], {tt, pi, 2 pi}, 
  {Axes->None, PlotStyle->Thickness[.0005], DisplayFunction->Identity}];

degree = Input["Enter the degree (between 1 and 15) of the polynomial curve connecting the two points."];
n = -1 + Input["Enter the number of frames desired. Values between 10 and 30 are best."]//N;

g = 980 2 pi /15; inc = pi/Sqrt[g]/(3 n);

f[t_] := {t, -2 (t/pi) ^ (1/degree)};

directionf[t_] := Block[{v}, 
If[(degree > 1 && t != 0) || degree == 1, (v = f'[t])/norm[v], {0,-1}]];

choice = InputString["Do you want the whole movie [movie] or just a single image [still]?"];

options = {PlotRange->{{-.2, 2 pi + .2}, {-2.2, .2}},
   AspectRatio->2.4/(2 pi + .4), DisplayFunction->$DisplayFunction};

box = {Thickness[.003], Line[{{2 pi + .2, .2}, {2 pi + .2, -2.2},
                             {-.2, -2.2}, {-.2, .2}, {2 pi + .2, .2}}]};
stopline = {Thickness[.002], Line[{{pi, -2.1}, {pi, -1.9}}]};
showstill :=
     Show[cycloidplot, plotf, Graphics[{box, list, stopline}], options];
showframe := Show[cycloidplot, plotf, Graphics[{Thickness[.0004],
	  Disk[{2 pi-posn[[1]], posn[[2]]}, radius],
	  Disk[posn1, radius], stopline}],                  options];
	  
newposn[{x_,y_}, t1_] := FindRoot[
	f[tt][[2]] - y == -(tt - x) Divide @@ d, {tt, Max[t1, 10.^-60]},
	MaxIterations->100];

plotf = ParametricPlot[f[t], {t, 0, pi}, Axes->None,
               PlotStyle->Thickness[.0005], DisplayFunction->Identity];
radius = .07; posn = posn1 = {0,0}; t1 = 0; d = directionf[0];
list = {Thickness[.0004],	Disk[{2 pi, 0}, radius], Disk[posn1, radius]};

If[choice == "movie", showframe];

Do[ t = pi i / (3 n);
	{x, y} = posn1 + d*inc*(Sqrt[-2 g posn1[[2]]] - g d[[2]] inc/2);
	
	(* Preceding line computes the point on the tangent line that the
		bead would roll to; the following line sends this position to 		newposn to get the corresponding point on the curve *)
		
	t1 = If[degree == 1, x, tt /. newposn[{x,y}, t1]];
	d = directionf[t1];
	posn  = cycloid[Min[t, pi - radius/2]];
	posn1 = f[Min[t1, pi - radius]];

	If[Mod[i, 3] == 0, If[choice == "movie", showframe,
	   list = Join[list, {Disk[{2 pi - posn[[1]], posn[[2]]}, radius],
	                      Disk[posn1, radius]}]]],
  {i, 3 n}];
  
If[choice != "movie", showstill];

Chapter 4

Cantor Set via Intervals

spawn1[{a_, b_}] := Block[{w = (b - a)/3}, {{a, a + w}, {b - w, b}}] /;
	           (Head[a] == Head[b] == Rational) || a==0 || b == 1;

spawn1[intervals_List] :=
        Flatten[Map[spawn1, intervals],1] /; Length[intervals[[1]]] > 1;

lines[n_Integer, l_List] :=
 			Map[Line[{{#[[1]], n + .5}, {#[[2]], n + .5}}] &, l];
level = 2;
Show[Graphics[{Thickness[.0005],
	Map[lines[level--, #] &, NestList[spawn1, {{0,1}}, level]]}],
	 	         Ticks->{N[Range[0,1,1/9],3], None}, Axes->Automatic]

Animation of Cantor functions

Clear[spawn, f];
spawn[{a_Rational, b_Rational}] :=
	Block[{w = (b - a)/3}, {{a - 2w, a - w}, {b + w, b + 2w}}];
spawn[intervals_List] :=
        Flatten[Map[spawn, intervals], 1] /; Length[intervals[[1]]] > 1;

removedintervals[n_] := Flatten[NestList[spawn, {{1/3, 2/3}}, n], 1];

connect[{a_, b_}] := Block[{y = f[a]}, Line[{{a, y}, {b, y}}]];

f[0]   = 0;
f[x_] := f[3x]/(r + 1)       /;  x <= 1/3;
f[x_] := 1 - r f[1 - x]      /;  2/3 <= x;

CantorFunction[parameter_] := (r = parameter;
			Show[Graphics[{Thickness[.001],
			Map[connect, removedintervals[6]]}],
				Axes->{0, 0},
				PlotRange->{{0,1}, {0,1}},
				Ticks->{N[Range[0., 1, 1/9], 3], Automatic}]);

Attributes[CantorFunction] = {Listable}

CantorFunction[{.01, .1, .3, .5, 1, 4, 10, 100}];

The Sierpinski Triangle

Clear[spawn]
spawn[{a_, b_, c_, a_}] :=
	Block[{ab = (a+b)/2, bc = (b+c)/2, ac = (a+c)/2},
	{{a, ab, ac, a}, {b, bc, ab, b}, {c, ac, bc, c}}];

spawn[l_List] := Flatten[Map[spawn, l], 1] /; Length[l[[1]]] > 1;

Show[Graphics[{Thickness[.001], 
  Line /@ Nest[spawn, {{0, 0}, {.5, Sqrt[3.]/2}, {1, 0}, {0, 0}}, 7]}],
	 AspectRatio->Sqrt[3]/2]

Chapter 5

Chaos Game Using a Universal Sequence

top = {.5, Sqrt[3.]/2}
f1[x_] := .5 x;     f2[x_] := .5(x + {1,0});     f3[x_] := .5(x + top)

list = Append[Table[0, {8-1}], 1];

Show[Graphics[Map[Point,
	NestList[(   list[[1]] = Mod[list[[1]] - list[[6]], 3];
		        list = RotateLeft[list];
		        {f1[#], f2[#], f3[#]}[[Last[list]+1]] ) &, 
		    {0,0}, 3^8-1]]],
PlotRange->{{0, 1}, {0, .87}}, Axes->Automatic, AspectRatio->.87,
Ticks->{{.25, .5, .75, 1}, {.25, .5, .75, 1}}]

Filled Julia Sets for the Real Quadratic Map

FilledJuliaSetReal::usage =
"FilledJuliaSetReal[r, meshx, meshy, x,  y, iters] shows the filled-in Julia set for r*z*(1-z), by using the Escape Time Algorithm on a meshx-by-meshy grid over the rectangle where x varies from 0 to x and y from 0 to y. Fourfold symmetry is then used to generate an image where x ranges from 0 to 1 and y from -y to y. Thus x should be set to 0.5 to get the entire set. The number of iterations used is given by iters. A Print statement in the subroutine orbitcheckreal acts as a progress checker."

orbitcheckreal[z_, r_, iters_] := (s = z; i = 0;
     If[Re[z] != xold, xold=Re[z]; Print[xold]];
     While[++i <= iters && Abs[s *= r (1-s)] < r0];
	If[i != iters + 1, {}, 
      a = Re[z]; b = Im[z]; Flatten[Outer[List, {a, 1-a}, {b, -b}], 1]])

FilledJuliaSetReal[r_, meshx_, meshy_, x_,  y_, iters_] := 
     (r0 = 1 + 1/Abs[r]; xold = 0;
     Show[Graphics[{PointSize[.002],
	  Point /@ Flatten[
	    Outer[orbitcheckreal[#1 + I #2, N[r], iters]&,
                      Range[0, x, x/meshx], Range[0, y, y/meshy]], 2]}],
	PlotRange->All, AspectRatio->y/x])

Tree Animations

Tree::usage = "Tree[depth, length, angle, ratio] uses the recursive
subroutine, branch, to generate and display a binary tree with small
round leaves at the end of each terminal branch. "

depth::usage = "The depth of the recursion of the tree function.  The
tree will have 2^depth terminal branches and leaves."

length::usage = "The original length of the trunk of the tree in the
tree function.  This argument determines the overall size of the tree."

angle::usage = "The whole angle, in degrees, between paired subbranches
of a single parent branch.  This argument determines how much the tree
spreads out."

ratio::usage = "The reduction ratio which determines how much smaller
a subbranch is than its parent branch.";

Attributes[Tree] = Listable

branch[depth_Integer, length_, angle_, ratio_] := (
	forward[length];
	left[angle/2];
	If[depth > 0, branch[depth - 1, length*ratio, angle, ratio],
			    AppendTo[points, X]];
	right[angle];
	If[depth > 0, branch[depth - 1, length*ratio, angle, ratio]] ;
	left[angle/2];
	back[length];)
		
Tree[depth_Integer, length_, angle_, ratio_] := (
	initialize[];
	points = {};  (* to be a list of dots at the branch-ends *)
	left[90];
	If[depth > 0, branch[depth, length, angle, ratio]];
	finished;
	Show[Graphics[{Line[path], PointSize[0.03], Point /@ points}],
		PlotRange->All])

Tree[5, 1, Range[45, 90, 2.5], .8];   (* Tree pushups *)

The following command generates a nice animation if the PlotRange setting in the definition of Tree is changed to 
{{-2.4, 2.4}, {0, 3.5}}. The sequence of frames should be animated in the forward direction

Tree[5, 1, 90, Range[.05, .8, .05]];  (* Tree growth *)

Chapter 6 — UltraTurtle and RecursiveTurtle3D placed in  Chapter 6 file.

Sierpinski and original Peano space-fillers

SierpinskiCurve[n_] := Show[RecursiveTurtle[
				{"X"-> "XF-F+F-XF+F+XF-F+F-X"}, "F+XF+F+XF",
				 n, 90, side = Sqrt[2.] 2^-(n+2),
			 	({temp = ((2^(n-1)6 + 2^(n) -2) side /Sqrt[2])};
(*stepsize*)		  {(temp+1)/2 - side/Sqrt[2], (1-temp)/2}),
(*start direction*)	Sqrt[{2,2}]/2,		
(*path length*)	(4^(n+2) - 1)/3	],
	Axes->{0, 0}, PlotRange->{{0, 1}, {0, 1}},
	Ticks->{{.25, .5, .75, 1}, {.25, .5, .75, 1}}]
		
OriginalPeanoCurve[n_] := Show[RecursiveTurtle[
 	{"X"-> "XFYFX+F+YFXFY-F-XFYFX", "Y"->"YFXFY-F-XFYFX+F+YFXFY"},
 	"X", n, 90, 3^-n, {3, 3}^-n / 2, 9^n],
 	Axes->{0,0}, PlotRange->{{0, 1}, {0, 1}},
 	Ticks->{{.33, .67, 1}, {.33, .67, 1}}]

Illustration of a Function from the Interval Onto the Square

CantorBijection::usage = "CantorBijection[t] gives the value of a bijection from the unit interval to the unit square based on even- and odd-indexed binary digits. The input and outputs are base-10 numbers in the unit interval and unit square, respectively.";

CantorImage::usage = "CantorImage[n] gives the path in the unit square obtained by using straight lines to connect the values of CantorBijection on {0, 1/n, 2/n, . . . , 1}.";

SetAttributes[{CantorBijection, CantorImage}, Listable];

base2digits[t_] := Block[{tt = N[t]},
		Map[If[2. ^ -# <= tt, tt -= 2. ^ -#; 1, 0] &,  Range[30]]]

CantorBijection[t_] := (
	temp = Flatten@Position[base2digits[t], 1];
	temp = {Select[temp, !EvenQ[#]&], Select[temp, EvenQ]};
	{Plus @@ (2^((-temp[[1]]-1)/2)),
	 Plus @@ (2^(-temp[[2]]/2)) } //N )

CantorBijection[1] = {0,1};

CantorImage[n_] :=
	Show[Graphics[Line[CantorBijection[Range[0, 1, 1/n]]]],
		PlotRange->{{0, 1}, {0, 1}}, AspectRatio->1,
		Axes->Automatic, AxesStyle->GrayLevel[1],
		Ticks->Range[{0, 0}, {1, 1}, {.25, .25}]]

CantorImage[13]  (* n should be a power of 2 *)

Chapter 8

Test for Whether Integers Are Pairwise Relatively Prime

PairwiseCoprime::usage = "PairwiseCoprime[mlist] returns True if the integers in mlist are pairwise relatively prime, False otherwise."

PairwiseCoprime[mlist_List] :=
	Length[mlist] == Length[Union[mlist]] &&
	And @@ Flatten[Outer[GCD[#1,#2] == 1 || #1 == #2 &, mlist, mlist]]

The Splitting Algorithm for Egyptian Fractions

EgyptianSplitting::usage = "EgyptionSplitting[q] uses the splitting method to write q as a sum of distinct unit fractions"

EgyptianSplitting[q_] := (
  list = Table[Denominator[q], {Numerator[q]}];
  While[Length[list] != Length[Union[list]], list =
	Sort[Flatten[
		Map[(t = list[[#]];
               If[{#} != First[Position[list, t]], {t+1, t(t+1)}, t]) &, 
 		   Range[Length[list]]]]]];
  list)

EgyptianSplitting[4/7]

Code to Generate Table 8.2

Map[Length[Digits[Last[#]]]&,
    EgyptianFraction[(10 ^ Range[9] - 1) / 10 ^ Range[9]]]

Chapter 9 routines placed in Chapter 9 file