Math 640: EXPERIMENTAL MATHEMATICS (GAMES!)
Spring 2013 (Rutgers University) Webpage
http://sites.math.rutgers.edu/~zeilberg/math640_13.html
Last Update: May 21, 2013
- Teacher:
Dr. Doron ZEILBERGER ("Dr. Z")
-
Classroom:
Allison Road Classroom Building
[Busch Campus], IML Room 116
-
Time: Mondays and Thursdays , period 3 (12:00noon-1:20pm)
-
"Textbooks": Links (see below).
-
Dr. Zeilberger's Office: Hill Center 704
-
Dr. Zeilberger's E-mail:
zeilberg at math dot rutgers dot edu (you MUST have MathIsFun in the message,
or ExpMathRocks)
-
Dr. Zeilberger's Office Hours (Spring 2013): MTh 10:30am-11:30am and by appointment
Description
Experimental Mathematics used to be considered an oxymoron,
but the future of mathematics is in that direction.
In addition to learning the
philosophy and methodology of this budding field, students will
become computer-algebra wizards, and that should be very helpful in whatever
mathematical specialty they are doing (or will do) research in.
We will first learn Maple, and how to program in it.
This semester we will focus on Games (of all kinds!), classical (the subject of so-called Game Theory
for which John Nash (and quite a few other people) got a Nobel), Combinatorial
(from Nim to Chess, via Backgammon), and Gambling (what Wall Street does).
In addition to the actual, very important content, students will
master the methodology of computer-generated and
computer-assisted research that is so crucial for their future.
There are no prerequisites, and no previous programming knowledge is assumed.
Also, no overlap with previous years.
The final projects will with high probability lead to publishable articles,
and with strictly positive probability, to seminal ones.
"Textbooks"
Added March 18, 2013: Pick a
final project
How to submit your homework
If you don't already have a public_html directory, go
to command-line Unix (or Linux), and type
-
mkdir public_html
-
chmod 755 public_html
-
cd public_html
-
mkdir em13
-
chmod 755 em13
-
cd em13
-
create the homework file, call it hw1.txt (and in the future, hw2.txt etc.)
and do
-
chmod 755 *
If you are a math grad student, and using your math account, you don't
have to Email me, but if you are using eden or another computer,
please Email me the url of your homework directory.
Diary and Homework
Preliminary homework assigned Jan. 21, 2013 (Due Thurs. Jan. 24, 2013)
-
Carefully read and understand
Dr. Z.'s two-page introduction to non-partizan combinatorial games,
and do all the exercises (by hand!)
Programs done on Thurs., Jan. 24, 2013
-
C1.txt ,
Contains procedures
-
f2(n): inputs an integer n, and outputs 0 or 1 according
to whether n pennies is a losing or winning position, repectively,
in a penny-removing game, where a legal move
consists of removing either one or two pennies, and
person who can't move (because there are no pennies left) is
the ultimate loser.
-
fk(k,n): inputs a positive integer k, and an
integer n, and outputs 0 or 1 according
to whether n pennies is a losing or winning position, repectively,
in a penny-removing game, where a legal move
consists of removing any positive number of pennies ≤k, and
person who can't move (because there are no pennies left) is
the ultimate loser.
-
fS(S,n): inputs a set of positive integers S, and an
integer n, and outputs 0 or 1 according
to whether n pennies is a losing or winning position, repectively,
in a penny-removing game, where a legal move
consists of removing a of pennies belonging to S, and
the person who can't move (because there are no pennies left) is
the ultimate loser. It is OK to leave the pile with a negative
number of pennies (but see homework below).
[Here is Edinah Gnang's
Sage version]
Homework for Thurs., Jan. 24, class (due Sunday Jan. 27, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw1.txt
-
(For novices only)
Read and do all the examples, plus make up similar ones,
in the first 30 pages of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Note: a few commands are no longer valid in today's Maple. The
most important is that " has been replaced by %.
-
(For novices only) Write a procedure S(n) that inputs
a positive integer n, and outputs the n-th Somos number,
defined by the recurrence
S(n)*S(n+4)=S(n+1)*S(n+3)+S(n+2)^2
with S(1)=1,S(2)=1,S(3)=1, S(4)=1.
What is S(100)?
-
(Mandatory for experts, optional for novices)
Write a procedure, FindUPi(L,i), that inputs
a list of numbers, and a positive integer i, and
finds out whether there exists a list M of length i, such that
L ends with at least three repeated copies of M, or else
returns FAIL.
For example
FindUPi([3,4,2,4,5,2,4,5,2,4,5],3);
should return [2,4,5], and
FindUPi([3,4,2,4,5,2,4,5,2,4,5],2);
should return FAIL.
-
(Mandatory for experts, optional for novices)
Write a procedure, FindUP(L), that inputs
a list of numbers, and
finds out whether there exists a list M, such that
L ends with at least three repeated copies of M, or else
returns FAIL.
For example
FindUP([3,4,2,4,5,2,4,5,2,4,5]);
should return [2,4,5], and
FindUP([1,2,3,4,5,4,5]);
should return FAIL.
-
[Optional for everyone, I don't know the answer]
Experiment with various S, and discover whether
[seq(fS(S,n), n=1..1000)]
is ultimately periodic. How large are the periods?
Can you predict from the size, or other features of S,
the length of the period?
Added Jan. 28, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Monday, Jan. 28, 2013
Contains the following procedures:
-
IsCyclic(G): inputs a directed graph (described as
a list of sets, andoutputs true (false) if it has
a cycle
-
fGi(G,i): inputs a directed graph G and
an integer i between 1 and nops(G) and outputs
0(resp. 1) if vertex i is a Losing (Winning) position
-
AM(G): adjacency matrix of G (written a matrix M)
-
WL(G): Winning list.
Inputs a directed G (representing a game)
and outputs a list of length nops(G) such that
L[i]=0 or 1 according to whether position (vertex) i
is a losing or winning position
-
RandGame(n,k): inputs pos. integers n and k
and outputs a random digraph that has the
property that out of vertex i all the labels
are less than i, and there are (usually) k outgoing
edges
[Added Jan. 29, 2013: Philip Benjamin proved that every acyclic digraph can be labeled to have
the above property, so this property is wlog. See his
Philip Benjamin's nice writeup]
[Here is Edinah Gnang's
Sage version]
Homework for Monday, Jan. 28, 2013 class (due Sunday Feb. 3, 2013
10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw2.txt.
[Note: please do them by Thurs., if possible.
You would get more out of the class if you do the homework as soon
as possible.]
-
(For novices only)
Read and do all the examples, plus make up similar ones,
in the pages 30-60 of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Note: a few commands are no longer valid in today's Maple. The
most important is that " has been replaced by %.
-
(Mandaotry for novices, optional for experts)
Write a program RanMat(n,k) that inputs a positive integer
n and a pos. integer k, and outputs a random matrix whose entries
are non-negative integers less than k.
-
(Mandaotry for novices, optional for experts)
Using RanMat(n,k), Write a program HowManySigulars(n,k,K)
that inputs n and k as above
and a large positive integer K, and picks
K randomly chosen matrices and finds out what fraction
have a determinant 0 mod k.
What did you get for HowManySigulars(3,5,10000)?
Run it several times and see whether your answers are close.
-
(Mandaotry for experts, optional but strongly recommended for novices)
Using procedures RandGame(n,k) and WL(G)
we did in
class,
write a program, AveLosing(n,k,K) that inputs also a large
positive integer K, and finds out the average number of losing
positions in K randomly chosen games (using RandGame(n,k) K times).
If K is big, and you run it many times, do you get similar answers?
Can you guess what the answer should be?
-
(Mandaotry for experts, optional for novices)
Write a program ProbWinStupid(G,i) that inputs a game given
as a directed graph, and a vertex i between 1 and nops(G)
and computes the probability of winning if both players
play completely randomly, by picking (if possible)
uniformly at random one of the legal available move.
Of course, if the vertex is a sink then the probability is
0, since
-
(Challenge, I don't know the answer. 5 dollars to be shared)
Try to conjecture closed form (if possible) expressions
for the above probability, call it W(n,k),
of winning, with the above random stupid play,
for the penny removing game where right now you have
n pennies and you can remove from 1 to k pennies at each turn.
[Clarification added Feb. 1, 2013: Matthew Russell wisely pointed out that there are two interprations
when the number of pennins left, i, is less than k. I meant that must remove up to min(i,k) pennies,
but you are welcome to do it where you are allwed to have the pile left with a negative number
of pennies, and then your opponet lost]
Added Feb. 4, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thurs., Jan. 31, 2013
Contains the following procedures:
-
W(G,i): same as fGi(G,i) from C2.txt,
but done more elegantly using an idea of
Matthew Russell (inspired by Pat Devlin
and Tim Naumovitz)
-
WL(L,R,i): inputs two digraphs L and R
(with nops(L)=nops(R), of course),
and a positive integer i, representing the positions,
where L[i] (R[i]) is the set of positions reachable
in ONE move by Left (Right).
Outputs 0(1) if it is a Losing (Winning) position for Left.
-
WR(L,R,i): inputs two digraphs L and R
(with nops(L)=nops(R), of course),
and a positive integer i, representing the positions,
where L[i] (R[i]) is the set of positions reachable
in ONE move by Left (Right).
Outputs 0(1) if it is a Losing (Winning) position for Right.
-
Nim1(a): inputs a non-neg.
integers, representing a Nim-1 position with
a counters, outputs whether it
is losing (0) or winning (1).
A legal move consist of
removing any number of pennies (up to the total)
from that pile
-
Nim2(Li): same as Nim1 but with two piles. Li
is a list of non-neg. integers of length 2.
-
Nim3(Li): same as Nim3 but with three piles. Li
is a list of non-neg. integers of length 3.
-
Nim(Li): inputs a list of any length of non-neg. integers
representing sizes of Nim piles (with nops(Li) piles), and
outputs 0(1) if it is a losing (winning) position.
[Here is Edinah Gnang's
Sage version]
Homework for Thursday, Jan. 31, 2013 class (due Sunday Feb. 3, 2013
10:00pm)
-
(For novices only)
Read and do all the examples, plus make up similar ones,
in the pages 60-90 of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Note: a few commands are no longer valid in today's Maple. The
most important is that " has been replaced by %.
-
(For everyone)
Using procedure Nim(Li) (Nim3(Li)) write a procedure
F(k,N) that inputs positive integers k and N and outputs
all the triples [k,a,b] such that [k,a,b] is a LOSING
position with 0<=a<=b<=N. Can you guess the set
F(k,infinity) for k=0,1,2,3,4, .. up to k=10.
-
(For everyone)
The Wythoff game is like 2-pile Nim with the extra
move that one can remove the SAME number of pennies
from BOTH piles. Write a procedure Wyt(a,b) that
inputs 0(1) according to whether [a,b] is a losing
(winning) position
-
Prove that for every non-negative integer n, there is exacty one Wythoff Losing
position [a,b] with b=a+n. Let's call it An so
[A0, A0], [A1, A1+1] , ..., [A2, A2+1] , ...
are the losing positions (with a<=b).
Using Wyt(a,b) write a procedure SeqA(n) that inputs a pos. integer n and
outputs the sequence
[A0,A1,A2, ..., An]
Is that sequence in Sloane?
Conjecture a value for the limit of An/n as n goes to infinity
-
(For experts only) Write a procedure kWyt(Li)
Generalizing Wythoff's game to a k-pile Wythoff game which is like k-pile Nim
and in addition one can take the same (positive) number of pennies from any TWO out of the k piles.
What is kWyt([3,4,5,6])?
Added Feb. 4, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Monday, Feb. 4, 2013
Contains, in addition to the previous procedures
the following procedures (listed in Help4();)
-
mex(S): inputs a set of non-neg. integers and outputs
the minimum of the set {0,1,2,3,...}\S
-
SG(G,i): inputs a digraph (representing a cycle-less
combinatorial game with positions 1, 2, ..., nops(G)
and G[i] is the subset of {1, ..., i-1} reachable
from position i in one legal move.
Outputs the value of the Sprague-Grundy function
at vertex i
-
SG2N(i,j): The S-G function of positon [i,j] in 2-pile Nim,
using the defining recurrence.
-
NimSum(L): the Nim sum of the list of non-neg. integers L,
doing it the fast way, using the theorem that
the Nim sum of a set of numbers can be computed
by converting them to binary and adding them
without carry. It uses the Xor function in the Maple
Bits package. This nice approach was suggested by
Philip Benjamin.
-
NimSumPat(L): Pat Devlin's rendition of NimSum(L)
-
SG1(m,n): Like SG2N(m,n), but done using the fast recurrene
SG(2*m,2*n)=2*SG(m,n)
SG(2*m,2*n+1)=2*SG(m,n)+1
SG(2*m+1,2*n)=2*SG(m,n)+1
SG(2*m+1,2*n+1)=2*SG(m,n):
Homework for Monday, Feb. 4, 2013 class (due Sunday Feb. 10, 2013
10:00pm)
-
(For novices only)
Read and do all the examples, plus make up similar ones,
in pages 91 to the very end of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Note: a few commands are no longer valid in today's Maple. The
most important is that " has been replaced by %.
-
[for everyone, purely human problem]
Starting from the recursive definition of the Sprague-Grundy
function for two-pile Nim (let's call it f(m,n))
f(m,n)=mex({f(m-1,n),f(m-2,n), ..., f(0,m),
f(m,n-1),f(m,n-2), ..., f(m,0)} ) ,
prove, by induction on m and n, the recurrences
f(2m,2n)=2f(m,n)
, f(2m+1,2n)=2f(m,n)+1
, f(2m,2n+1)=2f(m,n)+1
, f(2m+1,2n+1)=2f(m,n)
(that are equivalent to the fact that f(m,n) is the sum-without-carry
of m and n in binary)
-
[for everyone, purely human problem]
Prove that a position is Losing iff its Sprague-Grundy function is 0
-
[For everyone]
Write a program SGwyt(i,j) that inputs
two non-negative integers i and j and outputs the
Sprague-Grundy function of position (i,j) in Wythoff's
game, where a legal move consists of taking a positive number
of pennies from either pile, or the Same number from both.
-
Carefully read and understand
Chs. 7 and 8 of ONAG (alias the BIBLE)
[Added Feb. 7, 2013: Ch. 8 is hereby made optional. It
is more important to read Ch. 1]
-
[Mandatory for experts, optional to novices]
Write a program, WinningMove(G,i),
that inputs a game given as a digraph G, and a position i
(between 1 and nops(G)),
and outputs a winning move (a member of G[i])
if one exists, or FAIL if i
is a losing (alias P) position.
-
[Optional, open problem, as far as I know, 10 dollars]
The losing positions of Wythoff's game, alias the position with
Sprague-Grundy function value 0,
are well-known to be
[trunc(phi*n),trunc(phi*n)+n]
(and of course, by symmetry [trunc(phi*n)+n,trunc(phi*n)] )
where phi is the Golden ratio. Can you
find an analogous characterization for those that have
S-G function value 1?
-
[optional, famous open problem, 50 dollars]
Find a fast algoritm to compute SGwyt(i,j) analogous to
Nim-Sum. What is
SGwyt(13783204640276,6351234024524) ?
Added Feb. 11, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thursday, Feb. 7, 2013
Contains, in addition to the previous procedures
the following procedures:
-
Children(L,n): inputs a legal placement of
queens on a nops(L) by n board and finds
all the legal extensions
For example,
Children([2,5,1],6);
should give
{[2,5,1,4],[2,5,1,6]}
-
Qnk(n,k): the set of ways of placing k-non-attacking
queens on a k by n chess-board, represented as
[a1,a2,,,,, a_k] where
a1 is the column of the queen at the first row
a1 is the column of the queen at the second row etc.
-
RandDi(n,k): a random directed graph with n
vertices and each vertex has (usually) outdegree k
-
Sphere(G,i,k): inputs a directed graph G
(represented as a list of sets where G[i]
is the set of outgoing neighbors of vertex i)
and a vertex i (an integer between 1 and nops(G))
and a non-neg. integer k, outputs the
set of all vertices for which you can walk
from i to it with k steps but no less
-
Spheres(G,i): inputs a directed graph
and a vertex i (from 1 to nops(G))
and outputs the list of all concentric spheres
until it is empty
Homework for Thursday, Feb. 7, 2013 class (due Sunday Feb. 10, 2013
10:00pm)
-
Carefully read and understand
Ch. 1 of ONAG (alias the BIBLE)
(Note: Chapter 8 is not as important right now, so even though I
asked you to read it for Monday's homework, it is hereby made optional)
-
By modifying procedure Sphere(G,i,k) of
C5.txt ,
write a procedure SphereWithPaths(G,i,k) that outputs
a set of paths of minimal length k starting at i,
where a path is given a list of vertices,
but in case of multiple paths of length k ending at a given
vertex, you only have one of them.
For example
SphereWithPaths([{2,3},{4},{4},{}],1,2) ;
may output
{[1,2,4]}
or
{[1,3,4]}
(which one being a random choice of the machine)
-
a krook is a hybrid of a chess King and a chess rook,
i.e. it is like a Queen that can only move diagonally one unit.
Modify program Qnk(n,k) to write a program, Knk(n,k)
to generate the set of
ways of placing k non-attacking krooks on a k by n chess-board.
is the (start of the) sequence
[seq(nops(Knk(n,n),n=1..infinity)]
in Sloane?
-
a mook of order r
is a hybrid of a chess King and a chess rook,
i.e. it is like a Queen that can only move diagonally up to r units.
In particular a mook of order 1 is a krook.
Modify program Knk(n,k) to write a program, Mnk(n,k,r)
to generate the set of
ways of placing k non-attacking mooks of order r on a k by n chess-board.
(Of course Mnk(n,k,1) is the same as Knk(n,k) and Mnk(n,k,k) is the
same as Qnk(n,k), please check!).
For what values of r are the (starts of the) sequences
[seq(nops(Mnk(n,n,r),n=1..infinity)]
in Sloane?
-
[for experts, optional challenge for novices]
Use directed graphs to solve the problem of
n missionaries and n cannibals having to cross
a river with a row boat that can only contain
one or two persons, where no missionary gets eaten up.
(so at no time can either bank have more cannibals than
missionaries unless the number of missionaries in that
bank is zero)
-
Consider
peg-solitaire but played on a k by n board,
where a position is represented as a 0-1 matrix,
and we use the data structure of lists-of lists
(a list of length k where each entry is a list of
length n, so M[i][j]=1 iff the spot at the i-th
row and j-th column has a peg). Write a program
PegSoliChildren(M) that inputs such a position and
outputs the set of all the positions (with the number
of 1's reduced by 1) that are legally reachable from it.
-
[small challenge]
Is the four-by-four Peg Solitaire with a one hole at a corner
solvable?
-
[Big challenge, 5 dollars to be divided]
Extend this methodology (hopefully not
runing out of memory)
to solve the usual (English) Peg-Solitaire, with the usual board.
Added Feb. 11, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Mandatory homework due Feb. 14, 2013 (bring to class!, on paper!)
[Added Feb. 12, 2013: there was a typo in the previous version
(implicitplot was misspelled, sorry!]
-
Find the unique polynomial P, of the variables x and y,
that satisfy the following boundary value problem PDE:
Solve the PDE (Hint P is a polynomial in x and y of degree 6 in each variable, and total degree 6).
diff(P,x$2)+diff(P,y$2)=
36*x^4+72*x^2*y^2-48*x^2+36*y^4-48*y^2+12-2*y^3-6*x^2*y
subject to the boundary conditions on the unit square:
-
P(0,y)=(y-1)^3*(y+1)^3,
-
P(1,y)=y^3*(y-1)*(y^2+y+1)
-
P(x,0)=(x-1)^3*(x+1)^3
-
P(x,1)=x^2*(x-1)*(x+1)*(x^2+1)
-
Go to xmaple, and type
plots[implicitplot](P,x=-4..4,y=-4..4,axes=none, grid=[1000,1000]);
-
print it out, cut the shape neatly with a pair of scissors,
and write inside it
I LOVE EXPERIMENTAL MATHEMATICS!
-
Hand it in at the Feb. 14, 2013 class
[Another Hint: a generic polynomial of total degree 6 is add(add(a[i,j]*x^i*y^j,i=0..6-j),j=0..6)
in the set of unknowns {seq(seq(a[i,j],i=0..6-j),j=0..6)}. Let Maple set-up a system
of linear equations that incorporate the above conditions, and solve for the a[i,j]'s]
Programs done on Monday, Feb. 11, 2013
Contains the following procedures:
-
GamesD1(k): all the games created in day<=k
-
GamesD(k): all the games created in day k
-
RandGame(k,m): a (most probably not uniformaly at)
random game created in day <=k, where G=[L,R]
and where L and R are always of size <=m
-
Age(G):the age of G=[{L},{R}]
-
CanLeftWin(G): Can Left win the game G?
-
CanRightWin(G): Can Right win the game G?
-
WhatKind(G): Is the game, positive, negative, zero, or fuzzy?
Additional Homework for Monday Feb. 11, 2013 class (due Sunday Feb. 17, 2013
10:00pm)
Create a file hw6.txt and place in your em13 directory.
-
[Modified Feb. 14, 2013, thanks to Phil Benjamin, who won a dollar]
Write a program ConWay1(L,R,i) that inputs two
acyclic directed graphs on {1, ...n} (where n=nops(L)=nops(R))
representing a partizan game, where their union is also
acyclic ,
and an integer i between 1 and n, representing a game-position
(alias "game" in Conway's sense) and outputs it in Conway's format,
as a pair of sets of simpler "games" where the very bottoms are empty sets.
-
Using ConWay1(L,R,i) above write a one-line program, ConWay(L,R), that inputs
two directed graphs in the above format and outputs the list of length n:=nops(L) (=nops(R))
whose i-th entry is ConWay1(L,R)
-
Define the signature of a partizan game given as a pair of directed graphs (L,R) the expression
a*"Zero"+b*"Positive"+c*"Negative"+ d*"Fuzzy"
where a+b+c+d=n, and there are a "Zero" positions, b "Positive" positions etc.
Write a program Sig(L,R) that inputs a pair of directed graphs with the same
number of vertices and output its signature.
-
[modified 2/14/13 thanks to a remark of Phil Benjamin]
By using RandGame(n,k) in
C2.txt ,
2*K times, write a program
AveSig(n,k,K), that picks K random games and outputs their average signature.
What do you get for AveSig(10,3,100)? Do you get similar results all the time?
Added Feb. 18, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thursday, Feb. 14, 2013
Contains the following procedures:
-
GE(x,y): inputs two (genuine) (surreal) numbers
in John Horton Conway's sense and outputs true
iff x ≥ y
-
LE(x,y) is x ≤ y?
-
GT(x,y): x > y ?
-
LT(x,y): x < y ?
-
IsNumber(x): Is x a legal (surreal) number?
-
ADD(x,y): Inputs two numbers (or games) and outputs the number
(or game) x+y
-
SN(N): the surreal rendition of the pos. integer N
-
Minus(x): inputs a number (or game) x, and outputs -x
Homework for Thursday Feb. 14, 2013 class (due Sunday Feb. 17, 2013
10:00pm)
Create a file hw7.txt and place it in your em13 directory.
-
Recall that multiplication by an integer is repeated addition.
By using ADD, write a procedure
NaiveMUL(a,n)
that inputs a surreal number a, a usual pos. integer n,
and outputs a+a+ ...+a (repeated n times).
-
Write a procedure CheckNaiveMul(K) that checks that
NaiveMUL(SN(m),n) equals (in Conway's sense) SN(m*n) for
all m,n from 1 to K. Are they equal in the usual sense?
(say for K=10)
-
By reading pages 7-9 of
ONAG
figure out a (recursive) definition of
(1/2)n for integer n. Write a procedure
HalfToPower(n)
that inputs a positive integer n and outputs the surreal representation
of (1/2)n.
-
By comining HalfToPower(n) and NaiveMUL(a,n) check that
2^n * (1/2)^n =1
-
Write a procedure MUL(x,y) that implements the
definition of xy on page 5 of
ONAG.
Check that
MUL(SN(2),SN(2))
equals (in Conway's sense) SN(4). How long did it take your computer?
Check that
MUL(SN(3),SN(4))
equals (in Conway's sense) SN(12). How long did it take your computer?
Added Feb. 18, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Monday, Feb. 18, 2013
Note:
Sorry for messing up. The file C8.txt contains two extra procedures
that we did not do in class. RollLD (for Rolling a loaded die),
and RandGame1(n), a random "thinkless" game of size n.
Please study them carefully
Contains the following procedures:
-
Games1(n): the set of games of size n where always
always there is at most one option:
-
P1(n,x):
Inputs a positive integer n, and a letter x, and output
a linear combination of x["Zero"], x["Fuzzy"],
x["Postive"], x["Negative"],
whose coefficients indicat how many games of the
correponding kind are there of size n.
-
NuGames1(n): the number of games of size n where always
there is at most one option:
-
RollLD(L): inputs a list of positive integers, L,
rolls a nops(L)-faced loaded die, getting a positive integer
i between 1 and nops(L) (inclusive)
where the probability of getting i is L[i]/add(L[i],i=1..nops(L)):
For example, to roll a usual fair die do
RollLD([1,1,1,1,1,1]);
-
RandGame1(n): a uniformaly-at-random thinkless game of size n
Homework for Monday Feb. 18, 2013 class (due Sunday Feb. 24, 2013
10:00pm)
Create a file hw8.txt and place it in your em13 directory.
-
Study carefully the simplified NuGames1(n) procedure in the
revised version of
C8.txt ,
and convince yourself that it is OK.
-
Check empirically procedure RollLD(L), by doing
add(x[RollLD([1,1,1,1,1,1])],i=1..6000);
add(x[RollLD([1,2,3])],i=1..6000);
and quite a few other L.
-
Test the procedure
RandGame1(n):
(that I added after class), by doing
add(x[RandGame1(10)],i=1..NuGames1(10)*100);
and see if you get roughly 100 for each of the thinkless games
of size 10.
-
[Challenge]
A thinkless game of size n is either [{},{g1}] for some
thinkless game g1 of size n-1, or [{g1},{}],
or has the form
[{g1},{g2}]
for some thinkless games g1 and g2 with size(g1)+size(g2)=n-1.
That's how we constructed recursively Games1(n).
Define a level-2 game a game where at any instance both
players have either no option, exactly one option, or exactly
two options.
Write a (recursive) procedure, G2(n), that inputs a positive integer n
and outputs the set of all level-2 games of size n.
Hint: the format of a level-2 game is one of the following
-
[{},{}] (only if n=1)
-
[{},{g1}], [{},{g1,g2}] (g1 NOT g2) , [{g1,g2},{}] (g1 NOT g2) , [{g1},{g2}],
-
[{g1},{g2,g3}] (g2 NOT g3) , [{g2,g3},{g1}] (g2 NOT g3),
-
[{g1,g2},{g3,g4}] (g1 NOT g2 and g3 NOT g4)
where g1,g2, etc. are level-2 games whose sizes add up to n-1.
[Added Feb. 19, 2013: I thank Matthew Russell for correcting a previous error (origically I said "add up to n",
of course we have to give credit to the [] in the original game]
-
[Human Challenge, 5 dollars] We saw in class that the sequence
"number of zero thinkless games of size n" is in Sloane, but with
a different meaning. Prove that they are the same for all n.
-
[Human Challenge, 5 dollars] We saw in class that the sequence
"number of fuzzy thinkless games of size n" is in Sloane, but with
a different meaning. Prove that they are the same for all n.
Added Feb. 25, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thurs., Feb. 21, 2013
-
C9a.txt
[Finishing up Conway-style combinatorial games, with procedure G(n,r) debugged
by Nathaniel Shar who won 20 dollars to be donated in his honor to the OEIS foundation]
-
C9.txt
[Code written by Philip Benjamin],
[Added Feb. 22, 2013: here is
Nathaniel Shar's version ]
Today we did the simplest possible (not completely trivial) analog of the Bearoff stage
of Backgammon with two-faced ONE die, and a two-position home. As follows
-
"Bacterial" Backgammon
-
Two squares on board 1,2
-
Position: a,b # of counters on each square (zero or more)
-
Roll die: 1 or 2
-
Legal move:
-
Bearoff (google this!)
-
Input [a,b],i
-
i=1 --> [a,b] -> [a-1,b+1] or [a,b-1] (if possible)
-
i=2 --> [a,b] -> [a-1,b] or [0,b-1] otherwise (must do in order)
-
Win if [0,0]
-
LegalMoves(a,b,i): inputs non-negative integers a and b and
an integer i in {1,2} and outputs the set of all legal
positions reachable in one move in a bearoff mini Backgammon
where the square distance 2 units away has a counters
and the the square distance 1 units away has b counters
and the 2-faced die landed on i
-
f(a,b,i): inputs non-negative integers a and b and an integer
i in {1,2} and outputs the expected number of moves until finishing
using optimal play (minimizing the expected duration).
-
F(a,b): inputs non-negative integers a and b
outputs the expected number of moves (before rolling the die)
until finishing
mini-Bearoff Backgammon with 2 squares and a 2-faced
fair di, using optimal play (minimizing the expected duration).
-
BestMove(a,b,i): inputs non-negative integers a and b
and an integer i in {1,2} and outputs the pair [duration, BestMove]
Homework for Thursday Feb. 21, 2013 class (due Sunday Feb. 24, 2013
10:00pm)
Create a file hw9.txt and place it in your em13 directory.
-
Modify f(a,b,i) and F(a,b), calling them fS(a,b,i) and FS(a,b) respectively, where the goal is
to prolong the game as much as possible. Does the limit of FS(n,0)/F(n,0), as n goes to infinity
exist? It is does, what is it?
-
Find empirically constants c1 and c2 such that F(n,0) is asymptotically (and very close to!) c1*n+c2.
Can you explain why c1 is what it is?
-
Modify f(a,b,i) and F(a,b), calling them fR(a,b,i) and FR(a,b) respectively, where the strategy is
to pick a move at random (equally likely) whenver there is a choice of more than one legal moves.
Does the limit of FR(n,0)/F(n,0), as n goes to infinity exist? If it does, what is is?
-
Modify f(a,b,i) and F(a,b), calling them fG(a,b,i) and FG(a,b) respectively, where the strategy is
the "greedy one", trying to get out the pieces if necessary, so one would never
go from [a,b] to [a-1,b+1] if b>0 but would go to [a,b-1] thereby decreasing the number of
counters on the board.
Does the limit of FG(n,0)/F(n,0), as n goes to infinity exist? If it does, what is is?
-
[Challenge]
Modify f(a,b,i) and F(a,b) to play a bearoff game,
calling them fD(a,b,i,j) ([i,j]=[1,1],[1,2],[2,2]) and FD(a,b),
but with TWO dice (still fair two-faced ones)
with the backgammon convention that whenever there is a double you can go four times that number. In other words,
if if lended double 1, one has four "1" moves, and if it landed "double 2" then you have four "2" moves.
fD(a,b,i,j) is the expected optimal duration if the dice landed [i,j].
What can you say about the asymptotics of FD(n,0)?
-
[Human Challenge. 5 dollar donation to OEIS in your honor if you get it]
If using C9a.txt
you type:
read `C9a.txt`:
seq(nops(G(n,n)),n=1..9);
you would get, as output
1, 2, 5, 18, 66, 266, 1111, 4792, 21124
(using procedure G(n,r), kindly debugged by Nathaniel Shar) that tells you the total
number of Conway-style games of size n. It happens (conjecturally for now) to be the same as
A005753. Prove the equality rigorously.
-
[Human Challenge. 10 dollar donation to OEIS in your honor if you get it]
Find interpretations in terms of families of rooted trees, for the number of Conway-style games of size n that are
(a) zero-games (b) fuzzy games (c) positive games.
Added Feb. 25, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Monday, Feb. 25, 2013
Today we did general Bearoff Backgammon with an r-faced (fair) die
labeled(1, ...r ) and a field of size nops(L)
starting with L=[L[1], ..., L[nops(L)]] counters
where the L[i] counters at the ith-location are distance i from the end
Today's file is:
that contains the following procedures
-
LegalMoves(L,i): inputs a list of length nops(L)
and a pos. integer i telling you that you have the option
to move any counter i pieces forward, if possible
without waste. And if you must waste,
you must minimize it
-
F(L,r): inputs a list L of non-negative integers L and
a pos. integer r, and outputs the Expected number of
rolls until the end under OPTIMAL play (minimizing
the expected duration)
-
f(L,r,i): inputs L and r as above and i between 1 and r
and outputs the OPTIMAL move, followed by the
expected length if you do it
-
SimulateOptimalGameV(L,r):
simulates a Bearoff game with initial position L
and a fair r-faced die (labeled 1, ...r)
using pseudo-random number:Verbose version.
Outputs a story and the number of rolls.
Terse version. Only outputs the number of rolls.
Homework for Monday Feb. 25, 2013 class (due Sunday March 3, 2013
10:00pm)
Create a file hw10.txt and place it in your em13 directory.
-
Let's define the greedy strategy in Bearoff to be taking a piece out
of the board if possible, and if not, move with a piece closest
to the end (e.g. if the position is [0,0,0,1,0,1] and you rolled a 1,
you should go to [0,0,1,0,0,1] and not to [0,0,0,1,1,0]).
Modify all the procedures we did today, by naming them the same,
but with a G at the end (e.g. FG(L,r) instead of F(L,r) etc.)
-
Modify the procedures we did today, calling them the same names but
with a P at the end (for Pessimal, e.g. FP(L,r) instead of F(L,r) ),
that tries to maximize the
expected duration of the Bearoff game.
-
Modify the procedures we did today, calling them the same names but
with an R at the end (for Random e.g. FR(L,r)), that
picks a move uniformly at random among all the members of
LegalMoves(L,i) (of course if it only has one member, you must
do that move).
-
Compare the sequences
-
[seq(F([0$(r-1),n],r), n=1..10)]
-
[seq(FG([0$(r-1),n],r), n=1..10)]
-
[seq(FP([0$(r-1),n],r), n=1..10)]
-
[seq(FR([0$(r-1),n],r), n=1..10)]
For r=2, ...,7. (If Maple allows).
-
[Challenge, I don't know the answer]
conjecture linear expressions of the form c1(r)*n+c2(r) for
r between 2 and 7, that approximates
F([0$(r-1),n],r) for large n.
-
Learn how to play
Backgammon,
and practice by playing on-line.
-
Start reading the
few pages from Ch.1
and
Ch.2
from Karl Sigmund's
masterpiece
to be ready to study Game Theory, that we will probably start after Spring break.
Added March 4, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thurs., Feb. 28, 2013
Today's files are:
C11.txt contains:
-
LegalMoves(L,i): inputs a list of length nops(L)
and a pos. integer i telling you that the option
to move any counter i pieces towards the end
without waste. And if you must waste, minimize it
[taken from C10.txt]:
-
fW(L1,L2,r,i): The probability of L1 winning
if the roll was i, followed by the optimal move
that maximizes your probability of winning
-
FW(L1,L2,r): Probability of winning if your are about
to roll and your position is L1 and your opponent's
is L2 in a SINGLE r-faced fair die
-
[Added after class]
Exceptions3(N,r): Inputs a positive integer N and a pos. integer r in {1,2,3}, and
outputs all pairs of positions [[i1,i2,i3], [j1,j2,j3]] in Bearoff with house-size 3
and 3-faced fair ONE die, and i1+i2+i3=N, j1+j2+j3=N,
where the two strategies disagree, and by how much they lose
from their perspectives.
C11a.txt contains:
-
IsBearOff(POS,b,h):
are all the pieces of the player home yet?
LegalMoves(POS,b,h,i): all the legal moves
for player POS[1] (White)
with one-die roll i from position POS with
board-size b and home-size h
[Just started!]
Homework for Thurs. Feb. 28, 2013 class (due Sunday March 3, 2013
10:00pm)
Create a file hw11.txt and place it in your em13 directory.
-
Study the new procedure Exceptions3(N,r) carefully, and find out the smallest N for which
Exceptions3(N,1) is non-empty, and the smallest N (less than 15, if it exists) where
Exceptions3(N,2) is non-empty
-
Write the (easier analog), Exceptions2(N,r) for a house-size of 2 squares, and a 2-faced die.
What is the smallest N less than 20 for which Exceptions2(N,1) is non-empty?
-
Write the (harder analog), Exceptions4(N,r) for a house-size of 4 squares, and a 4-faced die.
What is the smallest N less than 7 for which Exceptions4(N,1) is non-empty?
-
Write a procedure PWG(L1,L2,r) that inputs lists L1, L2 denoting the positions of White and Black, and
a positive integer r indicating the die-size, and outputs the probability of White (who is about to move)
winning if both players follow the greedy strategy.
-
[Big Challenge, $20 donation to the OEIS foundation in honor of the people who will do it, it is OK
to work in teams, and share the burden, and divide up the task]
Finish procedure LegalMoves(POS,b,h,i) that we just started today for generalized Backgammon with
board-size b, house-size h, where the roll (of one die) was i.
Added March 4, 2013: See Patrick Devlin's brilliant
program.
Added March 4, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Added March 3, 2013: Here is
C10C11.txt ,
combining C10.txt and C11.txt and modifying Exceptions3(N,i) to only
include stictly positive losses of probability
(bug pointed out by Nathanaiel Shar)
Programs done on Monday, March 4, 2013
Today we started "gambling trees", and we learned how
to generate, uniformly at random a full binary tree
with a given number of leaves.
Today's file is:
It contains the following procedures:
-
BT(n): The set of all fulll binary tree
every vertex has either TWO children or ZERO children
(then it is called a leaf)
with n leaves, in the format
Tree=[] or [T1,T2] (T1 and T2 are smaller trees)
-
NuL(T): the number of leaves of T
-
DressUp(T,L): inputs a "naked" binary tree with
n (say) leaves and a list of numbers and
sticks the numbers in order from left to right
For example
DressUp([[],[[],[]]],[4,1,3]);
should be [[4],[[1],[3]]]
-
C(n): the number of full binary trees with n leaves
-
RollLD(L): inputs a list of NON-NEGATIVE
integers and outputs an integer from 1
nops(L) such that the prob. of getting i
is L[i]/convert(L,`+`):
-
RandBT(n): A uniformly random full buinary
tree with n leaves
Homework for Monday, March 4, 2013 class (due Sunday March 10, 2013
10:00pm)
Create a file hw12.txt and place it in your em13 directory.
-
Write a procedure Tn(n) that inputs a positive integer
n and outputs the set of all full ternary trees with n leaves,
where a vertex can either have no children, or exactly
three children.
-
Write a procedure C3(n) that computes nops(Tn(n)), fast
-
Write a procedure RandTT(n) that inputs a positive integer n
and outputs (uniformly at random) a ternary tree with n leaves.
[Hint: now it is a bit tricky to construct the loaded die, since
you need to consider all triples [a,b,c] that add up to n.
One way is to constuct a list of such triples and
create another list with the corresponding weights].
-
Write a procedure RandBTG(n,a) that inputs a positive integer n
and a positive integer a, and outputs a random gambling tree
where the values at the leaves are integers between 1 and a,
drawn independently and uniformly at random.
-
[Optional, big Challenge!, 2 dollars to the OEIS foundation]
Write a procedure RandGT(n,S) that inputs a positive
integer n and a set of positive integers S, and outputs,
uniformly at random, a tree with n leaves where the number of
children that a vertex may have is either zero, or a member of S.
For example, RandBT(n) that we did in class is RandGT(n,{2})
and RandTT(n) above is RandGT(n,{3}).
Added March 11, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thurs., March 7, 2013
Today we continued gambling trees. The program file is:
and it contains the following procedures
(it also contains the procedures of C12.txt):
-
StPete(n): The gambling tree for the St. Petersburg casino with n leaves
-
PGG(G,x): the prob. gen. function in x playing the gambling tree G
-
PGF(A,x): inputs a prob. dist in the format {[a,p]}
meaning the prob. of getting outcome a is p,
outputs the Prob. Gen. Function add(a[2]*x^a[1], a in A)
-
RandGT(n,L) inputs a pos. integer and an increasing
n-list of numbers outputs a randon gambling tree
whose leaves (final possible outcomes
in our casino) are a random permutation of L
-
ExpV(G): the expected gain playing the
gambling tree G
Homework for Thurs., March 7, 2013 class (due Sunday March 10, 2013,
10:00pm)
Create a file hw13.txt and place it in your em13 directory.
-
Write a program AveExp(L) that inputs a list L, of different numbers and outputs the
average of the expected value of DressUp(T,L) over all trees T with n leaves.
-
Write a program AveExpSeq(f,n,N) that inputs an expression f in n, and
outputs the list of length N, whose m-th entry is
AveExp([f(1), ..., f(m)])
-
[Challenge]
Can you conjecture explicit expressions in n for the following? (I don't know the answer)
-
AveExp([seq(i,i=1..n)]):
-
AveExp([seq(i^2,i=1..n)]):
-
AveExp([seq(2^i,i=1..n)]):
-
AveExp([seq(3^i,i=1..n)]):
-
[Version of March 10, 2013, 10:28 EDT, correcting a typo pointed out by Phil Benjamin,
m[r]:=(x d/dx)^r F(x) below was previously and erroneusly m[r]:=(x d/dx)^r P(x)]
The n-th moment of a random variable X is defined to be the expectation of X^n.
By elelmentary (discrete!) probability theory if P(x) is the prob. generating function
(so we must have P(1)=1), the average is P'(1) (let's call it a)). The
cetralized p.g.f is F(x):=P(x)/x^a. The r-th moment about the mean is
m[r]:=(x d/dx)^r F(x)
(where d/dx is the differentiation operator with respect to x).
In particular m[2] is called the variance and its square-root is called
the standard deviation. Finally the normalized moments are
N[r]:=m[r]/m[2]^(r/2)
Write a program Stat(P,x,R) that inputs an expression P in x, and a pos. integer R larger than 2, and checks that indeed P(1)=1,
and returns a list of length R, whose first entry is the average, its second entry is
the variance, and for r=3..R, the r-th entry is N[r].
-
What is
evalf(Stat(((1+x)/2)^10000,x,10));?
evalf(Stat(((1+x)/2)^100000,x,10));?
evalf(Stat(((1+x)/2)^1000000,x,10));?
Can you conjecture a trend?
-
What is Stat(FGF(StPete(n),x),x,10)? for n from 1 to 20? Can you guess a formula for (i) the average (easy)
(ii) variance (m[2]) (iii) N[r]? (I don't know the answer) in terms of n?
-
I unrecommend
Teddy Niemeic's primer
as being too infinitarian, so read it at your own risk, but please read the first few pages
(at least up to p. 302) of
William Boyce's
beautiful article on the intriguing Shepp Urn, that we will discuss next time
Added March 11, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Monday, March 11, 2013
Today we started to study the Shepp Urn, the file is
and it contains the following procedures
-
LarryE(m,p): input non.neg. integers m and p
representing a box (urn) with m (-1)dollar bills
and p 1-one dollar bills, outputs the
expected gain of the game (where at any given time)
you pick at random one or the other dollars bills
Pr. of picking a dollar bill is p/(m+p) and
Pr. of picking a -1 dollar bill m/(m+p)
And you quit whenever you want,
-
Larry(m,p): input non.neg. integers m and p
representing a box (urn) with m (-1)dollar bills
and p 1-one dollar bills, outputs the
expected gain of the game (where at any given time)
you pick at random one or the other dollars bills
Pr. of picking a dollar bill is p/(m+p) and
Pr. of picking a -1 dollar bill m/(m+p)
And you quit whenever you want
-
LarryTable(N):
-
beta(p): The largest m for which LarryE(m,p) is >0
-
LarrtPGF(m,p,x): The generalized "polynomial"
whose exponets are numbers (rational) and
where the coeff. of x^a is the prob. of exactly getting
a dollars when you end the game (but following the
stupid max. exp. strategy)
-
PrDist(m,p): the list of pairs [i,p(i)] where
the prob. of getting i is p(i)
-
Stat(P,x,K): inputs a prob. g.f. P in x and a pos. integer K
outputs a list of length K where the first entry if the
average, the second the variance, and the r-th is the so-called
alpha coefficiient m[r]/m[2]^(r/2) where m[r] is the
r-th moment about the mean
try:
evalf(Stat(((1+x)/2)^1000,x,10));
[Taken from homework problem, not done in class]
Homework for Monday, March 11, 2013 class (due Sunday March 24, 2013,
10:00pm)
Note:
Officially this homework is due after spring break, but it is
strongly recommended that you do it as soon as possible, so that
you would have a relaxing Spring break, and also while the
material is fresh in your mind.
Create a file hw14.txt and place it in your em13 directory.
-
According to Larry Shepp, (beta(p)-p)/sqrt(2*p) tends to a certain fixed
constant, that we may call the Shepp constant.
Estimate it.
-
Write a program PrWinBreakEvenLose(m,p) that inputs two positive
integers m and p and outputs a list of length 3 whose first
entry is the probability of exiting with a positive amount, zero amount,
and negative amouts, respectively. By examinining
PrWinBreakEvenLose(beta(p),p) for many p, decide whether it
is a good idea to play the Shepp "maximize expectation strategy", or
whether it is too risky, if you hate to lose.
-
Verify empirically the inequalities
proved in
William Boyce's
beautiful article for m,p up to 200.
-
In
William Boyce's
beautiful article, it is proved (Lemma 3.6) that V(1,p)=p^2/(p+1).
Conjecture expressions for V(2,p) and V(3,p), and if possible, prove
them by induction.
-
[Challenge!]
Modify LarryE(m,p) to write a procedure
LoveToWinButHateToLose(m,p,A,pW,B,pL,x)
where you keep playing as long
as your probability of winning more than A dollars is ≥ pW and
your probability of losing more than B dollars is ≤ pL and.
You should use the probability generating funcion with respect to x,
where the answer is 1 if you should stop (i.e. you can't maintain
both conditions).
-
Write a procedure betaApwBpL(p,A,pW,B,pL) that inputs a positive
integer p and outputs the smallest m such that you do not quit,
according to the above strategy with p +1 balls and m -1 balls.
-
Compare the sequence
[seq(betaApwBpL(p,A,pW,B,pL), p=1..40)]
for various A, pW, B, pL to that of
[seq(betaApwBpL(p,A,pW,B,pL), p=1..40)]
-
[Huge Challenge ($100 from Larry Shepp and another $100 from me to the
OEIS foundation). ]
Prove that Larry(m,p)=0 if and only if m=2 and p=1.
Added April 2, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thurs., March 14, 2013
Today we celebrated Pi Day! The file is
It contains the following procedures:
-
GL(N): Possibly one of the slowest ways to appx. Pi
(The Gregory-Leibnitz series, 4*arctan(1))
-
JG(N): The N-th partial sum of the amazing Jesus Guillera series for
128*sqrt(5)/Pi^2
(Possibly one of the fastest ways to appx. Pi)
-
WZpi(n): an OK (but not great) way to appx. Pi
-
CoinToss(N) the outcome of tossing a fair coin
randomly N times winning a dollar with Heads
and losing with Tails
-
MetaCoinToss(N,K): The number of times of breaking even
if you "toss a coin N times" every day for K days
-
IsUD(pi): Returns true if the permutation pi Up down
and false otherwise
-
AndreS(n): The number of Up-Down perms of length n
(the very slow way!)
-
AndreF(n): The number of Up-Down perms of length n
(a pretty fasat way!)
Homework for Thurs. March 14, 2013 class (due Sunday March 24, 2013,
10:00pm)
Note:
Officially this homework is due after spring break, but it is
strongly recommended that you do it as soon as possible, so that
you would have a relaxing Spring break, and also while the
material is fresh in your mind.
Create a file hw314159.txt and place it in your em13 directory.
-
Look through
Jesús Guillera's amazing PhD thesis and
pick trunc(2*Pi) of your favorite summations, and check them
empirically, (by writing procedures similiar to JG(N)),
label them JG1(N), ..., JG6(N),
and find out how many more digits they give for every turn.
-
Optimize JG(N) by writing a procedure JGfast(N), that
takes the summand (w/o the polynomial) part, let's call it
f(n), and once and for all (by hand or by Maple, finds
f(n)/f(n-1) in simplified form, let's call it g(n), and
then computes the successive terms by doing f(n)=f(n-1)*g(n).
By using time(), see how much faster is JGfast(1000) then
JG(N)
-
Implement the
Bailey-Borwein-Plouffe famous formula in the usual(decimal) notation,
by writing a procedure BBP(N) that uses the first N terms.
-
[Challenge], Write a procedure BPP2(n) that inputs a pos. integer n
and computes the n-th digit of Pi in BASE 16, WITHOUT COMPUTING
THE PREVIOUS DIGITS
-
Using
Sum(a^i/(8*i+j+1),i=0..infinity)=int(x^j/(1-x^8)),x=0..a)
(why?), rewrite the B-B-P summation as an integral of a certain
rational function from 0 to 16, and use Maple to prove it
(RIGOROUSLY!)
-
-
[Human part]
If you look at a typical Up-Down permutation of length n,
the largest entry, "n", must be
somewhere, and its location should be even (why?), say it
is 2*k (1<=k<=trunc(n/2)), the part to its left and the
part to its right are both Up-Down permutations of
length 2*k-1 and n-2*k, respectively, but since their set of
elements is NOT the minimal one, there are binomial(n-1,2*k-1)
ways of picking what entries amongst {1, ..., n-1} would go to
the left
(and the rest go to the right). Convince yourself that
we have the non-linear recurrence
A(n)=Sum(binomial(n-1,2*k-1)*A(2*k-1)*A(n-2*k),k=1..trunc(n/2))
[Added March 15, 2013: this corrects some typos from yesterday, sorry.]
-
Use this recurrence to write a procedure AndreFF(n)
(don't forget option remember)
How does it compare, speed-wise, with AndreF(n) we did in class,
for n=1000?, n=4000?
-
Write a procedure AndrePi(n) that computes 2*n*A[n-1]/A[n],
and see how close it gets to Pi for, n=10, 100, 1000.
-
Get ready to pick a project soon (I will post suggestions,
next week).
-
Pick a
final project
(must declare a project by March 31, 2013, 10:00pm.)
-
Happy Spring break.
Added April 2, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Monday, March 25, 2013
Today we did Wall-Street Gambling
It contains the following procedures:
-
N(x): the chance of being <=x if you belong to
a Normal population with mean 0 and variance 1
-
d1(S,t,r,sig,K,T): the quantity d1 featuring in the BS formula,
where
S=price of stock, t=time, r=interest rate for no-risk
goverment bond,
sig=volatility, K=strike price (how much
you can buy it for at time T),
T=end of game
-
d2(S,t,r,sig,K,T): d1(S,t,r,sig,K,T)-sig*sqrt(T-t):
-
C(S,t,r,sig,K,T): The value of a
European Call option according to the infamous
B-S formula.
S=price of stock, t=time, r=interest rate for no-risk
goverment bond, sig=volatility, K=strike price (how much
you can by it for at time T), T=end of game
-
PBS(): A completely rigorous proof of a formula that
won the Nobel (Memorial) prize (in Economics)
-
DC(S,t,r,sig,K,T,n): The appx. value (with n discrete
steps) using the Cox-Ross-Rubinstein Discrete Version
of B-S for a
European Call option according to the infamous
B-S formula.
S=price of stock, t=time, r=interest rate for no-risk
goverment bond, sig=volatility, K=strike price (how much
you can by it for at time T), T=end of game
Homework for March 25, 2013 class (due March 31, 2013,
10:00pm)
Create a file hw16.txt and place it in your em13 directory.
-
[5 dollars award for the first person to do it, DONE!]
Fix DC(S,t,r,sig,K,T,n) in such a way that
DC(S,t,r,sig,K,T,1000) would be VERY close to
C(S,t,r,sig,K,T)
[Added March 25, 2013: Joey Reichert won this prize, so it is cancelled]
-
Experiment with C(S,t,r,sig,K,T) and DC(S,t,r,sig,K,T,n), and figure out how big n has
to be approximate C(S,t,r,sig,K,T) well. Also, experiment with various sigmas
-
[Challenge, need to know some probability] Prove ("rigorously") using the Normal approximation the
binomial distribution that DC(S,t,r,sig,K,T,n) converges to C(S,t,r,sig,K,T) as n goes to infinity.
-
Read
Manuel Blum's classic paper on
Zero-Knowledge Proofs.
Added April 3, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thursday, March 28, 2013
Today we did the prover-verifier game as beautifully described in
Manuel Blum's classic
paper on Zero Knowledge Proof.
C17.txt
NOTE: THIS CODE IS FLAWED (thanks to Phil Benjamin. SEE HOMEWORK BELOW!)
Contains procedures:
-
RandHG(n,p): inputs a pos. integer n and a rational number
p and outputs a graph on n vertices that has a Hamiltionian
cycle, and the remaining edges are chosen independently
with prob. p
-
Options(G,pi,sig): Given a Hamiltonian graph G,
and a proof of its Hamiltonicity in the form of a
a permutation pi (n=nops(pi)) and,
a permutation sig, prepare both options:
-
either show your enemy the permutation sig and the image of the graph
under sig
-
Show the Hamiltonian cycle and the transformed graph (BUT not sig)
-
OneStep(G,pi): inputs a graph that you know is Hamiltonian,
together with a permutation pi on the vertices, one step
in the probabilistic protocol described in Manuel Blum's article
to convince a verifier that G is Hamiltonian, and that you know
a proof (i.e. you know the Hamiltonian cycle) WITHOUT giving
pi away.
-
VerifyOneStep(G,i,MrNakamura): Inputs a graph G, i=1 or i=2,
and a pair MrNakamura (whose content) depends on i, and outputs
true or false.
Homework for March 28, 2013 class (due March 31, 2013,10:00pm)
Create a file hw17.txt and place it in your em13 directory.
-
According to Phil Benjamin (see his
insightful message),
the code we did in class was not quite right. I agree. Part of the problem is the shortcut
of representing graphs as a set of edges rather than in terms of the n by n 0-1 matrix,
the adjacency matrix. Fix the code to conform exactly to what Manuel Blum told us to do.
-
After fixing it, finish it up, by writing a procedure VerifyManySteps(G,i,MrNakamura,NumberOfSteps),
where NumberOfSteps is the number of rounds.
-
Test it with several random graphs.
-
Look up procedure Finance[BlackScholesPrice] in Maple (with d=0) and test it against
C(S,t,r,sig,K,T) of
C16joey.txt.
-
Modify DC(S,t,r,sig,K,T,n) of
C16joey.txt
to write a procedure DCX(S,t,r,sig,K,T,n,X), that does not only give you the
expected gain (hence the "fair price" of the option to buy the stock at price K at time T, if
the price of the stock at time t is S, with the given r and sig), but gives you the generalized polynomial
(where the powers are any rational numbers), such that the coefficient of X^r is the
probability that you would make exactly r dollars from the deal.
-
Pick a
final project,
write a short summary/abstract, and put it in a file called proposal.txt in your secret em13 directory.
Added April 3, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Monday, April 1, 2013 Class
Initially,
Myron Scholes,
was scheduled to give a guest lecture in this class, but he couldn't make it due to a previous
commitment to play golf. Luckily, Edinah Gnang
graciously agreed to fill-in, and gave a fascinating introduction to SAGE, following his interesting
writeup.
Homework for April 1, 2013 class (due April 7, 2013,10:00pm)
Create a file hw18.txt and place it in your em13 directory.
-
Make up some neat examples of Sage code
-
[Optional, added April 4, 2013]
Do
Edi Gnang's assignment
Added April 8, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thursday, April 4, 2013
Today we started "real" Game Theory, in particular Nash Equilibrium.
It contains programs (some still not even started)
Homework for April 4, 2013 class (due April 7, 2013,10:00pm)
Create a file hw19.txt and place it in your em13 directory.
Homework for April 4, 2013 class (due April 7, 2013,10:00pm)
Create a file hw19.txt and place it in your em13 directory.
-
[Human reasoning]: Find analytic expressions for
BRA([[R,S],[T,P]],y) in terms of R,S,T,P,
finding the condition on y to make it
1,0, or "all x between 0 and 1".
Do the same for BRB([[R,T],[S,P]],x)
-
[Human reasoning]
Find conditions on R,T,S,P for which a Nash Equilibrium [x0,y0]
is
-
[1,1] (the happy case!)
-
[1,0]
-
[0,1]
-
[0,0]
-
Neither (i.e. it is not pure strategy). Find explicit expression
for when that happen
-
Test the above at random values of R,T,S,P
using our numerical BRA(A,y) and BRB(B,x)
Added April 8, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Monday, April 8, 2013
Today we continued "real" Game Theory, studying Nash Equilibria
for two players but for an arbitrary number of (finite!) options
for both players.
-
MaxF(F,x,n): inputs sumbol x, positive integer n
and a linear combination of x[1], ..., x[n]
outputs the set of i such that x[i]=1 makes F max
For example,
MaxF(3*x[1]+4*x[2]+4*x[3]+2*X[4],x,4);
should output {2,3}.
-
BRrow(A,y): The best response by the ROW player
(whose payoff matrix is the m by n matrix A)
to the strategy y of the column player.
The output is a set of pure strategies.
-
BRcol(B,x): The best response by the COL player
(whose payoff matrix is the m by n matrix B)
to the strategy x of the row player.
The output is a set of pure strategies
-
PNE(A,B): The set of pure Nash Equilibria of
a game where the ROW player's payoff matrix is
A and the COL player's payoff matrix is B.
A and B are given as lists of lists of the same size.
-
GenRanMat(m,n,K): a random m by n matrix with
entries from 0 to K
-
tra(A): the transpose of A (given as a list of lists)
Urgent Homework for April 8, 2013 class (due April 11, 2013,11:50am)
Read William Press and Freeman Dyosn's recent
seminal article about
the iterated Prisoner's dilemma.
[But according to Nathaniel Shar you should take the evolutionary claims with a
grain of salt.]
Usual Homework for April 8, 2013 class (due April 14, 2013,10:00pm)
Create a file hw20.txt and place it in your em13 directory.
-
Write a program GatherStat(n,K,L,t) that inputs positive integers
n, K, and L, and generates L times a random n by n matrices A,
whose entries are from 0 to K, and outputs a list whose
(i+1)-th item is the number of these that yield
i Nash Equilibria for (A,tra(A)). In particular
the first entry is the number of these matrices that
have no pure equilibria.
-
For modest n and K, but large L (say 1000), see if you can
conjecture the probability of having (i) no pure Nash equilibria
(ii) exactly one, etc.
-
Modify PNE(A,B) to write a procedure NE2(A,B,K) that
inputs 2 by 2 matrices A, and B, and
looks for not necessarily pure pairs [x,y] of Nash equilibria
where y is the form [i/K,1-i/K] for i between 0 and K.
-
Modify PNE2(A,B,K) to write a procedure NE3(A,B,K) that
inputs 3 by 3 matrices A, and B, and
looks for not necessarily pure Nash equilibria pairs [x,y] of Nash equilibria
where y is the form [i/K,j/K,k/K] for i,j,k between 0 and K
and, of course i+j+k=K
-
[Optional Challenge]
Write a procedure NE(A,B,K) that
inputs matrices A, and B of the same size, and
looks for not necessarily pure pairs [x,y] of Nash equilibria
where y is the form [i[1]/K,i[2]/K, ..., i[n]/K] for
i[1],i[2], .. i[n] between 0 and K
and, of course i[1]+ ...+i[n]=K
-
Write a procedure LetsCooporate(A,B,i,j) that inputs matrices A and
B of the same size, i, j such that [i,j] is a pure Nash equilibrium,
and outputs the (possibly empty) set of pairs [r,s] such that
A[r][s]>A[i][j] B[r][s]>B[i][j]
In other words, if the "prisoners" had a secret cell phone or
otherwise could communicate, they could propose to each other
to choose strategies that are not Nash Equilibria, but
are better for both of them (like to both cooporate in the famous
original Prisoners' Dillema).
Added April 15, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thurs., April 11, 2013
Today we started to implement the Iterated Prisoner's
Dillema following the Press-Dyson
mini-revolutionary article.
Contains procedures:
-
SS(T): inputs a stochastic matrix T
given as a list of lists, where T[j][i]
is the prob. of moving in one step from
state j to state i.Outputs a prob. vector v
of length nops(T), such that after googol days
your prob. of being at state is v[i]
i.e. the solution of v T= v
-
Apq(p,q): Given the prob. strategy of X, p
and of Y, q, finds the transition matrix
of the iterated Prisoner's Dillema
(or any two-person, two-option game).
-
sX(R,S,T,P,p,q): The ultimate expected pay-off per day of play
of player X with the paramers R=Rewared, S=Sucker, T=Temptation,
P=Punishment, and with the (probabilistic) strategies p and q for
player X and Y respectively.
-
sY(R,S,T,P,p,q): The ultimate expected pay-off per day of play
of player Y with the paramers R=Rewared, S=Sucker, T=Temptation,
P=Punishment, and with the (probabilistic) strategies p and q for
player X and Y respectively.
-
FixIncomeY(R,S,T,P,p,Amount):
Finds those p (in terms of parameters p[1], ..., p[4]
that forces player Y to always get Amount no matter
what q is
Homework for April 11, 2013 class (due April 14, 2013,10:00pm)
Create a file hw21.txt and place it in your em13 directory.
-
Modify procedure
FixIncomeY(R,S,T,P,p,Amount):
to write a procedur
FixRatioXY(R,S,T,P,p,Ratio):
that inputs numbers (or symbols!) R,S,T,P, a symbol p, and
a number Ratio, and outputs a general vector p=[p[1],p[2],p[3],p[4]]
such that
sX(R,S,T,P,p,q)=Ratio*sY(R,S,T,P,p,q):
For ALL q.
-
[A little challenging]
Using Loaded Dice, write a program to simulate a random walk
on a graph given by a transition matrix M, and see if the
fraction of the time in the long-run that random walker spends at each
site (vertex) is close to those predicted by SS(M).
It can start anywhere it wants.
Added April 15, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Monday, April 15, 2013
Today we started to implement gambling theory,
doing the classical Gambler's Ruin.
Contains procedures:
-
LC(p): simulating a loaded coin
that returns true with prob. p, p rational
-
GRgame(p,s,N):inputs prob. p of success (a rational number)
starting capital s, and exit capital N (if you are a winner)
once you have 0 dollars you must leave the casino.
Returns, true (false) if you won (or lost) followed
by the number of rounds it took
-
ManyGames(p,s,N,HowMany): repeats GRgame(p,s,N)
HowMany times, and returns the percentage of times
the gambler won, and the average duration if it won,
followed by the average duration if it lost followed
by the average duration
-
VW(p,N): inputs prob. of winning one round, p
and the exit capital N, outputs the list of
size N-1 whose i-th entry is the EXACT prob.
of winning if you currently posses i dollars
-
VD(p,N): inputs prob. of winning one round, p
and the exit capital N, outputs the List of
size N-1 whose i-th entry is the EXACT EXpected Duration
of winning if you currently posses i dollars
Homework for April 15, 2013 class (due April 21, 2013,10:00pm)
Create a file hw22.txt and place it in your em13 directory.
-
[Human Reasoning, No peeking!] Using discrete math 101
(the solution of a recurrence of the form
Ax(n-1)+Bx(n)+Cx(n+1)=0
is a linear combination
c1 αn+c2 βn
where α and β are roots of the auxiliary equation
A/x+B+Cx=0 )
find a formula, in terms of the symbols p,i,N,
for the probability of winning a Gambler's ruin game, with a loaded
coin with prob. of winning is p if your current capital is i and
the exit capital is N.
[In other words, find the exact values of the component of
VW(p,N)]
[Hint, to find c1 and c2 use the boundary
conditions x(0)=0, x(N)=1]
-
Check your formula against VW(p,N) for numeric N, (N between 2 and 100),
but symbolic p.
[Hint: since you would get radicals, you may need to use
the Maple command `simplify'. If this does not work, just
plug-in numeric values for p, and use evalf]
-
[Human Reasoning, No peeking!] Using discrete math 102
(the solution of a recurrence of the form
Ax(n-1)+Bx(n)+Cx(n+1)= D(n)
is a linear combination
c1 αn+c2 βn +PS(n)
where PS is a particular solution (similar-looking
to the right side, D(n))
of the above (inhomogeneous) linear recurrence ,
and where α and β are roots of the auxiliary equation
A/x+B+Cx=0 )
find a formula, in terms of the symbols p,i,N,
for the expected duration of a Gambler's ruin game, with a loaded
coin with prob. of winning is p if your current capital is i and
the exit capital is N.
[In other words, find the exact values of the component of
VD(p,N)]
[Hint, to find c1 and c2 use the boundary
conditions x(0)=0, x(N)=0]
-
Check your formula against VD(p,N) for numeric N, (N between 2 and 100),
but symbolic p.
[Hint: since you would get radicals, you may need to use
the Maple command `simplify'. If this does not work, just
plug-in numeric values for p, and use evalf]
-
Modify VD(p,N) to write a procedure VWD(p,N) that inputs a positive
integer N and a prob. of winning (a single coin toss)
p (numeric or symbolic) and outputs the expected duration of a game
conditioned on winning the game.
-
Use VDW(1/2,N) for various N to conjecture an explicit formula
(a polynomial in i and N) for the expected duration of a winning
game with a fair coin that starts with i dollars and ends with N dollars.
-
[Human and/or Computer reasoning]
Prove this formula rigorously.
Added April 22, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thurs., April 18, 2013
Today we continued to implement gambling theory,
extending the classical Gambler's Ruin to
following an arbitrary strategy and figuring out the
best one. (There are times to be timid and there are times to be bold).
Contains (in addition to the procedures of C22.txt) procedures:
-
GGgame(p,s,N,L):
Inputs
-
prob. p of success (a rational number) of a single coin toss
-
starting capital s
-
exit capital N
-
A list of integers L, such that 1<=L[i]<=min(i,N-i)
Outputs true (false) if you won (collecting N dollars)
or lost (leaving with 0 dollars), followed
by the number of rounds it took to finsih.
-
ManyGGgames(p,s,N,HowMany,L): repeats GGgame(p,s,N,L)
HowMany times, and returns the percentage of times
the gambler won, and the average duration if it won,
followed by the average duration if it lost followed
by the average duration
-
VGW(p,N,L): inputs prob. of winning one round, p
and the exit capital N, outputs the List of
size N-1 whose i-th entry is the EXACT prob.
of winning if you currently posses i dollars
following the strategy L
-
VGD(p,N,L): inputs prob. of winning one round, p
and the exit capital N, outputs the List of
size N-1 whose i-th entry is the EXACT EXpected Duration
of winning if you currently posses i dollars
and follow strategy L
-
IsAB(L1,L2): is the list of integers L1 larger or equal,
componentwise than the list of integers L2?
-
AllVecsLess(L): The SET of all the vectors M of length nops(L)
of pos. integers such that 1<=M[i]<=L[i]
-
AllStra(N): The set of all possible strategies with maximal capital N
-
BestStra(p,N): The set of best strategies with max. capital N
and prob. of winning a round p.
Homework for April 18, 2013 class (due April 21, 2013,10:00pm)
Create a file hw23.txt and place it in your em13 directory.
-
Given an expression f of a variable p, Maple can solve the
inequality
f ≥ 0
(and of course)
f ≥ g
(for g another expression in p) for p between 0 and 1 by using the command
solve({f>=0, 0<=p, p<=1},{p})
Use the command to write a program TimidVsBold(N,p) that inputs
a positive integer N and a SYMBOL p, and outputs a list of length N-1
whose i-th item is the set of p's for which the Timid Strategy
(always betting 1 dollar) is better than the Bold Stragegy
(betting min(i,N-i) dollars).
Find TimidVsBold(N,p) for N between 3 and 20.
-
Define the set of neighbors of a strategy, L, to be the set
of (legal!) strategies where all the components are the same
except at one entry, where it is off by one (so there must be
at most 2*nops(L) neighbors). Write a procedure
BestNeig(N,p,L,i) that inputs N, (numeric) prob. p, a stragegy L
(a legal list of integers of length N-1), and an integer i
between 1 and N-1, and outputs the best neighboring strategy
in the sense that the prob. of winning with that p and N, if
you currently have i dollars is highest. If L is better than all its
neighbors, return L.
-
By iterating BestNeig(N,p,L,i), write a procedure
BestStraFast(p,N,i)
that inputs p and N and i as above and outputs the
strategy that maximizes the probability of winning if you have i dollars.
(well, at least finds a local maximum that is hopefully a global one).
-
Experiment with BestStraFast(p,N,i) and convince yourself that
if p is less than 1/2 it is good to be bold, and if p is larger than 1/2,
it is good to be timid.
-
Read and understand
Karl Siegrist's engaging and lucid article
How to gamble if you must
[
errata (compiled by Phil Benjamin)].
Added April 22, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Monday, April 22, 2013
Today we concluded Gambling theory by finding out
how to gamble if you are in a hurry.
Contains (in addition to the procedures of C22.txt and C23. txt)
procedures:
-
TiS(N): the timid strategy
-
BoS(N): the bold strategy
-
KeS(N,f): The Kelly Strategy with fraction f of your capital
-
BestSDD(N,i,p,T): How much should you bet if you
currently have i dollars, but need to come up
with N dollars in ≤ T rounds (or else you would die)
and your chance of winning one round is p
Outputs the stake that would maximize your chances
of staying alive, followed by your chance of doing so
-
OptimalGamblingStory(N,i,p,T): Inputs as above,
simulate a gambling scenario following
the advice of BestSDD(N,i,p,T) to the happy or
bitter end, and tells you at each step your probability
of surviving. [Lightly edited and embellished after class]
Homework for April 22, 2013 class (due April 28, 2013,10:00pm)
Create a file hw24.txt and place it in your em13 directory.
-
Write a procedure
OptimalGamblingTerse(N,i,p,T) ,
that does not tell you the story but returns TWO outputs:
-
true (or false) if you came out alive (or dead)
-
A list of the intermediate probabilities (of length less than T
in the lucky case, and of length T in the sad case)
-
Write a procedure
ManyOptimalGambling(N,i,p,T,K) ,
that runs
OptimalGamblingTerse(N,i,p,T)
K times and returns TWO NUMBERS (in floating point)
-
the ratio of times the gambler came out alive
-
the "closest call", the lowest ever (intermediate) probability
of survival that ever showed up in the games that
ended happily.
-
Try out
ManyOptimalGambling(20,9,11/20,10,5000);
How close is the first output to
BestSDD(20,9,11/20)[2]; ?
What is the second output? (of course different students should
get different answers!)
-
Write a procedure
BestKellyFactor(N,i,p,TimeIsMoney,resolution)
where
-
N is the ultimate desired capital
-
i is your current capital
-
p is the probability of a succesfull one round
(usually larger than 1/2)
-
TimeIsMoney is the cost that you have to pay by playing a single
round
-
resolution is a small positive number (say 1/200)
that tries out all the possible Kelly strategies, starting with
f=resolution, and then f=2*resolution, ..., and so on,
all the way to the bold
scenario f=1, and finds the f that maximizes
(VGD(p,N,TiS(N))[i]-VGD(p,N,KeS(N,f)[i]))*TimeIsMoney
-(VGW(p,N,TiS(N))[i]-VGW(p,N,KeS(N,f)[i]))*N
In other words you maximize the expected saving from
exiting the casino sooner than you would with timid play,
take away the expected extra
loss due to your impatience and unwilling to play the
boring timid strategy.
[Note: If your probability of ultimately winning a gambling game where
you start out with i dollars is P (not to be confused with
the prob. of winning a single round p), then your
expected gain is (N-i)*P-i*(1-P)=N*P-i. In the above expression i cancells
out when we compare to strategies]
-
Write a procedure
WhatKellyCharges(p,N,i)
that finds out how much is TimeIsMoney when f is Kelly's recommended
value f=2*p-1
What are
-
WhatKellyCharges(11/20,20,10)?
-
WhatKellyCharges(11/20,40,21)?
-
WhatKellyCharges(11/20,400,150)?
Added April 29, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thurs., April 25, 2013
Today we were fortunate to have a guest lecture by
the greatest Maple guru in the world,
Frank Garvan
who showed us how to discover brand-new amazing analogs
to the number of partitions of n into distinct parts
of Ramanujan's famous congruences for the number of partitions of n.
Any specific one can be rigorously proved by computers, but
the (very probable) infinite family (the exact conjencture still
has to be made precise) still needs humans.
Contains many procedures. See the comments in the source code.
Frank Garvan produced the following
maple worksheet during the class.
Homework for April 25, 2013 class (due April 28, 2013,10:00pm)
Create a file hw25.txt and place it in your em13 directory.
-
Use procedure prodmake in
C25.txt
to conjecture a product formula for the exponential function.
Use human ingenuity to prove it.
-
Procedure vecptns(n) gives the set of triples ("vector partitions")
[Q,P1,P2]
where Q is a partition into distinct parts, and P1, P2 are partitions
such that all the parts of Q,P1,P2, add up to n.
Is the sequence "number of vector partitions on n" in the OEIS?
Is it for that reason? If it is for another reason, prove that
they are actually the same for all n, not just the first 30 terms.
-
[Added Aprl. 28, 2013]: Read and understand
Wikipedia's entry on the
centepede game
and Dr. Z.'s Opinion 70
Added April 29, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Monday, April 29, 2013
Today we developed generalizations of the
centepede game.
Contains procedures:
-
ElsterG(): Jon Elster's graph of the centepede game described in
Opinion 70, taken from his article mentioned there, using
our data structure
-
P1(G,i): The best move, followed by the pair
[payOff of I, payoff of II] in a game given
by the generalized directed graph G where
the sinks are labeled by payoof pairs,
If player I is moving
-
P2(G,i): The best move, followed by the pair
[payOff of I, payoff of II] in a game given
by the generalized directed graph G where
the sinks are labeled by payoof pairs,
If player II is moving
-
RandDi(n,k): a random directed graph with n
vertices and (usually) outdegree k
-
RandGame(n,k): inputs pos. integers n and k
and outputs a random digraph that has the
property that out of vertex i all the labels
are less i, and there are (usually) k outgoing
edges (taken from C2.txt)
-
RGG(n,k,Si,L): a random centepede-style
game with n vertices, roughly outdegree k
for non-sinks, and roughly Si sinks
and payoffs taken randomly to be integers from 0 to L
Homework for April 29, 2013 class (due May 5, 2013,10:00pm)
Create a file hw26.txt and place it in your em13 directory.
-
Write a program
AllBetterDeals(G,i)
That inputs a generalized game G (given in our format) and
a vertex i ( an integer between 1 and nops(G)) that
is not a sink, and finds ALL the smaller
subgames H (obtained by deleting some of the edges, i.e. replacing
some (or all) of the entries of G that are sets by proper subsets
and such that the payoffs for BOTH players are higher if
Player 1 starts at vertex i. You also want the improved payoffs.
So the output should be the a set of pairs {[H,[a,b]}
What is
AllBetterDeals(ElsterG(),1)?
[Note: Of course this is exponential time in the number of edges
so don't try it on too large graphs.]
-
Experiment with ten random graphs with n=6, k=2,
Si=3, and L=10 (using RGG(n,k,Si,L) ).
-
Write a program
OneBetterDeals(G,i,K)
That inputs as above, and tries to find, by
repeated random deletion of edges,
K times,
ONE subgame that is a better deal for BOTH players if
player 1 starts at vertex i. If no better deal is found, return FAIL.
-
Experiment with ten random graphs with n=20, k=3,
Si=4 and L=10 (using RGG(n,k,Si,L)) .
-
[Due May 2, 2013] Read and understand
Wikipedia's entry on the
Braess paradox.
and
Frank Kelly's beautiful and insightful
article on traffic networks
Added May 6, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Programs done on Thursday, May 2, 2013
Today we first wrote a linear-algebra solver that obviates
the need for Linear Algebra and then we started to figure
out how to compute Wardrop equilibria following
Frank Kelly's beautiful and insightful
article on traffic networks, but avoiding matrices.
C27.txt contains
-
MySolve(eq,var): inputs a seq of equations eq in a set
of variables var, and outputs a general solutioin or FAIL
For example MySolve({x+y=5, x-y=1},{x,y}); should return
{x=3,y=2}
C27a.txt (under constuction), contains
-
IsEdgeRoute(Edge,Route): Is edge Edge part of route Route?
-
PotF(J,R,DelayJ,y,x): the Traffic potential function with
edges J1 and Routes R1
(incomplete!)
Homework for May 2, 2013 class (due May 5, 2013,10:00pm)
Create a file hw27.txt and place it in your em13 directory.
-
[somewhat challenging]
Finish implementing the finding of Waldorp equilibria,
by writing a procedure
Waldorp(J,R,DelayJ,y,S,Sspec):
that inputs
-
J, a list of edges (in the form of pairs [i,j]),
-
R, a list of Routes (in the form [i1,i2,i3, ..., ir] where
i1 is the start, ir the end, and i2,i3, .. i(r-1) the
intermediate cities),
-
DelayJ (of the same length as the
list J), of delay expressions phrased in terms of the variable y.
-
a variable y (see above)
-
a list S of allowed pairs [start,end] (possibly with only one member)
-
A list Sspec, of the same size as S, indicating how many cars
have to move from start to end for any pair in S (correponding)
and outputs a list of length nops(R) that indicate the
flows along each route (corresponding to the list R), and
another list with the corresponding delay times for each route.
For example, for the wikipedia graph before the shortcut the
function call should be (calling START: 1, A:2, B:3, END: 4)
Waldorp([[1,2],[2,4],[1,3],[3,4]],
[[1,2,4],[1,3,4]],[T/100,45,45,T/100],T,[[1,4]],[4000]);
and with the "shortcut"
Waldorp([[1,2],[2,4],[1,3],[3,4],[2,3]],
[[1,2,4],[1,3,4],[1,2,3,4]],[T/100,45,45,T/100],T,[[1,4]],[4000]);
and the outputs should be (respectively)
[2000,2000],[65,65]
and
[0,0,4000],[0,0,80]
Make sure that you got it (with your procedure).
-
Find the right input and verify the output for the
example for Braess's paradox given in Frank Kelly's article
(before and after the "shortcut")
-
Construct more complicated traffic networks
(and delay functions) and see whether
you can find the Braess phenomenon.
-
[Challenge, 10 dollars to the OEIS] (modified May 2,8:25pm, 2013)
Edit the wikipedia entry
Braess paradox, to make it more readable.
Added May 6, 2013: See the nice
solutions by members of the class that kindly allowed them
to be made public.
Final Showdown: May 6, 2013
Today we had a final programming contest.
Compute the Grundy values, G[m,n] of the initial position of m by n
Chomp for
1 ≤ m,n ≤ N for as large N as possible. See whether the sequence
{G[i,i]} is in the OEIS, and for which i0, are the sequence {G[i0,i]} already
in the OEIS. During the weekend Shalosh found the following beginning
for the sequence {G[i,i]} (for i from 1 to 12)
[0, 2, 5, 6, 6, 13, 20, 13, 20, 19, 23, 41]
The three winners were:
Patrick, Matthew, and Tim each won a chocolate bar (none of whose pieces is poisonous, I hope).
They promised to submit these sequences to the OEIS, and as soon as they are accepted, I will
link to them here.
[Added May 10, 2013: Thank you, Pat Devlin, for contributing A225366
to the OEIS (and to Matthew Russell for approving it)]
Most of the other students finished it before the class was up. I am very proud of you!
Tomas Vyskocil came up with beautiful sage code.
Optional Challenge Homework for May 6, 2013 class
-
(10 dollars to the OEIS)
The sequence of Grundi values for 2 by n Chomp seems, empirically to be the same as
A001651 (all integers NOT divisible by 3) (I checked it up to n=100).
Prove that it is true for all n.
-
(30 dollars to the OEIS)
Find the Grundy value of the inital position of Chomp for a 3 by 100000 chocolate bar.
Added May 7, 2013: Explore the students' exciting
Final Projects.
Added May 21, 2013:
Here are the students'
evaluations.
Dr. Z.'s teaching page