Notebook for functional programming talk (March 2020)

Disclaimer: these notes presents a functional/rule-based view of Mathematica. Mathematica also supports iterative constructions.

Always press shift-return to evaluate. If you press return you will end up with a compound statement which may not evaluate as you expect. Compound statement are best written using parentheses and semicolons, e.g. (y:=3; y+1) returns 4. They are not often necessary.

If inexplicable errors continue: try quitting the kernel (Evaluation menu at bottom). This clears out all stored definitions. A new kernel is started automatically when you next evaluate something.

Fibonacci function

Fib[0] := 0

Fib[1] := 1

Fib[n_]:=Fib[n-1]+Fib[n-2]

Jumping ahead, Table is useful for getting lists.

?Table

Table[{n,Fib[n]},{n,0,10}]

Here // means `apply the next function’: if helpful, please think of it as like the pipe | in Linux. Try typing ?// to see Mathematica’s documentation for it

The TF functions below is very useful to define for printing 2-dimensional arrays.

TF[xs_]:=TableForm[xs,TableDepth->2]

Table[{n,Fib[n]},{n,0,10}]//Transpose//TF

Problem (1): meaningless results on non-natural numbers – see slide on patterns for /;

Module[{$RecursionLimit = 20},Fib[1.5]]

FibSafe[0]:=0; FibSafe[1] := 1; FibSafe[n_Integer/;n≥0] := FibSafe[n]=FibSafe[n-1]+FibSafe[n-2]

Now nothing matches on 1.5, so Mathematica just returns the input

{FibSafe[1.5],FibSafe[-3]}

Quiz: how should Fib be defined on negative numbers? How might it be defined on general real numbers?

FibByFormula[n_]:=Module[{Phi = (1+Sqrt[5])/2, Psi = (1-Sqrt[5])/2},1/Sqrt[5]*(Phi^n-Psi^n)]//Simplify

Table[{n,FibByFormula[n]},{n,-10,10}]//Transpose//TF

Plot[{Re[FibByFormula[x]],Im[FibByFormula[x]],Abs[FibByFormula[x]]},{x,-5,7},PlotStyle->{Red,Blue,Black}]

Problem (2): slow – see final slide for how memoization works in Mathematica

AbsoluteTiming[Fib[30]]

FibMemoized[0]:=0; FibMemoized[1] := 1; FibMemoized[n_] := FibMemoized[n]=FibMemoized[n-1]+FibMemoized[n-2]

AbsoluteTiming[FibMemoized[30]]

Collatz function using pattern guard /;

Collatz[x_]/;EvenQ[x]:=x/2

Collatz[x_]/;OddQ[x]:=3*x+1

?NestList

NestList[Collatz,23,15]

Exercise: define a function that returns the number of iterations to reach 1

CollatzIterations[1]:=0

CollatzIterations[x_]:=1+CollatzIterations[Collatz[x]]

Map[CollatzIterations,Range[1,100]]

Aside: Mathematica can be told to repeatedly apply a rule until there is no change. The pattern below matches x to the head (first element) of the list and xs to the rest. To type :> use :>

{23}//.{x_,xs___}/;x>1:>{Collatz[x],x,xs}

Here is CollatzIterations written using rules

{0,25}//.{i_,x_}/;x>1 :>{i+1,Collatz[x]}

