Math 640: EXPERIMENTAL MATHEMATICS
Spring 2010 (Rutgers University) Webpage
http://sites.math.rutgers.edu/~zeilberg/math640_10.html
Last Update: May 17, 2010.
- 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": classical articles, including
-
Dr. Zeilberger's Office: Hill Center 704 (732 445 2390 Ext. 1326)
-
Dr. Zeilberger's E-mail:
zeilberg at math dot rutgers dot edu (you MUST have MathIsFun in the message,
orExpMathRocks)
-
Dr. Zeilberger's Office Hours: 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 PRIMES. We will
learn how to tell whether any given integer is a prime in
polynomial time (thanks to AKS), really understand the
elementary proof of the Prime Number Theorem,
and, who knows?, may be one of you will prove the Goldbach conjecture?
But the actual content is not that important,
it is mastering the methodology of computer-generated and
computer-assisted research that is so crucial for your future.
There are no prerequisites, and no previous programming knowledge is assumed.
Also, very little overlap with previous years.
The final projects for this class may lead to journal publications.
Added March 22, 2010:
Pick a
project !
(First draft due May 5, 2010, 11:59:59pm).
Added April 7, 2010:
Read Kellen Myers translation of
Rene Taton's
tribute to Viggo Brun (1885-1978). One of the heroes of this course.
Diary and Homework
Programs done on Thurs., Jan. 21, 2010
-
jan21.txt ,
Contains procedures
-
Add1(x,y): inputs two integers in unary in the form
[1,1,1,...,1], and outputs their sum
-
Mul1(x,y): inputs two integers in unary in the form
[1,1,1,...,1], and outputs their product
-
The assignment of novices to experts.
Note that
Kellen Myers has an even
better implementation
, and he calls them unaryadd(x,y) and unarymult(x,y) .
Homework for Thurs., Jan. 21, class (due Jan. 25, 2010)
-
(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.
-
(For experts, optional for novices)
Write programs Sub1 for subtraction and Div1 for division
(outputting the quotient and remainder).
-
(For experts only)
Write (from scratch!) addition and multiplication programs
for positive integers represented in binary, as a list of 0's
and 1's (for example 5 is [1,0,1]).
-
(For everyone)
Read carefully, and understand, pages 225-230 (sections 1 and 2) of
Norman Levinson's paper,
Programs done on Monday, Jan. 25, 2010
-
jan25.txt ,
Contains procedures
-
SA(n): The list of primes ≤ n
-
IsDiv(n,L): inputs a pos. integer n and a list
of pos. integers L, and outputs true if
n is divisible by any of the members of L
-
OneStep(L): inputs the list of the first primes
and outputs the same list with one more prime that
comes right after
-
PrimeList(N): The list of the first N prime numbers, our way
-
PrimeListM(N): The list of the first N prime numbers, the Maple way
Homework for Monday, Jan. 25, class (due Jan. 28, 2010)
-
(For Novices only)
Read and do all the examples, plus make up similar ones,
in pages 30-60 of Frank Garvan's awesome Maple booklet.
-
(Optional for Novices, required for Experts)
Write a program , IF(n), that inputs a positive integer, and outputs
the prime-factorization p1k1p2k2p3k3 ...
as a list of pairs [[p1,k1],[p2,k2],[p3,k3], ... ].
For example, IF(600); should output [[2,3],[3,1],[5,2]] .
Compare your output with ifactor(N);
-
(for everyone)
Define π(n) to be nops(SA(x)); write a program that inputs
a pos. integer N and outputs the list of π(n)/(n/log(n))
for i=1..N. Do you see a trend?
-
Challenge( 5 dollars for the best entry): write a faster
program for PrimeList(N) ;
(Note: Brian Nakamura told me that the Sieve of Eratosthenes
(SA(N) in our program)
is pretty good for the first million primes, so try to use it!)
-
(For everyone)
Read (again!) carefully, and understand, pages 225-230 (sections 1 and 2) of
Norman Levinson's paper.
Programs done on Thursday, Jan. 28, 2010
-
Jan28.txt ,
Contains (in addition to those of Jan25.txt, listed in
HelpOld();) the new procedures (listed in Help();)
-
PNT(n): The ratio of #primes <=n and n/log(n)
(it should go to 1, if the the Prime Number Theorem is
correct)
-
EucSeq(n) The first n terms in the sequence of primes
implementing Euclid's "proof" that there are infinitely
many primes
-
TPL(n) The smaller sister of the first n pairs of twin-primes.
-
piJ(x): the number of primes less than x .
-
thetaJ(x): the sum of the ln of primes less than x .
-
Dennis Hou's entry,
Jan28DH.txt , in the contest for the fastest version of
PrimeList(n);
(for help, type: HelpDH(); );
-
Paul Raff's entry,
Jan28PR.txt , in the contest for the fastest version of
PrimeList(n);
(for help, type: HelpPF(); );
-
David John Wilson's entry
Jan28DW.txt , in the contest for the fastest version of
PrimeList(n);
(for help, type: HelpDW(); );
[Note: the lucky winner is Dennis Hou! He won five dollars!].
Homework for Thurs., Jan. 28, class (due Feb. 1, 2010)
-
(For Novices only)
Read and do all the examples, plus make up similar ones,
in pages 60-90 of Frank Garvan's awesome Maple booklet.
Hand-in (only!) two pages of random examples (that you made-up
yourself).
-
(For everyone) The Prime Number Theorem states that
the number of primes < x is asymptotic to
x/log(x) (i.e. their ratio goes to 1). Conjecture
an order-of-magnitude asymptotics for the number of
twin-prime pairs < x
-
(For experts, optional for novices):
Write a program that inputs positive integers N and k
and outputs the set of all arithmetical progressions of length k
consisting all of prime numbers ≤ N.
Using the above program, write another program (of one line)
that finds the number of such APs.
-
(For everyone)
Read (again!) carefully, and understand, pages 225-230 (sections 1 and 2) of
Norman Levinson's paper. Also do a first reading of
section 3 (p. 230-234).
Programs done on Monday, Feb. 1, 2010
-
feb1.txt
(as scribed by Emilie Hogan [thanks!]).
Contains
-
LamJ(n): the von Mangoldt function = log(p) if n=p^k and zero otherwise.
-
CheckFA(N): empirically checks the fundamental
theorem of arithmetic for the first N positive integers
in the equivalent version log(n)=sum(LamJ(i),i|n)
-
PsiJ(x): the psi(x) function (sum of log(p) over all prime powers <=x)
-
TJ(x): the sum of log(n) for n<=x (alias x!, if x is a positive integer)
-
Check210(N):Checks equation 2.10 in Levinson's paper
TJ(x)-add(PsiJ(x/n),n=1..x))=0
Homework for Monday, Feb. 1, class (due Feb. 4, 2010)
-
(For Novices only)
Read and do all the examples, plus make up similar ones,
in pages 90-120 of Frank Garvan's awesome Maple booklet.
Hand-in (only!) two pages of random examples (that you made-up
yourself).
-
(For everyone): using the package numtheory, write a program
that inputs a pos. integer n and outputs the sum of
mobius(n) over all the divisors of n. Make a conjecture about
these sums.
-
(For experts, a big challenge for novices)
[wording corrected Feb. 3, 2010, 5:48pm, thanks to Emilie Hogan]
Using the solve command, and by setting up a system
of N linear equations with the N unknowns
a[i] (i=1...N), find numbers a[i]'s that satisfy
the system of equations
{a[1]=1,seq( add(a[d], d in divisors(n))=0,n=2..N)}:
Can you conjecture an "explicit" form for the a[i]'s?
Do the a[i]'s depend on N?
-
(For experts, a big challenge for novices)
Consider all 2n subsets of {1,2, ..., n}.
Define a 2n by 2n matrix M
labelled by subsets such that M[T,S]=1 if S is a subset of T,
and 0 otherwise. Write a program that inputs a pos. integer
N and outputs the inverse of that matrix. Can you
conjecture a "closed-form" for the entries, let's call them
N[T,S].
-
(For everyone)
Define A(n) to be the sum of mobius(i) for i &le n.
Conjecture an asymptotic upper bound for |A(n)|
-
(Mandatory for experts, optional for novices):
Prove (rigorously!) the very weak inequality
|A(n)|<=100000000 n0.9999999999
-
(For everyone)
Read (again!) carefully, and understand,
section 3 (p. 230-234) of
Levinson's beautiful article.
Programs done on Thursday, Feb. 4, 2010
-
feb4.txt
Contains procedures:
-
muList(N): inputs a pos. integer N and outputs [seq(mobius(i),i=1..N)],
by using the defining property of the mobius function
(mu(1)=1, Sum(mu(i), i in divisors(n))=0 for n > 1)
-
muListD(N): inputs a pos. integer N and outputs
[seq(mu(i),i in divisors(N) ]:
-
IE(S,A): under construction, just ignore it for now.
Homework for Monday, Feb. 4, class (due Feb. 7, 2010)
-
(For everyone (but failry challenging for novices))
Recall that the mobius function μ(n) (given in the numtheory package
of Maple as mobius(n)) is defined as (-1)r if n is square-free
and a product of r distinct primes, and 0 otherwise.
Write a home-made program, mobius1(n) that looks at ifactors(n)[2],
and first decides whether it is square-free. If it is not then return 0,
if it is, then return (-1) to the power of the length of ifactors(n)[2].
-
(For everyone)
The analog of the Mobius function for sets is:
μ(S) := (-1)|S| ,
where
|S| is the number of elements of S. Write a program that inputs a set S
and outputs the sum of μ(T) where the sum ranges over all subsets T of S
(including S itself and the empty set).
you may use the function powerset(S) in combinat); . What do you get?
Can you make a conjecture?
By human means prove that conjecture.
-
Check Eq. (3.2) empirically by writing a program Check32(f,n,x), that inputs an expression f (representing a function of n),
a symbol n, and a real number x,
and outputs the difference between the left side and right side. Use the "int" command and the trunc command.
Check it out for a few f's and make sure that you get 0. For example Check32(n^3,n,30); should give 0.
-
Check Lemma 3.3., i.e. Eq. (3.7) by writing a program Check37(N) that finds the maximum of the absolute value
of the left side of (3.7) and ln(x) for x between 1 and N, and the winning x.
-
Check Lemma 3.4 empirically, by writing a program Check34(N), that
finds the largest value of |(ψ(x)-&pi(x)ln(x))/(xln(ln x)/ln x))| for x between 1 and N, and
the x where it is achieved.
-
(For everyone)
Read (again!) carefully, and understand,
everthing up to the end of section 3 (p.234) of
Levinson's beautiful article.
Programs done on Monday, Feb. 8, 2010
-
feb8.txt
Contains procedures:
-
Mob(K), the first K terms in the mobius sequence
-
M(n): The sum of mobius(i) for i between 1 and n
-
KickOutZ(L):
inputs a sequence of numbers L
and outputs the same sequence with 0's removed
(Warning: for amusement only!)
-
KOZ(L): (like KickOutZ(L) only much faster!)
inputs a sequence of numbers L
and outputs the same sequence with 0's removed
(using a suggestion of Emilie Hogan and Andrew Baxter)
-
KMob(n): The kicked-out-zero version of Mob(n)
-
RandSeq(n): a (uniform) random sequence of 1s and -1s
of length n
Homework for Monday, Feb. 8, class (due Feb. 11, 2010)
-
A random sequence of {1,-1} of length n has the properties that
-
There are roughly as many 1's as -1's
-
There are roughly as many [1,1], as [1,-1], as [-1,1], as [1,1]
-
There are roughly as many [1,1,1], as [1,1,-1], ..., as [-1,-1,-1]
Write a procedure RS1(L), that inputs a list of length n consisting of 1's and -1's and computes
-
a1:=(n/2-#1's)2+(n/2-#(-1)'s)2
-
a2:=(n/4-#[1,1]'s)2+(n/4-#[1,-1]'s)2+(n/4-#[-1,1]'s)2+(n/4-#[-1,-1]'s)2
-
a3:= (n/8-#[1,1,1]'s)2+ ....
-
Write a program RS(N,K) that repeats RS1(RandSeq(N)) K times, and finds the minimum of all the a1's,
the minimum of all the a2's, and the minimum of all the a3's.
-
Compare RS1(KOZ(Mob(N))) and RS(nops(KOZ(Mob(N)),100) for N from 100 to 100000, incremented by 100.
-
Another interesting "statistics" that a sequence (i.e. list), L[i] -1's and 1's has is the number of
times that "I was not in the red", i.e. the number of i's such that L[1]+...+L[i] > 0 or
L[1]+...+L[i]=0 and L[i]=-1. Write a program W(L) that computes this number for any such list L's of -1's and 1's.
-
Write a program that AS(N,K,x) that inputs positive integers N and K, and a variable x, and generates K RandSeq(N), computes
W(L) for each of these K random sequenecs, and adds up x^W(L).
-
Compare notes with W(KOZ(Mob(N)) for various N's.
Programs done on Thursday, Feb. 11, 2010
-
feb11.txt
Contains procedures:
-
L(n), The von Mangoldt function,log(p) if n is a p^k for some k, 0 otherwise
-
theta(x)=Sum(log(p),p<=x):
-
psi(x): Sum(L(n),n=1..x):
-
T(x)=Sum(log(i),i=1..trunc(x)); (=log trunc(x)!)
-
TestT1a(x): tests that T(x)=Sum(psi(x/i),i=1..x)
-
TestT1(x): tests that T(y)=Sum(psi(y/i),i=1..y), for y<=x
-
TestT2a(x): tests that psi(x)=Sum(mobius(k)*T(x/k),k=1..x);
-
TestT2(x): tests that psi(x)=Sum(mobius(k)*T(x/k),k=1..x) for y<=x
-
Check49a(x): The difference between the left and right sides
of (4.9) in Levinson's article with F=psi and G=T for a specific x
-
Check49(x): Checks that Check49a(y)=0 for all y<=x
Mandatory(!) Special Homework, due Feb. 14, 2010
Consider the list of lists
L:=[[8, 29, 104, 293, 680], [61, 112, 237, 496, 973], [216, 309, 504, 861, 1464], [557, 704, 989, 1472, 2237], [1192, 1405, 1800, 2437, 3400]]:
Let g=g(x,y) be the unique polynomial of degree 4 in both x and y,
such that
g(i,j)=L[i][j] 1 ≤ i,j ≤ 5
-
Use maple to find the polynomial g (as an expression in x and y, of degree 4 in x and degree 4 in y).
[Hint, write a generic polynomial in two variables x,y, of degree 4 in each, with 25 undetermined coefficients,
and use the solve command].
-
Using, xmpale, do:
plots[implicitplot](g,x=-6..6,y=-5..5,axes=none);
-
print it out.
-
Cut it out.
-
Write inside it: "I love Experimental Math".
-
Tape it on the door of Hill 704.
Homework for Thurs., Feb. 11, class (due Feb. 15, 2010)
-
Write a program Check49Abstract(x) that takes an abitrary function F(x) (left as a symbolic function!),
and defining G(x)=Sum(F(x/i),i=1..trunc(x)) (2.11), Checks Equation (4.9) of Levinson's paper for all
arguments less than x.
-
Defining R(x)=ψ(x)-x, check empirically Eq. (5.1) by writing a procedure
Check51(X) that inputs a positive integer X, and outputs the maximum of the
left side of (5.1) divided by x, amongst all x ≤ X. What is Check51(100)?
Check51(1000)?, how far can you go with a reasonable amount of time?
-
(Challenge!) Write a (numerical) procedure that inputs a real number ≥ 2 and
outputs a floating-point approximation to S(y) as defined by (5.3) (p. 238).
Hint: R(x)=psi(x)-x, and psi(x) is a piece-wise constant function that changes
every time you encounter a prime-power. Break the interval of integration [2,y] into
segments [2,3],[3,4],[4,5],[5,7],[7,8], ..., between two consecutive prime-powers
and add-up all these integrals.
-
Write a procedure Check57(x) that computes the largest ratio of the left-side
of (5.7) divided by y, for y between 2 and x.
-
Write a procedure L2(n) that computes Λ2(n) of (4.14).
What is Λ2(1000)?
-
Write a procedure Check513(x) that checks Eq. (5.13). It inputs a positive integer x, and find the largest K1
for y between 2 and x.
(i.e. the left side minus the first term of the right side divided by ylog(y)).
-
Read a first reading of the non-elementary, but nevertheless,
beautiful proof of PNT, due to to Don Newman, as told by Don Zagier
Programs done on Monday, Feb. 15, 2010
-
feb15.txt
Contains procedures:
-
Zetap(s,N):
The product of 1/(1-1/p^s) for the first N primes
in floating point, and s a complex number (in floating)
-
Zetas(s,N):
The sum of 1/n^s for n from 1 to N
in floating point, and s a complex number (in floating)
-
ACzeta(s,N): analytic continuation of PHI(s)-1/(s-1)
Homework for Monday, Feb. 15, 2010 class (due Feb. 18, 2010)
-
Write a program that inputs an expression f of t (defined for t ≥ 0) describing a bounded and locally integrable function, and
computes
g(z)=Int(f(t)*exp(-z*t),t=0..infinity)
where z is a complex number with positive Real Part given in floating point. Experiment with f(t)=1/(1+t^2) and try to plot
the real and imaginary parts of g(z).
-
Write an experimental numerical program that inputs a function f(t) and real numbers a and b (a < b), and
a positive integer K.
It outputs an approximate location of its minimum in that interval, and the value as follows.
Break-up the interval into K equal part, evaluate the function at the intermediate points
(a,a+(b-a)/K,a+2(b-a)/K,..., b) look at the minimum of these, and then between the two neighbors,
break the new (smaller) interval again, and keep doing it until you have located a very small interval
where that minimum is located.
-
Apply the above program to locate the first 10 zeroes of f(t):=abs(Zeta(1/2+I*t)).
-
Plot abs(Zeta(x+I*t)) for various x's (≠ 1/2) and t between 0 and 100, and convince yourself empirically that RH is right.
-
Read the
wikipedia article on the Fermat primality test
and translate into a Maple prorgram the pseudocode given there. Test it out.
-
Read the wikipedia article on the Miller-Rabin probabilistic primality test
and translate into a Maple prorgram the pseudocode given there. Test it out.
-
($5 award for the nicest proof!)
Find the nicest proof (combinatorial if possible) to Fermat's little theorem that ap-a is divisible by p for every a and every
prime integer p.
Programs done on Thursday, Feb. 18, 2010
-
feb18.txt
Contains procedures:
-
RE(a,n,p): computes (efficiently!) a^n mod p
-
REd(a,n,p): Like RE(a,n,p) but using Maple directly
-
IsComposite(n,a): Decides whether the pos. integer n is
composite using the witness a. If it says true, it is true
for sure. If it says false, than it is probably true.
-
IsPrime(n,k): The Miller-Rabin famous probabilstic
primality test. It inputs a pos. integer n and a pos. integer k
and performs IsComposite(n,a) k times for
k random a's. If they all say false then you
return true (with high prob.) otherwise false
Homework for Thursday, Feb. 18, 2010 (corrected Feb. 21, 20101) class (due Feb. 22, 2010)
-
Check how good is IsPrime(n,k) by writing a program
CheckInPrime(N,K,k) for N a large (but not too large) positive integer,
picks at random an integer between 1 and N, let's call it i,
and sees whether our IsPrime(i,k) agrees with Maple's built-in isprime(i)
[Added Feb. 21, 2010 (9:01am): The original version of the question stupidly asked
to look at IsPrime(ithprime(i),k) that is always true (there are no false negatives in
the Miller-Rabin test). Many thanks to Asya Pritsker for pointing this out!]
Do it for K times, and find out how many were false primes.
After you write it, try out
-
CheckInPrime(10000,10000,1);
-
CheckInPrime(10000,10000,2);
-
CheckInPrime(10000,10000,3);
(Note: if things are too slow, you can make the first two inputs smaller).
-
($5 for the best entry)
Write a Maple program that enters the date [x,y,z] in the American format Month/Day/Year
(Gregorian calendar) and outputs the day of the week (1=Sunday, ..., 6=Friday, 0=Saturday).
-
Write an algorithm for human mental math that computes the day of the week of a given date.
master it carefully. There will be a contest, and the person who can do it fastest
(in their head) for some randomly chosen dates, will win $5.
-
Look up the "Extended Euclidean algorithm", and write a program Egcd(a,b) that inputs
two positive integers a and b, and outputs the greatest common divisor, let's call it d,
and two (not necessarily positive) integers, let's call them e and f such that
e*a+f*b=d
The output should be in the form d,[e,f].
For example Ecgd(3,2); should return : 1, [1,-1].
-
Look up Euler's toitent function and write two programs ET1(a), ET2(a), both of which
compute the Euler toitent function φ(a), as follows.
- ET1(a) does it by actually counting how many integers between 1 and a are relatively prime to a, using igcd.
- ET2(a) does it by using the well-known formula featuring the prime factors.
Check that ET1(a)=ET2(a) for a between 1 and 10000.
-
Read a first reading of the
The Agrawal-Kayal-Saxena (AKS) masterpiece
-
(Riddle)
I overheard at a bookstore two people talking. One person says to the other:
look at that book, it has a number in its title. I just noticed that the product of the ages of my
three daughters equals to that number. All my daughters know how to read. Can you tell, for sure,
the ages of my daughters?
The other man replied, Yes, their ages are .... (but I didn't hear what he said), and the first man replied:
"yes", you got it!
Myself, being very smart, I figured out the ages of that man's daughters, even though I didn't hear what
the other man said. Can you?
Programs done on Monday, Feb. 22, 2010
-
feb22.txt
-
ModPol(p,x,h): inputs a polynomial p in the variable
x and another polynomial h, in x, finds p mod h
-
ord1(a,r): the smallest pos. integer k such that a^k mod r =1 (gcd(a,r)=1)
-
Check31(M) checks Lemma 3.1. in the AKS paper for m between 7 and M
-
AKS(n): the Agrawal-Kayal-Saxena famous primaity-testing algoritm (just started!, to be continued)
Homework for Monday, Feb. 22, 2010 class (due Feb. 25, 2010)
-
Finish procedure AKS(n).
-
Write a procedure FindAllPrimesAKS(M,N) that inputs pos. integers M and N ( M < N) and outputs
the list of prime integers between M and N, that uses the procedure AKS(n) that you have just
programmed.
-
Write a procedure FindAllPrimesMR(M,N,k) that inputs pos. integers M and N ( M < N) and outputs
the list of prime integers between M and N, that uses the procedure IsPrime(n,k) from
feb18.txt.
-
[Modified, Feb. 24, 2010, thanks to Kellen Myers, the previous inputs would take years]
Run FindAllPrimesAKS(1000,1100), and FindAllPrimesAKS(1000,1100,3). Did you get the
same output? How long did they each take?
Programs done on Thursday, Feb. 25, 2010
-
feb25.txt
Same as feb22.txt, but with AKS(n) completed, in other words
-
AKS(n): the Agrawal-Kayal-Saxena famous primaity-testing algoritm
Homework for Monday, Feb. 25, 2010 class (due March 1, 2010)
-
($5 for the fasted entry in Maple(!))
Try to optimize AKS(n), without "cheating", i.e. without
replacing any procedures or steps by probabilistic versions.
Call it MyAKS(n). How much faster can you have it run?
[Added March 1, 2010: congratulations to Josh Loftus for winnig the
prize!]
-
Read the first half of the proof of theorem 4.1
(up to end of page 4)
masterpiece very carefully, and write it,
by hand in your own words, spelling out everything.
-
Write a Maple program Check43(N) that empirically checks Lemma 4.3
for all n from 2 to N. Make sure that Check43(100) returns true.
-
Write a Maple program IsIntrospective(f,X,m,r,p) that inputs
a polynomial expression f in the variable X, an integer m,
and a prime number p, and outputs true iff m introspective for f(X).
Check that all the first 100 primes numbers are introspective
for X+a for all 1 ≤ a,r ≤p-1
-
Write a Maple program AKSgp(n,p,r) that inputs a composite integer
n, one of its prime divisors p, and an integer r relatively prime to
n (and hence to p), and outputs the set of elements of the group
G described in the last paragraph of page 4.
It should return FAIL, if the conditions for n,p, and r, are
not satisfied.
Programs done on Monday, March 1, 2010
-
mar1.txt
contains Josh Loftus' improved version of:
-
AKS(n): the Agrawal-Kayal-Saxena famous primaity-testing algoritm
Homework for March 1, 2010 class (due March 4, 2010)
-
Finish up the homework that you didn't do last time. It is
very important!
-
($5 for the fastest entry in Maple(!))
Try to optimize AKS(n), even on top of
what Josh Loftus did (e.g. do the exponentiation (X+a)n by
repeated squaring, and moding-out by n and by
Xr-1 at each step of the squaring.
Call it MyAKS(n).
Compare the running times, by typing
restart: read FileName: t0:=time(): {seq(MyAKS(n),n=1..10000)};time()-t0;
and
restart: read FileName: t0:=time(): {seq(AKS(n),n=1..10000)};time()-t0;
-
For a pair of two fair (cubical, i.e. 6-faced) dice, the probability
that the sum of the dots of the two faces that landed on top is 2 equals
1/36, that it is 3 equals, 2/36, that it is 4, equals 3/36,
..., that it is 11 equals 2/36, that it is 12 equals 1/36.
Construct two other, non-standard, dice (with positive number of
dots on each face, but you may have repeats) that would give the
exact same probability distribution as the above-mentioned
pair of fair dice, for the total number of dots on the top faces.
-
(A chellenge!) Generalize the above by
writing a Maple program, NSD(r,k) (for Non-Standard dice) that
inputs two positive integers r and k, and outputs all
r-tuples of non-standard k-faced dice that yield the same distribution
as r standard k-faced dice (with 1,2,3, ..., k dots on the faces,
and each face is equally likely to be on top).
Note: the previous problem was NSD(2,6);
-
One way to define the cyclotomic polynomial Cyc(n,t) is
the largest-degree monic polynomial that divides tn-1.
Write a program Cyc(n,t) that inputs a pos. integer n, and a variable
t, and outputs the n-th cylotomic polynomial, by following the
following steps
-
Use Maple's factor to find the irreducible factors of tn-1.
-
By loking at op(i,f) for i from 1 to nops(f) find the one with
the largest degree, and make it monic.
Compare your output with Maple's built-in command cyclotomic.
-
Read carefully the continuation of the proof of the validity of the AKS algorithm
(p.5). Try to understand every step, fill-in details, and rewrite it in your
own words.
Programs done on Thursday, March 4, 2010
-
mar4.txt
The programs of mar1.txt plus:
-
Mult11(p,a,x,r): p*x^a mod x^r-1 mod n
-
Mult1(p,q,x,n,r): Inputs two polynomials p and q
of degree <=r-1 and coefficients mod n
and outputs their product mod x^r-1 and mod n
-
Power1(p,x,m,r,n): p(x)^m mod (n, x^r-1)
(the "efficient" way)
-
NewAKS(n): like AKS(n), but with attempted "optimization"
(Note: bug in class corrected. In step 5 the variable
"asya" was defined incorrectly before)
-
kmAKS.txt
Kellen Myers improved version of AKS.
Homework for March 4, 2010 class (due March 8, 2010)
-
(real challenge, I don't know the answer, $5 for best entry)
Why is NewAKS(n) so much slower than AKS(n). It is supposed
to be faster! In particular, NewAKS(271) seems to take for ever.
-
(real challenge, I don't know the answer, $5 for best entry)
Resurrect NewAKS(n) to be faster than AKS(n).
-
Adapt the AKS(n) procedure to just output the integer r at
the end of step 2, that is used subsequently. Is that
sequence in Sloane?
-
According to Kellen Myers, my method of producing
the cylclotomic polynomials is not a good one.
He is probably right. Write two different good
programs, using different formulas/algorithms to
generate the cylotomic polynomial (you can consult
wikipedia or any source) and make sure that it coincides
with the output of cylotomic(n,x) in the numtheory package of maple
-
Write a Maple procedure that inputs a prime p and
an integer r and outputs a polynomial h(x) that is
an irreducible factor of cyclotomic(r,x) mod p.
Programs done on Monday, March 8, 2010
-
mar8.txt
contains the programs
-
DennisAKS.txt
Deenis Hou's AKS speed-up.
-
Gen(S,c):The multiplicative semi-group of integers modulo c
generated by the integers in S
-
IV(k,n): the set of vectors of non-negative
integers with k components that add-up to n
Homework for March 8, 2010 class (due March 11, 2010)
-
By using the methodology of Gen(S,c) or otherwise, write a Maple
program ScG(L,X,h,p) that inputs a positive integer L, a variable
X, an irreducible polynomial h in the variable X mod p, and a prime p,
and outputs the subgroup of Fp[X]/(h(X)) generated by
by X+1,X+2, ..., X+L. Test your program for small primes p and L
and h(X) an irreducible factor of cylcotomic(X,r) mod p for
various r's.
-
By mimicking the program IV(k,n), write a program,
NuIV(k,n) that computes recusively the number of elements
in IV(k,n). So replace union by adding etc.
Check the program by verifying that NuIV(k,n)=nops(IV(k,n))
for k between 0 and 8 and n between 0 and 10.
-
By using experimental mathematics, conjecture
a closed-form expression for NuIV(k,n).
-
Prove your conjecture (by induction or otherwise)
-
(challenge, $5 to be divided). Find a combinatorial proof
by establishing a bijection with a well-known problem
(i.e. an even more well-known one).
[Added March 11, 2010: Aron Samkoff and Kellen Myers shared the prize]
Programs done on Monday, March 8, 2010
Today we mainly finished up understanding the beautiful AKS proof.
The little programming we did is contained in:
-
mar11.txt
contains the programs
-
GC(n): the set of prime pairs [p,q], with p < q such p+q=2n
-
TestGC(N): The number of solutions of p+q=2*N divided by N/(log N)^2
Homework for March 11, 2010 class (due March 22, 2010)
-
Write a program AlmostPrimes(N,k) that inputs an integer N, and outputs all the integers less than or
equal to N that are the product of at most k primes.
-
Determine experimentally (empirically) an asymptotic formula for the number of almost k-primes
< N
-
Determine experimentally (empirically) an asymptotic formula (in n) for the number of
solutions of a+b=2n, where a and b are almost k-primes
-
Read, as carefully as you can
Viggo Brun's classic paper on the Goldbach problem
-
Write a Maple program N(Delta,D,x,aList,pList) that implelents
the function in p.101 of Viggo Brun's paper.
-
Have a good Spring Break. See you Monday, March 22, 2010.
Programs done on Monday, March 22, 2010
Today we mainly finished up understanding the beautiful AKS proof.
The little programming we did is contained in:
-
mar22.txt:
The Principle of Inclusion-Exclusion, contains procedures
-
PIE(S,P): inputs a set of objects S and
a list of "properties" P=[S1,S2, ..., Sk]
where S1 is the set of objects enjoying the first property ,
S2 is the set of objects enjoying the second property ,
...
Sk is the set of objects enjoying the k-th property ,
checks the Principle
of Inclusion Exclusion
|(1-S1)(1-S2)...(1-Sk)|=|S|-|S1|-|S2|-...
+|S1 intersect S2|+....
by computing it via the formula, and directly by
nops(S\(S1 union S2 union ... Sk))
-
RandSS(S): generates a random subset of S
where each element has prob. 1/2
-
TestPIE(N,M,k): Tests PIE for random N-element set with pos. integers ≤ M
and k properties
-
Bonn(S,P,r): checks level-r Bonferroni inequality
against the true value, yielding successive upper bounds and lower bounds
for the number of elements with no properties
Homework for March 22, 2010 class (due March 25, 2010)
-
A permutation pi of {1,...n} is a derangement if for every i between 1 and n,
pi[i] ≠ i . Write a program Der(n) that computes
the number of derangements in two different ways using PIE(S,P)
by defining S to be convert(permute(n),set) and P the list of sets of permutations
whose i-th element is the set of permutations with pi[i]=i.
(Note: this is a most inefficient way of doing it, but do it anyway).
-
By looking at PIE, derive a quick formula for Der(n).
-
What is the limit of Der(n)/n! ?
-
Read carefully
my three-page masterpiece
-
Write a program, Bonfferoni(S,P,k), that checks the identity version of the Bonfferoni formula
as descibed in Example 2 on p. 2.
(Hint: The left-over is a sum over subsets of k-elements {i1, i2, ..., ik}
of (-1)k times |S[i1] intersect S[i2]...intersect S[ik] intersect (S minus S[ik+1]) intersect (S minus S[ik+2])... intersect (S minus S[nops(P)])
-
For a set T={i1,i2, ... ik}, let xT denote x[i1]* ...*x[ik].
The Bonferonni identity follows from the algebraic identity
(1-x1)(1-x2) ...(1-xn)=Sum((-1)|T|xT, |T| < k)+
Sum((-1)|T|xT(1-x[max(T)+1])(1-x[max(T)+2])...(1-x[n)), |T|= k) .
Write a program Bo(n) that tests this algebraic identity empirically for a pos. integer n.
Programs done on Thurs, March 25, 2010
-
mar25.txt:
Starting the Brun game. Contains procedures
-
ER(x,r): the lower bound for the number of
integers <=x that are not divisible by the
first r primes, using the highly inefficient sieve
of Eratosthenes.
-
ERtest(x,r): the exact cound for the
number of integers <=x that are not divisible by the
first r primes, done directly, by constructing
the set in question.
-
AllWalls(S,r): the collection of all subsets
of {1, ..., r} that contain S
AllWals({1,3},4) should be e
{{1,3},{1,2,3},{1,3,4},{1,2,3,4}}
-
BrunGame1(x,r,S): Not yet complete!
(in fact, just started, to be completed next time)
It should input an integer x , a pos. integer r, and a set of
subsets of {1,...,r} such that we know that
NumberOfIntegers x not divisible by the
first r prime >x*Sum((-1)^|S1|/mul(ithrpime(i), i in S1),S1 in S)-nops(S)
and tries to remove free-loading wallss with "apex"
of even cardinality
Homework for March 25, 2010 class (due March 29, 2010)
-
($5) Find the factual mistake in
my masterpiece (Hint: it is in the Speculations section).
-
An implicant, denoted by [S,T], where S is a subset of T and T is a subset of {1, ...,n},
is the collection of all subsets of {1,...,n} that contain S and is a subset of T
(i.e. the "interval" in the Boolean lattice between S and T). For example
[{3,4},{3,4,7,8}] consists of the four sets {{3,4},{3,4,7},{3,4,8},{3,4,7,8}}.
Write a Maple program: Implicant(S,T) that outputs that interval.
-
(challenge!)
A Boolean function, f, on n variables can be viewed as a collection of subsets of {1,2,..,n}.
A prime-implicant (not to be confused with prime number)
is an implicant that is containded in f, and is not contained in
a larger implicant. Write a program PrimeImplicants(f,v,n) that inputs a "Boolean function" f
(i.e. a collection of subsets of {1,...,n}, a member v in f, and outputs all
prime-implicants containing v. (Note: this is always non-empty, but may-be the singelton {v}, in case
v is not part of a larger implicant.
-
(challenge!)
Write a program, QM(f,n) that "simplifies" a Boolean function f (given as above) and outputs
a collection of prime-implicants whose union is f. For example
SB({{1},{2},{1,2}}) should return {[{1},{1,2}],[{2},{1,2}]}.
You should do it "greedily" by first picking members of f that have only one prime-implicant,
and removing the members of that prime -mplicant, getting a smaller f.
When you run out of such members of f, pick members of f that have a unique largest prime-implicant,
and choose that prime-implicant in the constructed collection of prime-implicants, and remove all their members.
Finally, when you are left with the reduced f, where each member has several maximal-size prime-implicants,
choose any at random, until you cover everything. If you succeeded you (essentially) programmed the
Quine-McClusky algorithm. It is equivalent to expressing a set of vertices of the n-dimensional
unit cube as a union of "walls" of various dimensions, as economically as possible.
Programs done on Thurs, March 29, 2010
-
mar29.txt:
Contains procedures
-
BonnSieve(n,k): The Bonnferoni sieve of level k with n properties
-
ApplySieve(Se,P,Si): applies the sieve Si to estimate the number of elemenets
in Se minus the union of the P[i]'s
Special Homework for March 29, 2010 class (due March 30, 2010)
Use MathSciNet to compile a set of lists of length 4 of the format
{[FirstnameLastname,ErdosNumber,ZeilbergerNumber, ShelahNumber]}
for the list of names that you were given in class.
Usual Homework for March 29, 2010 class (due April 5, 2010)
-
Using the description in
my masterpiece, write a procedure BrunSieve to construct
the BrunSieve(k,N) first naively, and then, if possible cleverly.
April 1,2010
We had a visiting lecture by the famous Professor Andrew Baxter.
Prof. Baxter's Homework for April 1, 2010 class (due April 8, 2010)
- Write a procedure that uses tubeplot to plot the (p,q)-torus knot for a
given p and q, using the parameterization given at:
wikpedia's page on Torus knot
- Use textplot and the DrawYoung program we wrote in class to draw Young
Tableaux, as illustrated at
wikipedia's page on Young tableaux
- Use the DrawYoung program as a guide to write a procedure to draw 3-D
Young Diagrams for Plane Partitions (use the cube command in plottools).
See an example here.
Programs done on Monday, April 5, 2010
-
apr5.txt:
Contains the procedures of mar29.txt and
-
MakeSeP(n,k):inputs pos. integers n and k
and outputs the set {1,...,n} followed by
a list of k random subsets of {1,...n}
each with a random number of elements
-
BrunSieve(k,N) inputs the number of properties k
and a list N of weakly decrreasing pos. integers
N=[N1,N2,N3, Nr] with r<=k
outputs all subsets of {1, ..., k},
with exactly r elements written as
sets {i1,i2,..., ir} such that k>=i1>i2> ...>ir>=1 and i1<=N1 and i2<=N2, ...ir<=Nr
Homework for April 5, 2010 class (due April 8, 2010)
-
(Really easy!)
The procedure we did in class, BrunSieve(k,N), with k a pos. integer (number of properties)
and N a weakly decreasing set of integers N=[N1, ..., Nr]
only outputs the subsets of {1,...,k}, {i1,i2, ..., ir}, with k ≥ i1 > i2 > ... > ir ≥ 1 with
i1 ≤ N1, i2 ≤ N2, ..., ir ≤ Nr.
Using this as a subroutine, write a Maple program FullBrunSieve(k,N) that outputs all subsets
(including the crucial empty set) with ≤ r elements, such that if it is written {i[1],i[2], ..., i[s]}
with s ≤ r and k ≥ i[1] > i[2] > i[3] > ... > i[s] ≥ 1 we have
i[1] ≤ N1, ..., i[s] ≤ Ns .
-
Experiment with ApplySieve and convince yourself that if
N1=N2,N3=N4, ..., N[n-1]=N[n] (n even) you get an upper sieve
and with N2=N3,N3=N4, ..., N[n-1]=N[n] (n odd) you get a lower sieve.
-
Write a procedure Check12(x,r) that inputs a real positive number x, and
a pos. integer r, and computes both sides of Eq. (12) on p.114
of
Viggo Brun's classic paper with D=1, and check empirically for all
x &le 1000 and all r ≤ 10 that it is always true (provided that the condition
m+2 > σ=1/p1+...+1/pr holds ).
Programs done on Thursday, April 8, 2010
-
apr8.txt:
Contains the procedures of apr5.txt and
-
FullBrunSieve(k,N): the Full Brun Sieve with k properties and a list of weakly decreasing pos. integers N
-
HMPB(x,y): pi(y-1)-pi(x-1):end:
-
BrunSC(r,n,a): inputs a pos. integer r, a pos. integer n, and a real number a greater than 1, outputs
the Brun Staircase
Homework for April 8, 2010 class (due April 12, 2010)
-
Define the weight of a set of integers {i1,i2,i3, ...} to be x[i1]*x[i2]*x[i3]*...,
for example the weight of the set {1,5,7} is x[1]*x[5]*x[7].
Define the weight of a set of sets to be the sums of the weights of its members. For
example, the weight of {{1,2},{2,3}} is x[1]*x[2]+x[2]*x[3].
The elementary symmetric function ek(x1, ..., xn) may be defined
as the weight of the collection of all k-element subsets of {1, ...,n}.
Write a program ek(k,n,x) that inputs pos. integers k and n with 0 ≤ k ≤ n and
a symbol (letter) x, and outputs the elementary symmetric function
ek(x1, ..., xn).
By doing
- Naively, using choose(n,k)
- cleverly, by doing option remember and using the recurrence
ek(x1, ..., xn)=ek(x1, ..., x(n-1))+xn*ek-1(x1, ..., x(n-1))
(Explain why that recurrence is true)
Make sure that you get the same output
-
Write a procedure ekB(n,k,x,N) that inputs n,k,x as before, but also a list N of size k of decreasing pos. integers
N=[n1,n2, ..., nk], and outputs the weight enumerator of all k-element-subsets of {1,...,n}, let's call them
{i1,i2, ..., ik} with i1 > i2 > ... > ik and i1 ≤ n1, ..., ik ≤ nk
By doing it
- Using BrunSieve
- Cleverly, by discovering an appropriate recurrence analogous to the above.
Make sure that you get the same output.
-
Read
Dr. Z.'s short intro to Combinatorial Games
(We will take a break from Primes and do some more fun stuff).
Programs done on Monday, April 12, 2010
Today we investigated Eric Angelini's
glass of worms problem that inspired Sloane's sequence
A151986.
-
apr12.txt
constains procedures:
-
V1(L):Inputs a sequence of non-neg. integers
and performs one step n Eric Angelini's
(see A151986)'s glass worms game
-
IsOld(L): inputs a list and outputs
true if the last entry already showed up before
followed by the cycle. For example
IsOld([3,5,1,2,1]); reuturs true,[1,2,1] .
-
NuM(L): inputs a sequence of non-neg. integers and
outputs either the CYCLE, NumberOfMovesToGetThere
or FAIL, NumberOfMovesToGetThere
-
Comps(n,k): the set of lists of length k consisting of non-neg.
integers that add-up to n.
-
Enjoy
Kellen Myers neat worm page .
Homework for April 12, 2010 class (due April 15, 2010)
-
Write a program MyWormSeq(N) that inputs a pos. integer N and outputs
the sequence "the maximum number of steps needed to reach a cycle
for worms with n glasses and n worms (using Comps(n,n)).
-
(challenge, I don't know the answer, $5)
Modify V1 and/or the count so that things would agree with
Sloane's
A151986.
[Added April 15, 2010: both Kelln Myers and Dennis Hou shared the prize]
-
Write a program WorstWorm that inputs a pos. inetger n and outputs
the finite-sequence of lengh n that adds up to n that yields the
maximal number of steps until detecting a cycle.
-
(challenge, I don't know the answer, $5)
Modify V1 and/or the count so that things would agree with
Sloane's
A151987.
[Added April 15, 2010: both Kelln Myers and Dennis Hou shared the prize]
(Read
Dennis Hou's code)
-
Write a Maple program that inputs a string (list) consisting of 2's and 3's
and outputs its continuation according to the rule that the next term is
the curling number of the current sequence in the sense of
A116909.
-
Test the curling number conjecture for all strings of length ≤ 15 of 2's and 3's
-
($100) Prove the curling number conjecture .
Programs done on Thursday, April 15, 2010
Today we started Combinatorial Games!.
-
apr15.txt
contains the following procedures
(Note: I finished-up PRg(n,k) ).
-
Orbit1(G,i): the orbit of i in the directed graph G
-
IsCyclic(G): inputs a directed graph G in the
format ListOfSets and decides whether there exist cycles
where it is assumed that the vertices are called
1,2, ..., nops(ListOfSets) and ListOfSets[i] is the set of
vertices reachable from i in one more of the game.
It outputs true if the directed graph contains a cycle
(so it is not good for a game-graph, it has "perpetual check")
and false otherwise.
-
Losers(G):
it inputs a directed graph in the above format and outputs
the set of "losing vertices", in the artificial notation
of calling the vertices (alias game-positions) 1,2,3, ..., nops(G)
-
PRg(n,k): the game of removing up to k pennies
from a pile of n pennies. The output is the directed graph
in our data-structure for directed graphs G,
where G is a list of subsets of {1, ..., n} (where n=nops(G))
evertices are assumed to have names 1 through n
and G[i] is the set of vertices reachable from vertex i in one move
followed by the "key" the list of vertices in their "natural" names.
For example, try: PRg(10,3);
Homework for April 15, 2010 class (due April 19, 2010)
-
The new data structure for directed graphs that come from
combinatorial games would be a pair G,V where G is
a directed graph as we did before with the vertices
labelled by 1,2,3,..., nops(G) and V is a list of
"natural names", where V[i] is the natural name of
the vertex numbered i. To get used to it, look at
the output (for example) of PRg(10,3).
Write a procedure NaturalLosers(G,V) that inputs
a directed graph G, a list V of length nops(G),
consisting of the natural names, and returns the set
of losers using their natural names, rather than their
artificial "social-security" names.
-
The classical game of 2-pile Nim consists of two piles
having N1, N2 pennies. A legal move consists of
removing any number of pennies from ONE of the piles.
Write a Maple procedure Nim2(N1,N2) that inputs two
non-neg. integers and outputs the pair G,V where G is a directed
graph in the condensed form and V is the list of natural names
(that are the positions [i1,i2] with
0 ≤ i1 ≤ N1, 0 ≤ i2 ≤ N2
-
Usings NaturalLoses(G,V) try to conjecture the set of losing positions.
Using human reasoning (or otherwise) prove your conjecture.
-
The classical game of 3-pile Nim consists of three piles
having N1, N2, N3 pennies. A legal move consists of
removing any number of pennies from ONE of the piles.
Write a Maple procedure Nim3(N1,N2,N3) that inputs three
non neg. integers and outputs the pair G,V where G is a directed
graph in the condensed form and V is the list of natural names
(that are the positions [i1,i2,i3] with
0 ≤ i1 ≤ N1, 0 ≤ i2 ≤ N2, , 0 ≤ i3 ≤ N3
-
($5 w/o looking-it-up) Usings NaturalLoses(G,V) try to conjecture the set of losing positions.
-
Wythoff's game is like 2-pile Nim with the extra move of removing the SAME number of pennies
from both piles. Write a Maple procedure Wyt(N1,N2) that outputs the directed graph G and the
"key" V.
-
($5 w/o looking-it-up) Usings NaturalLoses(G,V) try to conjecture the set of losing positions.
Programs done on Monday, April 19, 2010
-
apr19.txt
Contains, in addition to the procedures of apr15.txt,
-
Losers(G,K): Given a directed graph G
representing a comb. game in condensened notation
where the vertices are labelled 1,2, ..., nops(G)
and a list K, where K[i] is the natural name
of vertex numbered i,
outputs the losers with their natural names
[Note: this is what we called NaturalLosers before]
-
Nim3(N1,N2,N3): the directed graph for three-pile Nim
where the piles have size ≤ N1 and ≤ N2 pennies
-
Grab3i(S,i): given a set of triplets looks
for those that contain an i in the last position
and returns the set w/o that i
Homework for April 19, 2010 class (due April 22, 2010)
-
Using Nim3(N1,N2,N3) and Losers
Formulate conjectures for the losing positions [a,b,i] with 3-pile Nim, if
i=1, i=2, i=3, i=4.
-
The winning status of a triple [a,b,c] is obviously symmetric with respect to a,b,c.
So it is enough to list the losing positions [a,b,c] with a ≥ b ≥ c. Write a
program Nim3F(N) that inputs a pos. integer N and outputs the directed graph whose vertices are
[a,b,c] with N ≥ a ≥ b ≥ c ≥ 0 and such that the outgoing neighbors of
[a,b,c] are the sorted vectors that go out of it in one move.
-
Recall that Wythoff's game is like 2-pile Nim with the extra move of removing the SAME number of pennies
from both piles.
Last time I asked you to
Write a Maple procedure Wyt(N1,N2) that outputs the directed graph G and the
"key" V. We also have symmetry here
Write a program NewSyt(N) that inputs a pos. integer N and outputs the
set of losing positions [a,b] with N ≥ a ≥ b ≥ 0, arranged in increasing
order. The output should start [[0,0],[2,1],[5,3], ...],
Let's call the k-th pair [A[k],B[k]]. Conjecture what A[k]-B[k] equals to.
-
Conjecture what the limit A[k]/B[k] equals to as k goes to infinity.
-
($5, w/o looking it up): prove your conjecture.
-
Write a program StupidGame(N,S),
that inputs a pos. integer N and a set of pos. integers S, and
outputs the directed graph of the game
where there is a pile of at most N pennies
and a legal move consists of removing a number of pennies in S.
For example, PRg(n,k) of last time is the same
as StupidGame(n,{1,2,..,k});
-
Find the losing positions of StupidGame(30,S) for various sets S.
Can you detect a pattern?
Programs done on Thursday, April 22, 2010
Today was a beautiful day, so we did experimental math
on the grass using our very primitive computers between
our shoulders. Here is what we did with pencil and
paper.
-
apr22.txt
-
mex(S): intputs a set of non-neg. integers and outputs the smallest
integer that is not in S
-
g(G): given a directed graph G in the data structure
where the set of vertices is {1, ..., nops(G)} and
G[i] is the set of vertices reachable from i in one move
outputs the list whose i-th entry is the value of the
value Grundy function at vertex i
-
g2(i,j): The Grundy function value of 2-pile Nim
where there are i pennies in the first pile and
j pennies in the second pile.
Homework for April 22, 2010 class (due April 26, 2010)
-
Prove, by human means (by induction) that
g2(2*i,2*j)=2*g2(i,j) ,
g2(2*i+1,2*j)=2*g2(i,j)+1 ,
g2(2*i,2*j+1)=2*g2(i,j)+1 ,
g2(2*i+1,2*j+1)=2*g2(i,j) .
-
To get the the Nim-sum of non.neg. integers i and j you convert them each to binary,
add them up without carry, and convert the answer back to the usual. For example
NimSum(3,5)=011+101=110=6.
Prove by induction, using the above recurrences that you have already proved,
that g2(i,j)=NimSum(i,j).
-
Recall that the Grundy function value of a directed graph without cycles is defined,
recursively by the property
g(v)=mex({g(v1), v1 in N(v)})
Given two different (or two copies of the same) games G1 and G2, the direct product of the games G1xG2 consists
of all pairs (v1,v2) and
(v1,v2)->(v1',v2) if v1->v1' in G1 and (v1,v2)->(v1,v2') if v2->v2's in G2
Prove that g(v1,v2)=NimSum(g(v1),g(v2))
-
Write a program g2W(i,j) that computes the Grundy function value of position (i,j) in Wythoff's game
(2-pile Nim with an extra option of removing the SAME number of pennies from both piles
-
A position is a winning position iff g2W(i,j)=0. Write a program
WinningPositions(n) that inputs a pos. integer N and outputs the first n winning pairs
[[A1,B1],[A2,B2], ...., [An,Bn]]
-
Using that data
-
Make conjecture about Bn-An
-
($5) Make a conjecture about the limit of An/n (without looking it up)
-
($10) Prove that conjecture (without looking it up)
-
Write a program g3W(i,j,k) that computes the Grundy function value of position (i,j,k) in generalized Wythoff's game
(3-pile Nim with an extra option of removing the SAME number of pennies from all three piles)
Programs done on Monday, April 26, 2010
-
apr26.txt
(transcribed by Emilie Hogan). Contains procedures
-
g2F(i,j): computes the nim sum of i and j
-
Wyt(i,j): computes the grundy function of Wytoff's game
with first pile of size i and second pile of size j
-
CD(N): inputs a positive integer N and outputs a list L such
that L[k+1] is the list of positions with Grundy function
exactly k, for all positions [i,j] with i+j<=N
-
Wyt0(m): the first m losing positions in Wythoff's game using
the formula [seq([trunc(i*phi),trunc(i*phi+i)],i=0..m)]
involing the Golden ration phi .
Homework for April 26, 2010 class (due April 29, 2010)
-
Write a program GenWyt(i,j,S) to compute the values of the Grundy function of
a generalized Wythoff game where one can remove any number of pennies
for one pile (like in 2-pile Nim) and remove k1 pennies from one file and k2 pennies
from another pile provided that |k1-k2| belongs to S.
Note that GenWyt(i,j,{}) is Nim2 and GenWyt(i,j,{0}) is Wythoff.
-
Challenge (I don't know the answer, $10). Find other cases of S where there is a "nice" description of
the losing positions.
-
Write a program GenWyt3(i,j,k,m) to compute the values of the Grundy function of
a generalized 3-pile Wythoff game with piles with i,j, and k pennies,
where one can remove any number of pennies
for one pile (like in 2-pile Nim) and remove k1 pennies from one file and k2 pennies
from another pile and k3 pennies from yet another pile provided that
k1+k2+k3 mod m =0.
Programs done on Thursday, April 29, 2010
Today we had a guest-lecturer,
Dr.
Christoph Koutschan, who taught us Mathematica and
the amazing
Buchberger Algorithm
-
apr29.nb
(Mathematica Notebook done in class)
Homework for April 29, 2010 class (due May 3, 2010)
-
Read and study, and experiment with
Dr.
Koutschan's Mathematica Notebook that implements the
Buchberger algorithm with the lexicographic order.
-
Write a Mathematica procedure Sym[n,x] that inputs a symbol x
and a pos. integer n, and outputs the Groebner bases of
the ideal generated by the polynomials
x[1]+...+x[n]
x[1]2+...+x[n]2
x[1]3+...+x[n]3
...
x[1]n+...+x[n]n
-
Make a conjecture about the form of this Groebner
basis for arbitrary n.
-
Prepare the first draft of your proect, due, 11:59:59pm, May 5, 2010.
Programs done on Monday, May 3, 2010
Today, the last day of classes, we investigated
an intriguing conjecture narrated to me, last Friday
by Professor
Roger Nussbaum
and Dr.
Bas Lemmens,
see
Professor Nussbaum's message.
The conjecture that we investigated, that while a special case
of the general 2n conjecture of Nussbaum,
seems to be indicative of the general problem, and is also
still open. It says that any Boolean transformation,
when OR is replaced by max and AND is replaced by min,
when applied to the box {0,1/2,1}n with 3n
vertices, the largest cycle has still length 2n.
-
may3.txt.
Contains procedures
-
randB(n,m,k,x):a random Boolean function in disjunctive
normal form ("sum" of "products"} of Boolean literals
in the variables x[1], ..., x[n]
n=number of variables
m is the number of clauses
k is the number of terms in each clause
with the data structre of the form (e.g.)
{{x[1],1-x[3],x[4]},{ x[2],1-x[3],x[5]}}
-
Ev(f,n,x,v): evaluates a Boolean function f
in x[1], ..., x[n] with the above data structre
at the point x[1]=v[1], ..., x[n]=v[n]
Since we use min and max, it makes sense for all inputs,
not just 0-1 vectors
-
randT(n,m,k,x): a random Boolean transformation, i.e.
a list of n Boolean functions
-
EvT(T,n,x,v): evaluates a Boolean transformation
in x[1], ..., x[n] with the above data structre
at the point x[1]=v[1], ..., x[n]=v[n]
-
Cycle(T,x,n,v):
The first cycle encountered in the orbit of
the transformation T starting at v
-
MaxCyc(n,m,k,K): finds the largest length of
orbits of Boolean transformations applied to random
starting points in {0,1/2,1}^n K times, followed
by the winner.
Homework for May 3, 2010 class (due May 3, 2110)
-
(5$) Prove that 2n is sharp (won by Kellen Myers,
with a brilliant proof: find the Boolean transformation that
implements mapping [x1,...,xn] (viewed as an integer written
in base 2) to "[x1,...,xn]+1" .
-
Prove ($100) or disprove ($25) the Nussbaum 2n
conjecture specialized to Boolean transformations
acting on {0,1/2,1}n
-
Don't forget to Email me your project, by May 5, 2010,
23:59:59, or else you will fail this class, no exceptions!
I prefer that you will put the file in your own website,
and just Email me the url.
-
Have a great summer.
Here are the students'
Final Projects,
that will be continually updated (indefinitely). Experimental Mathematics never ends!
Added May 17, 2010:
Here are the students'
evaluations.
Dr. Z.'s teaching page