Supplementary Materials

Here again is the distance function:

distance[pt1_, pt2_] := Sqrt[Total[(((pt1 - pt2))^2)]]

points = Table[Random[], {5}, {2}]

Probably the simplest way to solve the problem is to use the built-in Subsets function to obtain all two-element subsets of our point set.

pairs = Subsets[points, {2}]

Length[%]

Here then is the new diameter function using Subsets.

diameter[lis_] := 
      Max[Apply[
          distance[(((#1 - #2))^2)] &, 
Subsets[
            lis, {2}], {1}]]

And here is a version that uses pure functions and avoids the need for the auxiliary distance function.

diameter[lis_] := 
      Max[Apply[
          Sqrt[Total[(((#1 - #2))^2)]] &, 
Subsets[
            lis, {2}], {1}]]

Finally here are some timing tests to check the efficiency of this improved code against that given in the solutions section of the book.

diameterOrig[lis_] := 
      Max[Apply[
          Sqrt[Apply[
                Plus, (((#1 - #2))^2)]] &, 
Distribute[{lis, lis}, List], {1}]]

points = Table[Random[], {1000}, {2}];

Timing[diameter[points]]

Timing[diameterOrig[points]]

From a comparison of (n^2) with (({{n}, {2}})) we would expect about a twofold increase in computational speed and we are pretty close to that.

(n^2)/Binomial[n, 2] // Simplify

%[[1]]/%%[[1]]

Actually, for large point-sets (and more so for point sets in higher dimensions), this naive approach to computing the diameter can bog down very quickly due to the quadratic nature of the underlying computation. An important result in computational geometry (Hocking and Young 1961) states that the diameter of a point set is equivalent to the diameter of its convex hull. Since very efficient and fast routines have been developed for computing the convex hull, it seems reasonable to try and incorporate this fact into a computation of the diameter.

This loads the package containing the ConvexHull definitions.

<< DiscreteMath`ComputationalGeometry`

The ConvexHull function returns a list of the positions from our list of those points that are in the convex hull.

convexhull = ConvexHull[points]

Here are the explicit points in the convex hull.

points[[convexhull]]

And here is the diameter of the convex hull, agreeing with the earlier value.

diameter[%]

Putting these pieces together, here is a function we can use for timing comparisons.

diameterCH[lis_] := 
      Module[{convexhull}, 
convexhull = 
          ConvexHull[lis]; 
diameter[
          lis[[convexhull]]]]

Timing[diameterCH[points]]

Let’s try a larger example now.

points = Table[Random[], {2500}, {2}];

Timing[diameter[points]]

Timing[diameterCH[points]]

Pretty fast! In fact, using the convex hull, we can compute the diameters of quite large point-sets and see how the diameter for points in the unit square approaches the limiting value of (Sqrt[2]).

points = Table[Random[], {(10^6)}, {2}];

Timing[diameterCH[points]]

A functional approach

With slight modifications, here is the function as presented in the text on page 110. In this implementation we give an integer n as argument to Survivor and compute the Josephus iteration for the list of numbers 1 through n.

Survivor[n_] := Nest[Rest[RotateLeft[#]] &, Range[n], n - 1]

Survivor[10]

Here is a listing of the survivors, each row representing a round of execution (beginning with the original orientation).

SurvivorList[n_] := 
      NestList[Rest[RotateLeft[#]] &, Range[n], n - 1]

SurvivorList[10] // TableForm

The following functions generalize this problem by considering groups of n people, executing every m-th.

Survivor::"usage" = "Survivor[n,m] determines who will be the \
last to live of a group of n people, numbered 1 to n if every remaining mth \
person were eliminated until only one survives.  The function outputs the \
lone survivor after the requisite executions.";

SurvivorList::"usage" = "SurvivorList[n,m] determines who will \
be the last to live of a group of n people, numbered 1 to n if every \
remaining mth person were eliminated until only one survives.  The output \
consists of successive lists  containing the survivors after each \
execution.";

Survivor[n_, m_] := 
      Nest[Drop[RotateLeft[#], m - 1] &, Range[n], n - 1]

Survivor[10, 2]

SurvivorList[n_, m_] := 
      NestList[Drop[RotateLeft[#], m - 1] &, Range[n], n - 1]

SurvivorList[10, 2] // TableForm

A recurrence relation

In the case of an odd number of people, person number 1 is to die immediately after the nth. A bit of thought should convince you that the following recurrence defines the Josephus function for n people, assuming that every second person is eliminated.

SurvivorRec[1] = 1;
SurvivorRec[n_?EvenQ] := 
      2  SurvivorRec[(n/2)] - 1 /; 
        n ≥ 1
SurvivorRec[n_?OddQ] := 
      2  SurvivorRec[((n - 1)/2)] + 1 /; n ≥ 1

Here now is a table containing n and the Josephus number of n, 1 ≤ n ≤ 10:

TableForm[Table[{i, SurvivorRec[i]}, {i, 1, 10}]]

This is a very efficient recurrence because it reduces n by a factor of 2 each time through. Note how quick it is in computing SurvivorRec[(10^9)].

Timing[SurvivorRec[(10^9)]]

Looks like the 926,258,177th person will survive.

A closed form

Hint: Consider any number n in the form n = 2 + q where 2¸ is the largest power of 2 less than or equal to n, and q is what’s left over.

f[n_] := n - (2^Floor)[N[Log[2, n]]]

f[15]

Survivor2[n_] := 2 (n - (2^Floor)[N[Log[2, n]]]) + 1

Survivor2[100]

Survivor2[(10^9)]

Cyclic shifts

The Josephus number of 100 is 73.

Josephus[100, 2]

First let’s put 100 in base 2,

BaseForm[100, 2]

and then get a list of its digits.

IntegerDigits[100, 2]

Do a binary shift of one unit,

l = RotateLeft[IntegerDigits[100, 2], 1]

and convert back from base 2 to base 10:

r = Reverse[l]

Dot[r, (2^Range)[0, Length[r] - 1]]

SurvivorBinary[n_] := 
      Module[{p}, 
p = 
          Reverse[RotateLeft[IntegerDigits[n, 2], 1]]; 
Dot[
          p, (2^Range)[0, Length[p] - 1]]
]

SurvivorBinary[100]

Timing[SurvivorBinary[(10^9)]]

Collatz[n_?EvenQ] := (n/2)

Collatz[n_?OddQ] := 3  n + 1

CollatzSequence[n_] := 
NestWhileList[Collatz, 
        n, # ≠ 1 &]

CollatzSequence[11]

An alternative to this approach uses a fourth argument to NestWhileList in such a way that it will detect cycles in a sequence. The form to use is the following:

NestWhileList[f,expr,UnsameQ,All]

NestWhileList will iterate the function f, starting with initial value expr, and continue applying f until a repeated output appears. Here is the implementation for CollatzSequence.

CollatzSequence[n_] := 
NestWhileList[Collatz, n, 
        UnsameQ, All]

CollatzSequence[27]

Note that this implementation outputs the repeated number in the sequence before terminating. You could do some post-processing of the output (Drop[…,-1] or Most[…]) if you wished to omit that number.

First we identify the directory where the notebook(s) that we wish to work with are located.

\($BookDirectory = 
      ToFileName[{$BaseDirectory, "Applications", "IPM3", "Documentation", "English"}]\)

Here is the code to count the number of cells in a particular notebook of a given style. The CountStyles function performs the following operations: open notebook (NotebookOpen), read it into the kernel (NotebookGet), and then count the number of cells of a particular style (Length[Cases[…]]). The open notebooks are closed using a front end token (FrontEndTokenExecute[…,Close]. Tokens allow you to carry out actions on the command line that are normally performed through the menu system. We did not cover them in the book, but the interested reader should consult the documentation on Front End Tokens.

CountStyles[nb_String, style_String] := 
      Module[{nbExp, nbb}, 
nbb = 
          NotebookOpen[
            ToFileName[{$BookDirectory}, nb]]; 
nbExp = 
          NotebookGet[nbb]; 
FrontEndTokenExecute[
          nbb, "Close"]; 
        Length[Cases[nbExp, Cell[__, style, ___], 
            Infinity]]
]

Now let’s use CountStyles on one notebook from $BookDirectory to count the number of occurrences of the ExerciseItem cell style. Note: ExerciseItem is a style that is particular to the IPM3Styles.nb style sheet and is not part of any style sheet in the regular Mathematica distribution.

CountStyles["12Packages.nb", "ExerciseItem"]

In order to count cell styles for all the notebooks in a particular directory, we will first set the current working directory.

SetDirectory[$BookDirectory]

Next we pick out all those files for which we want to count cell styles.

filelist = Take[FileNames[], {3, 13}]

Map CountStyles over the files from filelist (this may take a while depending upon how many notebooks you are processing and how large they are).

Map[CountStyles[#, "ExerciseItem"] &, filelist]

Add 'em up.

Total[%]

So there are 242 cells of the ExerciseItem style in the given directory; in other words, 242 exercises in total in this directory of notebooks.

First we create a path to the style sheet we wish to work with.

\(styleSheet = 
      ToFileName[{$InstallationDirectory, "SystemFiles", "FrontEnd", "StyleSheets"}, "Default.nb"]\)

The function menuStyleNames creates a list of the style names from this particular style sheet. Basically it strips out all of the StyleData information leaving us with the style name itself.

menuStyleNames = 
Cases[Get[styleSheet], 
        Cell[StyleData[str_String], opts___] :> 
          str /; ! MemberQ[{opts}, 
                StyleMenuListing -> None], ∞]

You might want to omit certain styles from your palette. This can be done using DeleteCases, giving an explicit list of those styles that you don’t want on your palette when it is done.

\(menuStyleNames = 
      DeleteCases[
        menuStyleNames, "Outline1" | "Outline2" | "Outline3" | "Outline4"]\)

Here is the code that will create the palette from a list of menu style names. Note the CellTag so that we can refer to this cell when we generate the free-floating palette below. We have also added some formatting of the text on the buttons (FontFamily, FontSize, etc.)

\(styleButtonCell[lis_] := 
Cell[
        BoxData[GridBox[\({ButtonBox[#, 
                    ButtonData :> #]} &\)\  /@ \ 
              lis, 
RowSpacings\  -> \ 
              0, 
ColumnSpacings\  -> \ 
              0]], 
"Text", 
FontFamily -> "Helvetica", 
        FontColor -> Purple, 
FontSize -> 9, 
        FontWeight -> "Bold", 
CellTags\  -> \ "style palette", 
ButtonBoxOptions\  -> \ {
ButtonEvaluator\  -> \ 
              None, 
ButtonFunction :> \((FrontEndExecute[{
FrontEndToken[
                        InputNotebook[], "Style", #2]}] &)\)
}]\)

Turn the menuStyleNames into a palette.

styleButtonCell[menuStyleNames] // CellPrint

({{Title}, {Subtitle}, {Subsubtitle}, {Section}, {Subsection}, {Subsubsection}, {Text}, {SmallText}, {Input}, {Output}, {InlineFormula}, {DisplayFormula}, {Program}})

Here we use a front end token to generate the palette. Tokens essentially give you a way to programmatically perform actions that are normally associated with menu items. In this case, the FrontEndTokenExecute["GeneratePalette"] is equivalent to using the menu item File⊳Generate Palette from Selection.

NotebookLocate["style palette"];
FrontEndTokenExecute["GeneratePalette"]

Apply[Plus, {a, b, c, d, e}]

Total[{a, b, c, d, e}]

Fold[Plus, 0, {a, b, c, d, e}]

But what if we needed to accumulate the partial sums? (This is particularly useful for creating random walks, for example.) One way to do this would be to use FoldList instead of Fold above.

FoldList[Plus, 0, {a, b, c, d, e}]

We would then need to toss out the first element.

Rest[%]

Let us create a function, Accumulate, that puts these few pieces together.

Accumulate[lis_List] := Rest[FoldList[Plus, 0, lis]]

Accumulate[{a, b, c, d, e}]

Now suppose we wished to generalize this to higher dimensions. (Imagine that we were using this function to create random walks in two and three dimensions.) So in other words, starting with a list like {{a,b},{c,d},{e,f}}, we would like to end up with the following partial sums.

{{a,b}, {a+c,b+d}, {a+c+e,b+d+f}}

We can follow what we did in one dimension but we need to adjust the initial expression that we start with.

Accumulate[lis_List] := 
      Rest[FoldList[Plus, Table[0, {Length[lis[[1]]]}], lis]]

Accumulate[{{a, b}, {c, d}, {e, f}}]

In fact, this generalizes just fine to higher dimensions as well.

Accumulate[{{a, b, c}, {d, e, f}, {g, h, i}}]

Unfortunately, this breaks the code for one-dimensional lists.

Accumulate[{a, b, c, d, e}]

We can fix that by adding a condition for the case when the argument to Accumulate is a one-dimensional list. We end up with two rules – one for the one-dimensional case, and another for higher-dimensional cases.

Clear[Accumulate]

Accumulate[lis_List /; Length[Dimensions[lis]] == 1] := 
      Rest[FoldList[Plus, 0, lis]]

Accumulate[lis_List] := 
      Rest[FoldList[Plus, Table[0, {Length[lis[[1]]]}], lis]]

Accumulate[{a, b, c, d, e}]

Accumulate[{{a, b, c}, {d, e, f}, {g, h, i}}]

randomPermutation[lis_] := 
      Module[{res = {}, x, 
          l2 = lis}, 
Do[
x = 
            Part[l2, 
              Random[Integer, {1, Length[l2]}]]; 
res = 
            Append[res, x]; 
l2 = 
            Complement[l2, {x}], 
{i, 1, 
            Length[lis]}]; 
res]

randomPermutation[Range[10]]

randomPermutation[{a, b, c, d, e}]

Here is another procedural approach that avoids Append by making the list assignment inside the Do loop.

randomPermutation2[lis_List] := 
      Block[{r, l = lis, len = Length[lis]}, 
Do[
          r = {i, Random[Integer, {i, len}]}; 
l[[r]] = 
            l[[Reverse[r]]], 
{i, 
            len}]; 
l]

randomPermutation2[{a, b, c, d, e}]

The first approach is very slow due to the fact that Append creates copies of the list it is working with and does so every time through the loop.

lis = Range[(10^4)];

Timing[randomPermutation[lis];]

Timing[randomPermutation2[lis];]

Here is a functional approach that uses the Ordering function. First we create a list of random numbers of the same length as our list that we will be working with (see below).

Table[Random[], {5}]

This gives the positions in the above list where each element would appear after Sort was run on this list.

ord = Ordering[%]

And here is our permutation on a short list.

lis = {a, b, c, d, e};
Thread[Part[lis, ord]]

Putting the pieces together…

randomPermutationF[lis_List] := 
      Thread[Part[lis, Ordering[Table[Random[], {Length[lis]}]]]]

randomPermutationF[Range[10]]

Let’s do some timing comparisons of this functional approach with the second procedural approach above.

lis = Range[(10^5)];

Timing[randomPermutation2[lis];]

Timing[randomPermutationF[lis];]

lis = {a, b, c, d, e};

Part[lis, Random[Integer, {1, Length[lis]}]]

Or, using shorthand notation for Part:

lis[[Random[Integer, {1, Length[lis]}]]]

Here then is a function, RandomChoice, that takes an arbitrary list as an argument and chooses an element at random from that list.

RandomChoice[lis_List] := 
      Part[lis, Random[Integer, {1, Length[lis]}]]

RandomChoice[{a, 3, "bob", π, f[x]}]