Table[{n,0,n}//.{n_,i_,x_}/;x>1 :>{n,i+1,Collatz[x]},{n,1,20}]//Transpose//TF

Here are two version of Fib written using rules

FibByPattern[n_]:=n/.{0->0,1->1,n:>FibByPattern[n-1]+FibByPattern[n-2]}

FibByRepeatedPattern[n_]:={n,0,1}//.{{0,a_,b_}:>a,{c_,a_,b_}:>{c-1,b,a+b}}

Table[{n,FibByPattern[n],FibByRepeatedPattern[n]},{n,0,10}]//Transpose//TF

Most specific pattern wins

g[1]:=1

g[x]:=x+1

g[y_]:=y/2

g[{x_,y_}]/;EvenQ[y]:=y/2

Patterns are useful for quickly restructuring data

Table[{n,Fib[n]},{n,0,10}]/.{x_,y_}:>{x+1,2*y}

Digit sum

Digits[x_]/;x<10:={x}

Digits[x_]:=Reverse[{Mod[x,10]}~Join~Digits[Quotient[x,10]]]

Shorthand: f@@xs means apply f with arguments taken from the list xs: e.g. Plus @@ {1,2,3} becomes Plus[1,2,3], evaluating to 6.

DigitSum[x_]:=Plus@@Digits[x]

We use the `apply pattern until no change’ trick from the Collatz function: you could also use a recursion with two cases, by conditioning on x < 10.

RepDigitSum[x_]:=x//.{y_:>DigitSum[y]}

Quiz (for the very brave): why does RepDigitSum fail if y is replaced with x above? Is this a bug or a feature?

SquareDigitSumsUpTo[n_]:=Table[{m,RepDigitSum[m^2]},{m,1,n}]

SquareDigitSumsUpTo[20]//Transpose//TF

GroupIn (just to make printing nicer)

GroupIn[k_,xs_]:={xs}/;Length[xs]≤k

GroupIn[k_,xs_]:=Join[{xs[[1;;k]]},GroupIn[k,xs[[k+1;;]]]]

Quiz: explain the periodic pattern, i.e. why is the repeated sum of the digits  the same for n^2 and (n + 9)^2? [Hint: this is easier than it might seem: try evaluating Table[RepDigitSum[x], {x,1,100}] to see what happens without the square.]

GroupIn[9,#[[2]]&/@SquareDigitSumsUpTo[45]]//TF

Using Mathematica to explore binomial identities

Table[Binomial[n,m],{n,0,10},{m,0,n}]//TF

Alternating sum along rows, for instance 1-10+45-120 = -84.

Table[Sum[(-1)^k*Binomial[n,k],{k,0,m}],{n,0,10},{m,0,n}]//TF

Mathematica has powerful algorithms due to Gosper, Zeilberger, ... that allow it to simplify many sums involving binomial coefficients

Sum[(-1)^k*Binomial[n,k],{k,0,m}]

Quiz: is there a simple formula for the sum of the first m entries in each row of Pascal’s Triangle? What does Mathematica give?

Toy RSA-Cryptosystem

p:=NextPrime[2^32+2^31]; q:=NextPrime[2^32+2^30]; M := p*q; eBad := 2;

Quiz: why not use N for pq? [Hint: Mathematica specific. Similarly there is now a built-in ‘Encrypt’ functions, so we use a different name below.] Why is using 2 as the encryption exponent, i.e. x -> x^2 bad?

e:=11;d := PowerMod[11,-1,(p-1)*(q-1)];

ToyEncrypt[x_]:=Mod[x^e,M]

ToyDecryptBad[y_]:=Mod[y^d,M]

ToyEncrypt[1000]

Quiz: now try decrypting: the problem is that Mathematica computes y^d (a huge number) before doing any reduction mod M

ToyDecryptBad[ToyEncrypt[1000]]

A better implementation uses PowerMod, as used to find d.

ToyDecrypt[y_] := PowerMod[y,d,M]

ToyDecrypt[ToyEncrypt[1000]]

Quiz: evaluate in your head

(1000+M)//Encrypt//Decrypt

Exercise: how to write PowerMod using just Mod

PM[x_,0,n_]:= 1

PM[x_,e_,n_]/;EvenQ[e]:= PM[Mod[x^2,n],e/2,n]

PM[x_,e_,n_]/;OddQ[e] := Mod[x*PM[Mod[x^2,n],(e-1)/2,n],n]

ToyDecryptPM[y_]:=PM[y,d,M]

ToyDecryptPM[ToyEncrypt[1000]]

The definition of PM might be more readable using a variable to store Mod[x^2,n]. Below we use a compound expression (first definition; second definition).

PM[x_,0,n_]:= 1

PM[x_,e_,n_]/;EvenQ[e]:= (y=Mod[x^2,n];PM[Mod[y,n],e/2,n])

PM[x_,e_,n_]/;OddQ[e] := (y=Mod[x^2,n];Mod[x*PM[Mod[y,n],(e-1)/2,n],n])

ToyDecryptPM[ToyEncrypt[1000]]

Evaluating PM changes y: this can be seen below for 2^3 = 1 mod 7 and 2^4 = 2 mod 7

{{PM[2,3,7],y},{PM[2,4,7],y}}

This kind of side-effect can lead to nasty bugs. So we use ‘Module’ to make y a local variable to PM. (Mathematica has other scoping constructions, but in almost all cases, Module is the best one to use.)

Clear[y]; Clear[PM]

PM[x_,0,n_]:=1

PM[x_,e_,n_]:=Module[{y = Mod[x^2,n]},
If[EvenQ[e],PM[y,e/2,n],Mod[x*PM[y,(e-1)/2,n],n]]]

ToyDecryptPM[ToyEncrypt[1000]]

As expected, y is now unset

y

Map: Encrypting many numbers

ToyEncryptList[xs_]:=Map[ToyEncrypt,xs]

ToyDecryptList[ys_] :=Map[ToyDecrypt,ys]

ToyEncryptList[{1,2,3,20,21,22,30,31,32}]

ToyDecryptList[ToyEncryptList[{1,2,3,20,21,22,30,31,32}]]

Map can be replaced with /@. (You will see this, and the related @@ apply operator a lot in code written by experts.) For instance

ToyEncrypt/@{1,2,3,20,21,22,30,31,32}

Map is sometimes unnecessary in Mathematica: functions declared as Listable are mapped automatically. (This is useful, but also a potential source of confusion and bugs.)

{1,2}*3

ToyEncryptL=ToyEncrypt;SetAttributes[ToyEncryptL,Listable]

ToyDecryptL=ToyDecrypt;SetAttributes[ToyDecryptL,Listable]

ToyDecryptL[ToyEncryptL[{1,2,3,20,21,22,30,31,32}]]

Extension to strings: ToCharacterCode (a Listable function) converts a string to a list of numbers (using ASCII as it happens)

WordToNumbers[w_]:=ToCharacterCode[w]

NumbersToWord[xs_]:=FromCharacterCode[xs]

WordToNumbers["HELLO"]

NumbersToWord[{72,69,76,76,79}]

ToyEncryptWord[w_]:=Map[ToyEncrypt,WordToNumbers[w]]

ToyDecryptToWord[xs_]:=NumbersToWord[Map[ToyDecrypt,xs]]

ToyEncryptWord["HELLO FUNCTIONAL PROGRAMMERS"]

ToyDecryptToWord[ToyEncryptWord["HELLO FUNCTIONAL PROGRAMMERS"]]

Further Map/Reduce and recursion problems

True if list contains only odd numbers

AllOdd[xs_]:=Length[Select[xs,EvenQ]]==0

{AllOdd[{1,3,5}],AllOdd[{1,2,3}]}

CountList (implements the built in function Tally)

Tally[{1,1,2,2,1,3,1}]

The first pair in the output is {First[xs], number of elements of xs equal to First[xs]}. This is clearly a function of interest in its own right, so we start by implementing it.

FirstPair[xs_] := {First[xs],Length[Select[xs,Equal[First[xs],#]&]]}

Then define CountList on the empty list, and by recursion on a general list

CountList[{}]:={}

CountList[xs_]:=Join[{FirstPair[xs]},
CountList[Select[xs,Not[Equal[First[xs],#]]&]]]

CountList[{1,5,2,1,2,1}]

Note there is a nasty trap here: CountList doesn’t work as expected when the input contains lists rather than just numbers

CountList[{{3, 4}, 1, 2}]

{{{3, 4}, 1}}

CountList[{1,2,{3,4},{3,4},1}]

The problem is that Equal refuses to evaluate when the arguments are not obviously comparable. For example

Equal[{1,2},3]

To avoid this, replace Equal with SameQ

CountList[xs_]:=Join[{FirstPair[xs]},
CountList[Select[xs,Not[SameQ[First[xs],#]]&]]]

CountList[{1,5,2,1,2,1,{3,4},1,{3,4},2}]

Asymptotics of the sum of ϕ(k)/k for k between 1 and n

f[n_]:=Apply[Plus,Table[EulerPhi[k]/k,{k,1,n}]]

Sum is more natural here, but this definition shows how to do it using only Apply and Table

fAlt[n_]:=Sum[Table[EulerPhi[k]/k,{k,1,n}]]

TableForm[Table[{n,f[n],N[f[n]],N[f[n]/n]},{n,1,10}]]

TableForm[Table[{n,N[f[n]],N[f[n]/n]},{n,Select[Range[1,10000],Mod[#,1000]==0&]}]]

This suggests that the sum of ϕ(k)/k for k from 1 up to n is asymptotically cn, for some constant c about 0.608. In fact this is the case, with c = 6/π^2. (Search for ‘coprime probability’ to find out why.)

N[6/Pi^2]

Mergesort

First split the list into two

HalfLists[xs_]:={Take[xs,Floor[Length[xs]/2]],Drop[xs,Floor[Length[xs]/2]]}

Function to merge two sorted lists

MergeLists[{},zs_] := zs

MergeLists[ys_,{}] := ys

MergeLists[ys_,zs_]:=If[First[ys]≤ First[zs],
Join[{First[ys]},MergeLists[Drop[ys,1],zs]],
Join[{First[zs]},MergeLists[ys,Drop[zs,1]]]]

Now ready for Mergesort: the ‘base cases’ are lists of length 0 (arise as the left half list when splitting a list of length 1), and lists of length 1.

MergeSort[{}]:={}

MergeSort[{x_}]:={x}

MergeSort[xs_]:= MergeLists[
MergeSort[HalfLists[xs][[1]]],
MergeSort[HalfLists[xs][[2]]]
]

MergeSort[{4,3,5,6,7,3,2,1}]

Mergesort sorts a list of length n in time O(n log n). Further exercise: implement Quicksort in Mathematica

Comparing lists under lexicographic order (using patterns)

LexGreater[{},{}] := True

LexGreater[{},_] := False

LexGreater[_,{}]:=True

LexGreater[xs_,ys_] := Which[
First[xs]> First[ys],True,
First[xs]< First[ys],False,
First[xs]==First[ys],LexGreater[Drop[xs,1],Drop[ys,1]]]

LexGreater[{1,3,5,2},{1,3,4,4}]

Plotting the derivative of a function (hard for annoying reasons)

The obvious way to write this is

PlotDerivativeBad[f_] := Plot[D[f[x],x],{x,-3,3}]

This fails!

PlotDerivativeBad[Sin]

D[Sin[x],x]

The problem is that Plot[g[x],{x,-3,3}] evaluates g by substituting real numbers for g. Here we have g[x] := D[Sin[x],x], and so Mathematica substitutes real real numbers, such as -2.87743 for x, giving the meaningless D[Sin[-2.87743],-2.87743].

To fix this problem, we need D[Sin[x],x] to already have been evaluated, before it is used in Plot. So this is a time when we need Set (=) rather than SetDelayed (:=)

Clear[PlotDerivative]

PlotDerivative[f_]:=(g[x_] = D[f[x],x];Plot[g[x],{x,-3,3}])

PlotDerivative[Sin]

To avoid potentially changing the definition of g, it would be even better to write

PlotDerivativeSafe[f_]:=Module[{g},g[x_] = D[f[x],x];Plot[g[x],{x,-3,3}]]

PlotDerivativeSafe[Sin]

Another solution uses patterns to change the variable under Mathematica’s nose

PlotDerivativeAlt[f_] := Plot[D[f[x],x]/.{x->y},{y,-3,3}]

PlotDerivativeAlt[Sin]

Derangements: let d_n be the number of derangements of {1,...,n}

To split up IsDerangement, we first write a function that determines the fixed points of a permutation xs, using pattern matching to take those i such that xs[[i]] = xs

FixedPoints[xs_] := Map[First,Cases[Table[{i,xs[[i]]},{i,1,Length[xs]}],{i_,i_}]]

IsDerangement[xs_] := FixedPoints[xs]=={}

Derangements[n_] := Select[Permutations[Range[1,n]],IsDerangement]

For example

Derangements[4]

NumberOfDerangements[n_] := Length[Derangements[n]]

Table[{n,NumberOfDerangements[n]},{n,0,8}]

Table[{n,N[NumberOfDerangements[n]/Factorial[n]]},{n,0,8}]

ListPlot[Table[NumberOfDerangements[n]/Factorial[n],{n,1,8}]]

In fact d_n/n! converges to 1/e

N[Exp[-1]]

ComposePermutations[xs_,ys_] := Table[ys[[xs[[i]]]],{i,1,Length[xs]}]

ProductsOfDerangements[n_] :=
Module[{ds := Derangements[n]},
Flatten[Table[ComposePermutations[ds[[i]],ds[[j]]],{i,1,Length[ds]},{j,1,Length[ds]}],1]]

ProductsOfDerangements[4]

ProductProbability[n_] := Length[Select[ProductsOfDerangements[n],IsDerangement]]/
NumberOfDerangements[n]^2

Table[{n,ProductProbability[n]},{n,2,5}]

A more efficient program would pick one derangement in each fixed-point free conjugacy class of S_n.

All lists with specified entries: shows that nesting isn’t always what one wants

Table[i+j,{i,1,10},{j,1,5}]

So to run multiple iterators without creating a list of lists, need to join up lists

Join@@Table[i+j,{i,1,10},{j,1,5}]

Ls[_,0]:={{}}

Ls[xs_,n_]:=Join@@Table[{x}~Join~l,{x,xs},{l,Ls[xs,n-1]}]

Ls[{1,2,3},3]

This is built into Mathematica as tuples

Tuples[{1,2,3},3]