EXPERIMENTAL MATHEMATICS 2022
P. ZINN-JUSTIN
Earlier versions by J. de Gier, A. Garbali, A. Ghitza
1
2 P. ZINN-JUSTIN
1. Introduction
1.1. Goals. The aim of this course is to provide insights in some of the tools that are
available to modern researchers for developing new mathematics. Mathematical results p-
resented in textbooks follow from logical chains of arguments. In practice, however, very
often an important new theorem is found after extensively “playing around” with a problem,
conjecturing the answer and then finding a logical chain of arguments.
This way of making progress was already used by Gauss, but in more modern times com-
puters have played a large role in solving problems that would otherwise be intractable.
This course aims to develop students’ mathematical inventiveness and imagination in math-
ematical experimentation using modern tools such as symbolic manipulation packages. More
precisely the experimentation with the use of symbolic packages (or any programming soft-
ware) allows us to
(1) Gain insights about the problem
(2) Discover patterns and relationships
(3) Use graphical visualisation
(4) Test conjectures
(5) Explore a result to see if it worth a formal proof
(6) Suggest proofs
(7) Replace lengthy derivations with computer based derivations
(8) Confirm results
1.2. The programming languages. In this subject you will be required to write programs
in either Macaulay2 or Mathematica. This subject is however not a programming course,
and you will have to develop basic programming skills yourself with practice, either at
home or during our computer lab sessions. Fortunately, learning how to use Macaulay2 or
Mathematica is not hard, and you will have plenty of time to get familiar with the software
and the languages. You can run Mathematica natively on your own computer, or on the
lab computers; Macaulay2 is best run in the browser at https://www.unimelb-macaulay2.
cloud.edu.au/. You are encouraged to turn to the documentation whenever you need to
work out how to do something or what a particular function does: for Mathematica you can
go to
https://www.wolfram.com/language/fast-introduction-for-math-students/en/
whereas for Macaulay2 you can browse the documentation at
https://www.unimelb-macaulay2.cloud.edu.au/usr/share/doc/Macaulay2/Macaulay2Doc/
html/index.html
Macaulay2 and Mathematica are very similar languages; Here we briefly describe some of
the common features and differences:
• Both are interpreted languages, which means you have an interactive session: in
Macaulay2, you can type
i1 : 6*7
o1 = 42
and obtain the answer immediately; and the same in Mathematica:
6*7
• Both are procedural programming languages: in Macaulay2
i2 : for i from 1 to 5 do print i^2
EXPERIMENTAL MATHEMATICS 2022 3
1
4
9
16
25
vs in Mathematica
For[i = 1, i <= 5, i++, Print[i^2]]
• They both have a modicum of functional programming:
i3 : apply({1,7,12}, x -> x^2)
o3 = {1, 49, 144}
o3 : List
vs
Map[Function[x, x^2], {1, 7, 12}]
Note that where Macaulay2 uses (), Mathematica uses [].
• Macaulay2, being a more modern language, also has class inheritance: (often associ-
ated to object-oriented programming)
i4 : f = method();
i5 : f Number := x -> print "this is a number (but not a real number)";
i6 : f RR := x -> print "this is a real number";
i7 : f 3.14 -- this calls f RR
this is a real number
i8 : f (5+ii) -- this calls f Number
this is a number (but not a real number)
Note the use of method to define a function whose behaviour depends on the class of
its argument (here, RR is a descendant of Number), of ; to suppress output, and of
-- for comments.
4 P. ZINN-JUSTIN
2. Integer sequences
In this section we will study some problems where involving integer sequences (an)n≥0;
in particular, one would like to know a “closed form” formula for an. This will give us an
excuse to learn our programming languages.
2.1. Guessing integer sequences. A valuable online resource is the The On-Line Ency-
clopedia of Integer Sequences at
https://oeis.org/
This database contains many sequences, and information about sequences, that have oc-
curred in the mathematics literature. You can search for any (finite) sequence, and if it
happens to be (a subsequence of a sequence) in the OEIS, you will receive information about
it.
Example 1. Consider the sequence (an)n≥1,
an = 1, 2, 6, 24, 120, 720, 5040, . . .
If you haven’t already recognised this sequence, you might try to recognise the sequence
(bn)n≥1 of ratios bn = an+1/an,
bn = 2, 3, 4, 5, 6, 7, . . .
Guessing that bn = n+ 1, we discover that an+1 = (n+ 1)an, and thus that an = n!.
Example 2. A very classical problem in enumerative combinatorics is to count triangulations
of an n-gon:
It leads to the sequence
cn = 1, 2, 5, 14, 42, 132, 429, 1430, . . .
This may also look familiar, but let’s consider again the sequence of ratios
bn = 2,
5
2
,
14
5
, 3,
22
7
,
13
4
,
10
3
, . . .
which you may be able to recognise as
bn =
2 + 4n
2 + n
, n = 1, 2, . . .
It follows that the original sequence were the Catalan numbers
cn =
1
n+ 1
(
2n
n
)
.
There are plenty more interpretations of these numbers, see https://oeis.org/A000108.
One of their first appearances is in the work of Minggatu (early 18th century), who found
the curious formula
2 sinα− sin(2α) =
∞∑
n=0
cn
4n
sin2n+3 α
They were also investigated by Euler, and a hundred years later by Catalan!
EXPERIMENTAL MATHEMATICS 2022 5
What makes these examples easy is that the ratios bn are simple rational expressions of
the form
(1) bn =
αn+ β
α′n+ β′
.
Let’s look at another example:
Example 3. Consider the sequence
(2) an = 1, 2, 7, 42, 429, 7436, 218348, 10850216, 911835460, . . .
These numbers count n × n alternating sign matrices, which we will encounter in the
tutorial. Trying again to recognise the sequence of ratios bn = an+1/an
bn = 2,
7
2
, 6,
143
14
,
52
3
,
323
11
,
646
13
, . . . ,
by fitting it to the form (1) does not lead to success.
Here, a variety of strategies can be employed:
• One can try to factor the numbers an:
i1 : L={1,2,7,42,429,7436,218348,10850216, 911835460};
i2 : apply(L,factor)
o2 = {1, 2, 7, 2 · 3 · 7, 3 · 11 · 13, 2211 · 132, 2213217 · 19, 2313 · 172192, 225 · 17219323}
o2 : List
or in Mathematica
FactorInteger[{1,2,7,42,429,7436,218348,10850216, 911835460}]
The good news is that despite being very big, these numbers have small prime factors!
This suggests that there should exist a “nice” product formula for them. But which?
• One can try to push further the idea of taking ratios. How about defining cn =
bn+1/bn, so that
cn =
7
4
,
12
7
,
143
84
,
56
33
,
969
572
,
22
13
,
115
68
. . .
While this is still not of the form (1), it turns out that cn is a ratio of two quadratics
in n:
(3) cn =
αn2 + βn+ γ
α′n2 + β′n+ γ′
.
Exercise 1. Find the parameters in (3).
• Finally, one can cheat and ask OEIS! In Macaulay2, this is easily implemented as:
i3 : oeis L
o3 = {0. A005130 Robbins numbers: a(n) = Product_{k=0..n-1} (3k+1)!/(n+k)!;
also the number of descending plane partitions whose parts do not exceed n;
also the number of n X n alternating sign matrices (ASM’s).}
o3 : NumberedVerticalList
which conveniently gives us a product formula for an:
(4) an =
n−1∏
k=0
(3k + 1)!
(n+ k)!
.
6 P. ZINN-JUSTIN
2.2. The Collatz conjecture. A different type of integer sequences occurs in the so-called
Collatz conjecture, which asks about the following procedure: given an integer n, divide it
by 2 if it’s even, or else multiply by 3 and add 1. If we repeat this over and over, does the
process always reach 1?
We can experiment by computer: first let’s define one step of the procedure:
i4 : collatz1 = n -> if even n then n//2 else 3*n+1
o4 = collatz1
o4 : FunctionClosure
Next we can ask for the whole sequence until it hits 1: (which it should eventually, if the
conjecture is true!)
i5 : collatz = n -> while n != 1 list n = collatz1 n
o5 = collatz
o5 : FunctionClosure
In Mathematica, it would look like
collatz1[n_] := If[EvenQ[n], n/2, 3 n + 1];
collatz[n_] := NestWhileList[collatz1, n, # != 1 &]
Let’s try:
i6 : collatz(12345)
o6 = {37 036, 18 518, 9 259, 27 778, 13 889, 41 668, 20 834, 10 417, 31 252, 15 626, 7 813,
23 440, 11 720, 5 860, 2 930, 1 465, 4 396, 2 198, 1 099, 3 298, 1 649, 4 948, 2 474, 1 237, 3 712, 1 856,
928, 464, 232, 116, 58, 29, 88, 44, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1}
o6 : List
It stopped, conjecture checked!
In order to understand a process it is often useful to tabulate the results of applying it
many times. One feature of the Collatz process is how many steps it takes to get to 1. We
can tabulate this statistic for the first 25 values of n with the function tally, as follows:
i7 : T = tally for n from 1 to 30 list length collatz n
o7 = Tally{0 ⇒ 1, 1 ⇒ 1, 2 ⇒ 1, 3 ⇒ 1, 4 ⇒ 1, 5 ⇒ 1, 6 ⇒ 1, 7 ⇒ 3,
8 ⇒ 1, 9 ⇒ 2, 10 ⇒ 2, 12 ⇒ 1, 14 ⇒ 1, 15 ⇒ 2, 16 ⇒ 1, 17 ⇒ 2,
18 ⇒ 3, 19 ⇒ 1, 20 ⇒ 2, 23 ⇒ 1, 111 ⇒ 1}
o7 : Tally
An entry of the form 18 ⇒ 3 in the result means that a Collatz sequence of length 18
was seen 3 times. If you want to use kewl functional programming, you could write instead
tally apply(1..30, length @@ collatz), with the exact same result.
Tutorial
2.3. Lines in the plane. What is the maximal number of regions R(n) that you can obtain
by drawing n lines in the plane R2?
We’ll explore this by hand first. Clearly
R(0) = 1, R(1) = 2, R(2) = 4.
Compute R(3).
EXPERIMENTAL MATHEMATICS 2022 7
Convince yourself that the following recurrence relation holds:
R(0) = 1
R(n) = R(n− 1) + n if n ≥ 1
Now implement it in either Macaulay2 or Mathematica!
Solution.
i8 : R = n -> if n == 0 then 1 else R(n-1)+n
o8 = R
o8 : FunctionClosure
If your function is called R, then you should be able to display the first few terms using
i9 : L=apply(20,R)
o9 = {1, 2, 4, 7, 11, 16, 22, 29, 37, 46, 56, 67, 79, 92, 106, 121, 137, 154, 172, 191}
o9 : List
or
L=Table[R[n], {n, 0, 19}]
and check that it matches the terms we already knew.
Very satisfying. How about a formula for R(n) though? We could think about it and
figure it out fairly quickly in this example, but let’s instead search for the first few terms
in the Online Encyclopedia of Integer Sequences (OEIS), either directly at oeis.org or using
Macaulay2’s built-in oeis.
Solution.
i10 : oeis L
o10 =

0. A000124 Central polygonal numbers (the Lazy Caterer’s sequence):
n(n+1)/2 + 1;
or, maximal number of pieces formed when slicing a pancake with n cuts.
1. A152947 a(n) = 1 + (n-2)*(n-1)/2.
2. A025739 Index of 9^n within sequence of numbers of form 9^i*10^j.

o10 : NumberedVerticalList
Obviously the first entry is the correct one.
Check that the closed form given by OEIS for R(n) matches the recurrence relation. Go one dimension up and play with it: what is the maximal number of regions
that can be obtained from n planes in R3?
As frivolous as the topic may seem (slicing a pancake with n cuts, indeed), it is an active
area of research. Look up hyperplane arrangements on the web, or [Sta12, Section 3.11].
2.4. The Fibonacci sequence. The famous Fibonacci sequence is given by the recurrence
relation
F0 = 1
F1 = 1
Fn = Fn−1 + Fn−2
8 P. ZINN-JUSTIN
We’ve learned in the previous exercise how to implement recurrence relations; define a func-
tion F that does that for the Fibanacci sequence.
Solution.
i11 : F = n -> if n==0 or n==1 then 1 else F(n-1)+F(n-2);
Now try to run it to compute say the first 35 terms, and let’s time the result:
i12 : time apply(35,F)
-- used 31.007 seconds
o12 = {1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1 597, 2 584, 4 181, 6 765, 10 946, 17 711, 28 657, 46 368, 75 025, 121 393, 196 418, 317 811, 514 229, 832 040, 1 346 269, 2 178 309, 3 524 578, 5 702 887, 9 227 465}
o12 : List
or Timing[Table[F[n],{n,0,34}]] in Mathematica.
Can you figure out why it is so slow? If you did, then you realize you need to modify your
program so it doesn’t keep recomputing the same entries of Fn! In Macaulay2, all you need
to do is to add memoize in front of the definition of F (or you can do it after the fact by
writing
i13 : F=memoize F;
or in Mathematica, by writing F[n_]:=F[n]=... where the dots stand for your previous
definition of F. Try that and time again your calculation of F (0), . . . , F (34); is it better
now?
i14 : time apply(35,F)
-- used 0.000046394 seconds
o14 = {1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1 597, 2 584, 4 181, 6 765, 10 946, 17 711, 28 657, 46 368, 75 025, 121 393, 196 418, 317 811, 514 229, 832 040, 1 346 269, 2 178 309, 3 524 578, 5 702 887, 9 227 465}
o14 : List
It is well-known that the following formula holds
(5) Fn =
φn+1+ − φn+1−
φ+ − φ− φ± =
1±√5
2
Check that the expression in the r.h.s. satisfies the recurrence above, thus proving the equal-
ity. Then implement this new formula for Fn by replacing φ± with their numerical values.
Check that it agrees with the old one on examples.
Solution.
i15 : phiplus = (1+sqrt 5)/2
o15 = 1.61803398874989
o15 : R (of precision 53)
i16 : phiminus = (1-sqrt 5)/2
o16 = −.618033988749895
o16 : R (of precision 53)
i17 : F’ = n -> (phiplus^(n+1)-phiminus^(n+1))/(phiplus-phiminus)
o17 = F’
o17 : FunctionClosure
i18 : apply(30, i-> round(F’(i+2)-F’(i+1)-F’ i))
o18 = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
o18 : List
EXPERIMENTAL MATHEMATICS 2022 9 Implement (5) exactly rather than numerically. Show that you have exact equal-
ity between the two different definitions of Fn.
Solution.
i19 : R=ZZ[phi]/(phi^2-phi-1);
i20 : F’ = n -> (phi^(n+1)-(-phi)^(-n-1))/(phi+phi^-1);
i21 : apply(30, i-> F’(i+2)-F’(i+1)-F’ i)
o21 = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
o21 : List
It is known that a positive integer m is part of the Fibonacci sequence if and only if at
least one of 5m2 + 4 or 5m2− 4 is a perfect square. Implement two functions that take m as
input: one that directly checks if m = Fn for some n, and one that uses the condition above.
Check that they agree.
Solution. For the first test we use the fact that Fn is monotonic to stop the search as soon
as the value of Fn exceeds m:
i22 : test1 = m -> (
n:=0;
while F n < m do n=n+1;
F n == m)
o22 = test1
o22 : FunctionClosure
To test if a number is a perfect square, we can use the exact same strategy:
i23 : isSquare = m -> (
n:=0;
while n^2 < m do n=n+1;
n^2 == m)
o23 = isSquare
o23 : FunctionClosure
or can one use the characterization that all exponents of the factorization of m must be even:
i24 : isSquare = m -> all(factor m, even @@ last)
o24 = isSquare
o24 : FunctionClosure
Then one has
i25 : test2 = m -> isSquare(5*m^2+4) or isSquare(5*m^2-4)
o25 = test2
o25 : FunctionClosure
and sure enough,
i26 : apply(1..30,test1) == apply(1..30,test2)
o26 = true
10 P. ZINN-JUSTIN
2.5. The number of Alternating Sign Matrices. Alternating Sign Matrices (ASMs) are
matrices with
(i) entries taken from {−1, 0, 1},
(ii) row- and column-sums are equal to 1,
(iii) in each row and each column, the non-zero entries alternate in sign.
For example, there are seven 3× 3 ASMs1 0 00 1 0
0 0 1
 ,
0 1 01 0 0
0 0 1
 ,
1 0 00 0 1
0 1 0
 ,
0 1 00 0 1
1 0 0
 ,
0 0 11 0 0
0 1 0
 ,
0 0 10 1 0
1 0 0
 ,
0 1 01 −1 1
0 1 0

We want to know how many ASMs of a given size there are. In order to enumerate ASMs
efficiently, we first transform them by taking partial column sums for each column yielding,
e.g., 
0 0 1 0 0
0 0 0 1 0
1 0 0 −1 1
0 1 −1 1 0
0 0 1 0 0
 →

0 0 1 0 0
0 0 1 1 0
1 0 1 0 1
1 1 0 1 1
1 1 1 1 1

Because of the alternating condition, the new matrix has only 0’s and 1’s, and because each
columns starts and ends with a 1, the bottom row has all ones. Next we record for each row
the location of the 1’s. In the example above we get
3
3 4
1 3 5
1 2 4 5
1 2 3 4 5
These triangles are called monotone triangles, and are defined by the properties that
entries in each row are strictly increasing, and every entry is weakly between the two entries
right below it.
Monotone triangles are easy to enumerate using a recurrence scheme. Let us consider
the slightly more general problem of finding the number F (a1, . . . , an) of monotone triangles
whose bottom row is a1, . . . , an. Given the bottom row, it follows from the conditions above
that the second-row from the bottom, let’s call it b1, . . . , bn−1, has to satisfy the following
restrictions:
a1 ≤ b1 ≤ a2 ≤ b2 ≤ . . . ≤ bn−1 ≤ an,(6)
b1 < b2 < . . . bn−1.(7)
The number F (a1, . . . , an) is then recursively given by
(8) F (a1, . . . , an) =
∑
(b1,...,bn−1)
F (b1, . . . , bn−1),
where the summation is over all b’s that satisfy (6) and (7).
EXPERIMENTAL MATHEMATICS 2022 11
After this lengthy introduction, let’s implement all this.
First, write a function, that, given an input a = (a1, . . . , an), generates all possible config-
urations of b = (b1, . . . , bn−1) satisfying (6) and (7). If this is too hard, you can look at the
suggested solutions in tutorial1.m2 or tutorial1.nb, and try to understand them:
i27 : -- this function returns, given (a_k), the set of sequences (b_k)
-- such that a_0 <= b_0 <= a_1 ... <= b_(n-2) <= a_(n-1)
T = a -> (
if #a == 0 then {}
else if #a == 1 then {{}}
else join apply(a_0..a_1, b0 -> apply( -- first fix b_0
if b0 < a_1 then T drop(a,1) else T(prepend(b0+1,drop(a,2))),
-- find other b’s being careful that if b_0=a_1, a_1 is excluded
bs -> prepend(b0,bs) -- insert b_0 in front of the other b’s
)
)
)
o27 = T
o27 : FunctionClosure
Next, implement (8) and compute F (1, 2, . . . , n), the number of n × n ASMs. Compare
your answer with the sequence given in the lectures.
Solution.
i28 : F = memoize ( a -> if length a == 0 then 1 else sum(T a,F) )
o28 = F
o28 : FunctionClosure
i29 : A = n -> F (1..n)
o29 = A
o29 : FunctionClosure
i30 : apply(8,A)
o30 = {1, 1, 2, 7, 42, 429, 7 436, 218 348}
o30 : List If you’re more ambitious, you can write a function that generates (rather than
enumerates) all the ASMs (or monotone triangles) of a given size!
Horizontally symmetric ASMs. Horizontally symmetric ASMs (HSASMs) of size (2n +
1) × (2n + 1) are ASMs which are symmetric with respect to reflection in the horizontal
symmetry axis. There is, for example, only one HSASM of size 3. Show that HSASMs of
size (2n + 1) × (2n + 1) are enumerated by F (1, 3, 5, . . . , 2n + 1). Compute the number of
HSASMs and derive a closed form formula following the same steps as in the lecture.
Solution.
i31 : H = n -> F apply(1..n,i->2*i)
o31 = H
o31 : FunctionClosure
i32 : apply(8,H)
12 P. ZINN-JUSTIN
o32 = {1, 1, 3, 26, 646, 45 885, 9 304 650, 5 382 618 660}
o32 : List
i33 : oeis oo
o33 =

0. A005156 Number of alternating sign 2n+1 X 2n+1 matrices
symmetric about the vertical axis (VSASM’s); also
2n X 2n off-diagonally symmetric alternating sign matrices (OSASM’s).

o33 : NumberedVerticalList
2.6. Pascal mod 2. The Pascal triangle is the triangle of binomial coefficients. We’re
not going to bother to define binomials inductively, but rather use the built-in functions
binomial or Binomial. We’re interested in these coefficients mod 2; something like
i34 : binomial(6,4)%2
o34 = 1
or
Mod[Binomial[6,4],2]
Now create a nested list L of binomials mod 2 (let’s say, up to n = 31); the output should
look like {{1}, {1, 1}, {1, 0, 1}, {1, 1, 1, 1}, . . .}.
Solution.
i35 : L=apply(32,i->apply(i+1,j->binomial(i,j)%2))
o35 = {{1} , {1, 1} , {1, 0, 1} , {1, 1, 1, 1} , {1, 0, 0, 0, 1} , {1, 1, 0, 0, 1, 1} , {1, 0, 1, 0, 1, 0, 1} , {1, 1, 1, 1, 1, 1, 1, 1} , {1, 0, 0, 0, 0, 0, 0, 0, 1} , {1, 1, 0, 0, 0, 0, 0, 0, 1, 1} , {1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1} , {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1} , {1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1} , {1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1} , {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1} , {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} , {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1} , {1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1} , {1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1} , {1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1} , {1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1} , {1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1} , {1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1} , {1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1} , {1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1} , {1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1} , {1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1} , {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1} , {1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1} , {1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1} , {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1} , {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}}
o35 : List
Finally, draw the triangle as a picture where black/white squares represent 0s and 1s. In
Macaulay2, this can be achieved with
i36 : needsPackage "VectorGraphics"; matrixPlot L
EXPERIMENTAL MATHEMATICS 2022 13
o37 =
o37 : GraphicsList
Mathematica has fancier drawing abilities, and there are quite a few options to produce the
picture, see e.g.
https://www.wolfram.com/language/fast-introduction-for-math-students/en/plots-in-2d/
and
https://www.wolfram.com/language/fast-introduction-for-math-students/en/more-plots-in-2d/.
You can for example use Image or ArrayPlot.
In any case, does the picture look familiar? Try again with other values of the modulo,
e.g., mod 3.
14 P. ZINN-JUSTIN
3. Gauss’s agM
From Gauss’s diary, entry dated 30 May 1799:
Terminum medium arithmetico–geometricum inter 1 et
√
2 esse = pi
$
usque
ad figuram undecimam comprobavimus, qua re demonstrata prorsus novus
campus in analysis certo aperietur.
For those as Latin-challenged as I am, this translates to
We have established that the arithmetic–geometric mean between 1 and
√
2 is
= pi
$
to the eleventh decimal place; the demonstration of this fact will surely
open an entirely new field of analysis.
This probably raises more questions than it answers. So let’s look at it in more detail.
3.1. Two sequences, one limit. Consider the two sequences of real numbers (an) and (bn)
defined by the recursions
a0 =
√
2, an+1 =
an + bn
2
b0 = 1, bn+1 =
√
anbn
This looks a little strange, as the two definitions are intertwined (and what’s up with the
initial values
√
2 and 1?), but we can recognise some familiar features, at least locally. The
recursive rule defining an+1 takes the arithmetic mean of the numbers an and bn. And the
recursive rule defining bn+1 takes the geometric mean of the same two numbers.
I wonder how the two sequences behave as n increases. If I were Gauss, I would compute
the first few terms by hand (to 11 decimal places, indeed). Having the luxury of living in
the era of funny cat YouTube videos (and incidental infrastructure), I’ll instead call upon
Macaulay2:
i1 : a = memoize ( n -> (a(n-1)+b(n-1))/2, { 0 => sqrt 2 } );
i2 : b = memoize ( n -> sqrt(a(n-1)*b(n-1)), { 0 => 1 } );
Having defined the two sequences, I can now ask for values:
i3 : apply(5,a)
o3 = {1.41421, 1.20711, 1.19816, 1.19814, 1.19814}
o3 : List
i4 : a 5
o4 = 1.19814023473559
o4 : R (of precision 53)
One might guess now that (an) converges to a number close to 1.19814023473559. How
about (bn)?
i5 : apply(6,b)
o5 = {1, 1.18921, 1.19812, 1.19814, 1.19814, 1.19814}
o5 : List
This also seems to converge, and to the same limit? Let’s be bold and proclaim:
Proposition 1. Given any starting values a0 = α ≥ b0 = β ∈ R>0, the sequences (an) and
(bn) both converge to the same limit M(α, β).
EXPERIMENTAL MATHEMATICS 2022 15
Proof. More of a sketch, really.
There are two things to prove: that the sequences converge, and that they have the same
limit1.
Here is one approach to this, courtesy of Wikipedia:
• bn ≤ an for all n. This is clear for n = 0 and true by simple algebraic manipulation
in general (the arithmetic mean of two numbers is never smaller than their geometric
mean).
• Use the previous part to note that bn+1 ≥ bn, so the sequence (bn) is non-decreasing.
• Note that bn ≤ α for all n ≥ 0, so the sequence (bn) is bounded above. Hence it has
a limit L.
• Note that bn ≥ β > 0 for all n so L ≥ β > 0, in particular L 6= 0.
• Note that
an =
b2n+1
bn
,
so using a theorem about the limit of the quotient of two convergent sequences we
conclude that (an) converges to
L2
L
= L.

The real number M(α, β) is called the arithmetic-geometric mean (agM) of α and β. The
above numerical experiment leads us to believe that
M(
√
2, 1) = 1.19814023473559 . . .
But Gauss’s diary entry said more about this value. What’s that about?
3.2. An elliptic integral. Consider the lemniscate (of Bernoulli), a plane curve given in
polar coordinates (r, θ) by the equation
r2 = cos(2θ).
It would be nice to graph it, wouldn’t it? As far as I know, neither Macaulay2 nor
Mathematica have a built-in command for implicit polar plots. After a bit of algebraic
manipulation with x = r cos θ, y = r sin θ, I get the implicit Cartesian equation
(x2 + y2)2 = x2 − y2.
This is more amenable to plotting:
i6 : needsPackage "VectorGraphics";
i7 : R=RR[x,y];
i8 : P=(x^2+y^2)^2-(x^2-y^2)
o8 = x4 + 2x2y2 + y4 − x2 + y2
o8 : R
i9 : plot2d(P,{[-1.2,-1.2],[1.2,1.2]})
1It would not be sufficient to only show that the sequence of differences (an − bn) converges to 0
16 P. ZINN-JUSTIN
o9 =
-1
-1
0
0
1
1
-1-1
00
11
o9 : GraphicsList
which produces the pretty infinity sign in the picture.
Or you can follow a Mathematica lead from StackExchange:
https://mathematica.stackexchange.com/a/549
ContourPlot[
Evaluate@With[
{r = Sqrt[x^2 + y^2],
theta = ArcTan[x, y]},
r^2 - Cos[2*theta] == 0
],
{x, -1, 1}, {y, -1, 1}
]
(Digression: there are plenty of pretty pictures (and Mathematica code) related to the
lemniscate in the paper [LS10].)
The question is: what is the arclength of this curve?
The standard arclength formula in polar coordinates gives
4
∫ pi/4
0
(
r2 +
(
dr
dθ
)2)1/2
dθ = 4
∫ pi/4
0
(cos(2θ))−1/2 dθ.
If we introduce a new variable α with the property that cos(2θ) = cos2 ϕ, the integral
becomes
4
∫ pi/2
0
(
1 + cos2 ϕ
)−1/2
dϕ = 4
∫ pi/2
0
(
2 cos2 ϕ+ sin2 ϕ
)−1/2
dϕ.
One-half of this value is what Gauss denoted by $ in his diary entry.
EXPERIMENTAL MATHEMATICS 2022 17
We shouldn’t just believe him like this. Let’s check his calculations by estimating the
integral numerically. In Macaulay2:
i10 : pi/integrate(x->2/sqrt(2*(cos x)^2+(sin x)^2),0,pi/2)
o10 = 1.19814023473538
o10 : R (of precision 53)
In Mathematica:
Pi/NIntegrate[2/Sqrt[2*Cos[t]^2+Sin[t]^2],{t,0,Pi/2},
{WorkingPrecision->15}]
Fine, Gauss seems to have gotten his decimals right. More generally, he proved
Theorem 1. If α ≥ β > 0 then∫ pi/2
0
(
α2 cos2 ϕ+ β2 sin2 ϕ
)−1/2
dϕ =
pi
2M(α, β)
(These are examples of so-called elliptic integrals, and the arithmetic-geometric mean
is a very efficient way of approximating them. An elliptic integral is an integral of the
form
∫
dt F (t,
√
P (t)) where F is a rational function, and P is a polynomial of degree 3 or
4. What’s elliptic about it? Finding the arc length of an ellipse falls in this category of
integrals. The curve y2 = P (t) is known as an elliptic curve. Elliptic curves play a key role
in a number of areas of mathematics, in particular number theory, and have applications in
cryptography!)
Proof. Write
I(α, β) =
∫ pi/2
0
(
α2 cos2 ϕ+ β2 sin2 ϕ
)−1/2
dϕ
We introduce a variable ϕ′ with the property that
(9) sinϕ =
2α sinϕ′
α + β + (α− β) sin2 ϕ′
Then some secret magic sauce (in Gauss’s words “after the development has been done
correctly, it will be seen”2):(
α2 cos2 ϕ+ β2 sin2 ϕ
)−1/2
dϕ =
(
a21 cos
2 ϕ′ + b21 sin
2 ϕ′
)−1/2
dϕ′
which leads us to the conclusion
I(a0 = α, b0 = β) = I(a1, b1).
A moment’s further thought brings us to continue this as
I(α, β) = I(a1, b1) = I(a2, b2) = · · · = I(an, bn) = . . .
2Jacobi must have been as unimpressed as us by Gauss weaseling out of this, because he wrote down a
couple of intermediate steps: Given the relation (9) between ϕ and ϕ′, show that
cosϕ =
(2 cosϕ′)
(
a21 cos
2 ϕ′ + b21 sin
2 ϕ′
)1/2
α+ β + (α− β) sin2 ϕ′
and then that (
α2 cos2 ϕ+ β2 sin2 ϕ
)1/2
= α
α+ β − (α− β) sin2 ϕ′
α+ β + (α− β) sin2 ϕ′
After this and implicitly differentiating (9), it’s all smooth sailing to Gauss’s equality of differentials.
18 P. ZINN-JUSTIN
and some judicious (real) analysis allows us to pass to the limit and get
I(α, β) = I(M(α, β),M(α, β)) =
pi
2M(α, β)
as the integral I(c, c) is a piece of cake.
3.3. Built-in agM commands. You might wonder if the arithmetic-geometric mean is
already implemented in the software we’re playing with.
But of course! In Macaulay2, it is agm:
i11 : agm(1,sqrt 2)
o11 = 1.19814023473559
o11 : R (of precision 53)
In Mathematica, it is ArithmeticGeometricMean.
3.4. Relation to hypergeometric series. Because of the homogeneity relationM(tα, tβ) =
tM(α, β), we can assume without loss of generality α = 1, 0 < β < 1. Setting β =
√
1− z
and expanding in power series in z under the integral results in
(10)
1
M(1,
√
1− z) =
2
pi
I(1,
√
1− z) =
∞∑
n=0
(
2n
n
)2 ( z
16
)n
Such a series is a special case of the hypergeometric series 2F1(a, b; c; z) defined by
(11) 2F1(a, b; c; z) =
∞∑
n=0
(a)n(b)n
(c)nn!
zn,
where (a)0 = 1 and for n a positive integer
(a)n = a(a+ 1) · · · (a+ n− 1),
It is not hard to check that the series in (10) corresponds to a = b = 1/2, c = 1, so that
1
M(1,
√
1− z) = 2F1(1/2, 1/2; 1; z)
We shall encounter hypergeometric series again later in the course.
3.5. Recognizing the constant. Can M(1,
√
2) be expressed in terms of known functions?
Yes! Recall the Γ function, defined by Γ(z) =
∫∞
0
xz−1e−xdx for Re z > 0. Then
M(1,
√
2) =
Γ(3/4)
Γ(5/4)
√
pi
2
Want to know more? There is plenty more where this came from, namely David Cox’s
articles [Cox85] (the shorter version) and [Cox84] (the longer one).
Ah, I can’t resist mentioning one more thing, which I learned from Cox’s papers. In 1973,
Salamin discovered the formula
pi =
2M(
√
2, 1)2
1−∑∞n=1 2n+1(a2n − b2n)
which, as you can see, uses an infinite sum involving the terms an and bn of the sequence we
started with, as well as the common limit M(
√
2, 1) of these two sequences. The existence
of such a formula for pi is weird enough, but it turns out [Sal76] that it’s actually a pretty
efficient way to compute lots of digits of pi (in that the number of significant digits doubles
at each step).
EXPERIMENTAL MATHEMATICS 2022 19
Tutorial
3.6. Numerical precision. Both Macaulay2 and Mathematica support arbitrary precision
numerical calculations. In Macaulay2, the desired numerical precision can be specified with
numeric:
i12 : numeric_100 pi
o12 = 3.14159265358979323846264338328
o12 : R (of precision 100)
For numbers, one can use the following shorthand:
i13 : 1/3p100
o13 = .333333333333333333333333333333
o13 : R (of precision 100)
In each case the precision is measured in binary digits (a precision of 100 is roughly 30 digits
in base 10).
Or one can simply (as one of you suggested!) set once and for all the default precision:
i14 : defaultPrecision=100;
and then all subsequent calculations will be done using this precision.
For Mathematica, there is some information about numeric precision and interval arith-
metic at
https://reference.wolfram.com/language/tutorial/Numbers.html
First check the last identity of the lecture up to 30 digits: that
M(1,
√
2) =
Γ(3/4)
Γ(5/4)
√
pi
2
In both Macaulay2 and Mathematica, the Γ function is Gamma.
Solution.
i15 : agm(1,sqrt 2)
o15 = 1.19814023473559220743992249228
o15 : R (of precision 100)
i16 : sqrt(pi/4)*Gamma(3/4)/Gamma(5/4)
o16 = 1.19814023473559220743992249228
o16 : R (of precision 100)
Now consider the real number
α =
pi10
ζ(10)
where ζ(s) =
∑∞
n=1
1
ns
is implemented as zeta in Macaulay2 (resp. Zeta in Mathematica).
Compute α to default precision and guess what the exact value of α might be.
Increase the precision a few times to see if your guess persists.
Solution.
20 P. ZINN-JUSTIN
i17 : pi^10/zeta(10)
o17 = 93554.9999999999999999999999998
o17 : R (of precision 100)
i18 : (numeric_200 pi)^10/zeta(10p200)
o18 = 93555
o18 : R (of precision 200)
Consider the real number
β =
(
epi
√
163 − 744
)1/3
Compute β to default precision and venture a guess about the value of β.
Increase the working precision and recompute β. Do you need to adjust your guess? If
not, maybe increase the precision some more.
Solution.
i19 : (exp(pi*sqrt 163)-744)^(1/3)
o19 = 640320.000000000000000000000002
o19 : R (of precision 100)
i20 : defaultPrecision=200;
i21 : (exp(pi*sqrt 163)-744)^(1/3)
o21 = 640319.999999999999999999999999390317352319470126502835539026
o21 : R (of precision 200)
Think about what strategies you might employ to check whether what you are seeing is
just a numerical glitch.
A good approach to this type of question relies on interval arithmetic.
3.7. Interval arithmetic. This is a continuation of the previous exercise.
In Macaulay2, one can write e.g.,
i22 : defaultPrecision=20;
i23 : mypi = interval(3.14,3.15)
o23 =
[
3.14, 3.15
]
o23 : R (of precision 20)
i24 : (exp(sqrt 163 * mypi)-744)^(1/3)
o24 =
[
635979, 663660
]
o24 : R (of precision 20)
If the interval does not contain the expected value, this means there is inequality (but if it
does, we can’t conclude anything yet!).
Because approximating constants is so common, there is in fact a built-in function:
i25 : mypi=numericInterval(100,pi)
o25 =
[
3.14159265358979323846264338328, 3.14159265358979323846264338328
]
o25 : R (of precision 100)
Note that calculations are still performed with some given precision; if you want the
interval to be smaller and smaller, you need to make sure every constant has the prescribed
precision!
EXPERIMENTAL MATHEMATICS 2022 21
Experiment until you have a definite conclusion about the value of β defined in the previous
exercise.
i26 : defaultPrecision=120;
i27 : mypi=numericInterval(100,pi);
i28 : (exp(sqrt 163 * mypi)-744)^(1/3)
o28 =
[
640319.999999999999999999999998, 640320.000000000000000000000009
]
o28 : R (of precision 100)
i29 : mypi=numericInterval(120,pi);
i30 : (exp(sqrt 163 * mypi)-744)^(1/3)
o30 =
[
640319.999999999999999999999999390303, 640319.999999999999999999999999390324
]
o30 : R (of precision 120)
Given some accuracy bound ε, implement your own function myagm that returns the agM
of two positive real numbers a and b in interval arithmetic, where the width of the interval
is less than ε. You can reuse the functions a and b defined in the lectures.
Solution. Ignoring the suggestion, let’s reprogram directly:
i31 : defaultPrecision=53;
i32 : eps=1e-10;
i33 : myagm = (a,b) -> (
while abs(a-b)>eps do (a,b)=((a+b)/2,sqrt(a*b));
interval(a,b)
)
o33 = myagm
o33 : FunctionClosure
3.8. Continued fractions. In the first exercise, we’ve found that pi
10
ζ(10)
is an integer. This
leads to the obvious question: given an integer k ≥ 2, is pik
ζ(k)
an integer? If not, is it at least
a rational number?
Experiment numerically, and formulate a conjecture.
Solution. It seems likely that pi
k
ζ(k)
is a rational number (but not always an integer) when k
is even; nothing obvious appears for k odd.
We shall use continued fractions to find the best rational approximation of a real number.
A continued fraction expansion of a positive real number β is an expression of the form
β = a0 +
1
a1 +
1
a2 +
1
a3 +
1
. . .
where a0 ∈ Z≥0, a1, a2, · · · ∈ Z>0.
22 P. ZINN-JUSTIN
Given β, the continued fraction expansion is easily computed by the recursion
β0 = β
an = bβnc for n ≥ 0
βn =
1
βn−1 − an−1 for n ≥ 1
The truncation at an of a continued fraction is called its n-th convergent
a0 +
1
a1 +
1
a2 +
1
· · ·+ 1
an
=
pn
qn
If you harbour any nostalgia for the good ol’ Real Analysis days, try your hand at this:
prove that the sequence (pn/qn) converges as n→∞. One possible approach is to first show
that
pn−1qn − pnqn−1 = (−1)n
then conclude that the sequence (pn/qn) is Cauchy.
Write a function that, given n and a real number β, returns pn/qn (and optionally prints
the sequence a0, . . . , an). b. . .c is implemented as floor, resp. Floor. Try it on the golden
ratio 1+
√
5
2
, the square roots of small integers, e, pi. . . What do you observe? Can you prove
a statement about numbers with periodic continued fraction expansion?
Solution. Here’s one possible implementation:
i34 : defaultPrecision=100;
i35 : contfrac = (x,n) -> (
L:=for i from 1 to n
list a=floor x
do if xprint L;
s:=0; scan(reverse L,a->s=1/(s+a));
1/s
)
o35 = contfrac
o35 : FunctionClosure
Experimenting, we see that solutions of quadratic equations with integer coefficients, such
as
i36 : contfrac((1+sqrt 5)/2,20)
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
o36 = 10 946
6 765
o36 : Q
i37 : contfrac(sqrt 42,20)
{6, 2, 12, 2, 12, 2, 12, 2, 12, 2, 12, 2, 12, 2, 12, 2, 12, 2, 12, 2}
o37 = 69 544 807 113 937
10 730 996 710 148
o37 : Q
EXPERIMENTAL MATHEMATICS 2022 23
have (eventually) periodic continued fraction; this is not hard to prove directly. One can
discern a pattern in say
i38 : contfrac(exp 1,20)
{2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1}
o38 = 28 245 729
10 391 023
o38 : Q
but not in
i39 : contfrac(1*pi,20) -- 1* due to current bug in Macaulay2, soon to be fixed
{3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2}
o39 = 14 885 392 687
4 738 167 652
o39 : Q
A number β is rational if and only if its continued fraction expansion is finite. Experiment
on pi
k
ζ(k)
, minding accuracy issues.
Solution.
i40 : L=apply(10, k -> contfrac(pi^(2*k)/zeta(2*k),100))
{−2}
{6}
{90}
{945}
{9 450}
{93 554, 1}
{924 041, 1, 3, 1, 2, 2, 1, 14}
{9 121 612, 2}
{90 030 844, 1, 30, 2, 4, 1, 2, 1, 1, 1}
{888 579 011, 9, 1, 1, 3, 1, 1, 2, 1, 1, 4, 10}
o40 =
{−2, 6, 90, 945, 9 450, 93 555, 638 512 875
691
, 18 243 225
2
, 325 641 566 250
3 617
, 38 979 295 480 125
43 867
}
o40 : List
Hopefully, at this stage, you’ve found a certain sequence of numerators/denominators.
Ask the OEIS for help recognizing them!
Solution.
i41 : oeis apply(L,numerator)
o41 = {0. A002432 Denominators of zeta(2*n)/Pi^(2*n). }
o41 : NumberedVerticalList
i42 : oeis apply(L,denominator)
o42 = {0. A046988 Numerators of zeta(2*n)/Pi^(2*n). }
o42 : NumberedVerticalList
3.9. agM and Ramanujan identity. Consider the function
A(q) =
(∑
n∈Z
qn
2
)2
for |q| < 1.
Implement A numerically to some given accuracy.
24 P. ZINN-JUSTIN
Solution. One possible implementation is
i43 : defaultPrecision=53;
i44 : eps=1e-10;
i45 : A = q -> (
s:=1;
i:=1;
while (x:=q^(i^2))>eps do (s=s+2*x; i=i+1);
s^2
)
o45 = A
o45 : FunctionClosure
Now compute A(e−pi). Do you recognize this number? (Hint: it is closely related to
M(1,
√
2))
Solution.
i46 : A(exp(-pi))
o46 = 1.18034059901381
o46 : R (of precision 53)
i47 : A(exp(-pi))*agm(1,sqrt 2)/sqrt 2
o47 = .999999999998065
o47 : R (of precision 53)
This is a special case of an identity of Ramanujan, namely that if q is defined by
q = exp
(
−pi M(1,
√
x)
M(1,
√
1− x)
)
for 0 < x < 1, then A(q) = 1/M(1,
√
x). (can you see why it’s a special case?)
EXPERIMENTAL MATHEMATICS 2022 25
4. Constant recognition
In the previous chapter, and in particular the tutorial, we have started investigating a
common task in computer-assisted mathematics: obtaining a numerical approximation to
a number we are interested in. and then “recognising” this; i.e., whether this number is
of a special type (e.g. rational, or algebraic, or a simple combination of other numbers we
like such as pi, e, M(
√
2, 1))? There are surprisingly robust ways of approaching such an
ill-defined question, as we’ll see now.
4.1. Finding the integer relation. In the tutorial, we introduced the Riemann zeta-
function, which is a function of a complex variable s defined, for <(s) > 1, by the infinite
series
ζ(s) =
∞∑
n=1
1
ns
Despite having Riemann’s (1826–1866) name attached to it, particular aspects had been
considered two hundred years before Riemann. As one such example, Mengoli asked in 1650
for the value
ζ(2) =
1
12
+
1
22
+
1
32
= . . .
What Mengoli meant by “value” was not an approximation such as
i1 : zeta 2
o1 = 1.64493406684823
o1 : R (of precision 53)
which he could compute himself (maybe to slightly fewer decimals). It was the exact value,
in closed form, a notion that is more metaphysical than mathematical. He was basically
saying “I want a simple and pretty answer involving terms that I know already, and a proof
that the answer is correct.” This became known as the Basel problem and Euler solved it
in 1741.
In the tutorial, we “guessed” that ζ(2) = pi
2
6
:
i2 : zeta 2 * 6 / pi^2
o2 = 1
o2 : R (of precision 53)
Let’s take some of the previous discussion and twist it around a little bit. Suppose I consider
the real numbers ζ(2) and pi2 and I ask: are there integers a and b such that
(12) a ζ(2) + b pi2 = 0?
With the 20/20 hindsight afforded by the previous section, we can confidently answer:
Yes, take a = 6 and b = −1. But the point is that in Equation (12) we organised this in
the form of a linear relation with integer coefficients between the two numbers ζ(2) and pi2.
And one can consider linear relations involving more than two numbers.
More precisely, given a vector x = (x1, x2, . . . , xn) ∈ Rn, we define an integer relation for
x to be a nonzero vector m = (m1,m2, . . . ,mn) ∈ Zn with integer entries such that
m · x = m1x1 +m2x2 + · · ·+mnxn = 0.
Of course, there is no guarantee that such a magical m exists. (Take, for instance, x =
(1,
√
2).) This leads us to state the
26 P. ZINN-JUSTIN
Integer relation problem: Given x ∈ Rn, either find a “small” integer relation m for
x or prove that no such “small” integer relation exists.
You may find the rather imprecise use of the adjective “small” distasteful. In that case,
there is a more assertive version of the problem that fixes a bound 2k and asks for m ∈ Zn
with ‖m‖ ≤ 2n+k or a proof that there is no m ∈ Zn with ‖m‖ < 2k.
There is yet another version that is most attractive in practice, where we take the input
vector x ∈ Zn as well. The following example shows how this might come about.
Example 4. Let
x1 = arctan(1) = 0.785398 . . .
x2 = arctan(1/5) = 0.197395 . . .
x3 = arctan(1/239) = 0.004184 . . .
Can we find an integer relation m for x = (x1, x2, x3):
m1x1 +m2x2 +m3x3 = 0?
We turn this into a question about integers by fixing a multiplier A, say A = 106, and
considering
m1bAx1c+m2bAx2c+m3bAx3c ≈ 0,
that is
785398m1 + 197395m2 + 4184m3 ≈ 0.
There are several approaches to finding m1,m2,m3 and we will be looking in more detail at
one of them, the LLL algorithm. For now I will just say that this algorithm takes the matrix
i3 : M = matrix { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 }, { 785398, 197395, 4184 } }
o3 =

1 0 0
0 1 0
0 0 1
785 398 197 395 4 184

o3 : Matrix Z4 ←− Z3
and returns the matrix
i4 : LLL M
o4 =

1 −13 −52
−4 58 203
1 −296 184
2 272 345

o4 : Matrix Z4 ←− Z3
(LLL is also implemented in Mathematica as LatticeReduce, except Mathematica never
follows normal mathematical conventions, so all matrices need to be transposed.) The first
column of this matrix will be telling us that
785398 · 1 + 197395 · (−4) + 4184 · 1 = 2,
or, going back to our original real numbers, that
0.785398 · 1 + 0.197395 · (−4) + 0.004184 · 1 = 2
106
≈ 0.
In other words, m = (1,−4, 1) seems to be an integer relation.
EXPERIMENTAL MATHEMATICS 2022 27
Indeed, Machin found this formula in 1706:
arctan(1) = 4 arctan(1/5)− arctan(1/239).
It works remarkably well as a rapidly convergent approximation to pi = 4 arctan(1).
Obvious questions remain. What is this LLL algorithm? What exactly is it doing to that
matrix? And how come the first column of the resulting matrix gives us the integer relation
we have been looking for?
A quick answer is that LLL is a lattice reduction algorithm, and that the integer rela-
tion problem can be solved via lattice reduction. To make sense of this, we need to know
something about lattices.
4.2. Lattice reduction. Fix a natural number n and let V be an n-dimensional real vector
space endowed with an inner product u ·v. (You may think of V as being Rn with the usual
dot product.)
A lattice in V is a subset L ⊂ V such that there exists a basis {b1,b2, . . . ,bn} of V with
L = SpanZ{b1,b2, . . . ,bn} = {r1b1 + r2b2 + · · ·+ rnbn | r1, r2, . . . , rn ∈ Z}.
We say that {b1, . . . ,bn} is a Z-basis for the lattice L, and we call n the rank of L. The
determinant d(L) of L is the determinant of the matrix with columns b1, . . . ,bn. This turns
out to be independent of the choice of basis.
It would be good to recall the setup of Gram–Schmidt orthogonalisation. Let b1, . . . ,bn
be a basis for V . For i = 1, . . . , n let Vi = SpanR{b1, . . . ,bi}.
The Gram–Schmidt process returns vectors b∗1, . . . ,b
∗
n ∈ V and scalars µij ∈ R for 1 ≤
j < i ≤ n, defined inductively by
b∗i = projV ⊥i−1(bi) = bi − projVi−1(bi) = bi −
i−1∑
j=1
µijb
∗
j
µij =
bi · b∗j
b∗j · b∗j
with b∗1 = b1.
Note that, for all i = 2, . . . , n,
Vi−1 = SpanR{b1, . . . ,bi−1} = SpanR{b∗1, . . . ,b∗i−1}
and b∗1, . . . ,b
∗
n is an orthogonal basis of V .
Let’s go back to the setup of a lattice L ⊂ V . We say that a basis c1, . . . , cn for L is
(LLL-)reduced if
|µij| ≤ 1
2
for all 1 ≤ j < i ≤ n
‖c∗i + µi,i−1c∗i−1‖2 ≥
3
4
‖c∗i−1‖2 for all 1 < i ≤ n.
Here 3/4 can be replaced by anything in the open interval (1/4, 1).
The advantage of a reduced basis is that its vectors are in some sense small:
28 P. ZINN-JUSTIN
Proposition 2. If {c1, . . . , cn} is a reduced basis of a lattice L, then
‖cj‖ ≤ 2(i−1)/2‖c∗i ‖ for all 1 ≤ j ≤ i ≤ n
d(L) ≤
∏
i
‖ci‖ ≤ 2n(n−1)/4d(L)
‖c1‖ ≤ 2(n−1)/4d(L)1/n
‖c1‖ ≤ 2(n−1)/2‖x‖ for all x ∈ L− {0}
It is worth dwelling a little on this last inequality. If the constant (once L is fixed)
multiplier 2n−1 were not there, we would have that c1 is a shortest nonzero vector in L. This
may seem like a desirable outcome (and indeed many problems rely on finding a shortest
vector), but it has the great disadvantage that its time complexity is exponential in the
dimension n. A reduced basis gives up some control on the size of c1, to the extent shown in
the last inequality of the Proposition. What is gained however is that this can be computed
in time polynomial in the dimension n.
The latter is achieved by the LLL reduction algorithm. Its name is derived from its
authors, Arjen Lenstra, Hendrik Lenstra, and La´szlo´ Lova´sz. The algorithm starts with some
given basis of L and iteratively modifies it to achieve a reduced one. It is not particularly
complicated, and the exposition in the original paper [LLL82] is quite good.
Example 5. Going back to the situation of Example 4, the lattice is Z-spanned by the columns
of the matrix
M =

1 0 0
0 1 0
0 0 1
785398 197395 4184

The norms of the three basis vectors are roughly
785398.0, 197395.0, 4184.0
The LLL algorithm returns the matrix
1 −13 −52
−4 58 203
1 −296 184
2 272 345

whose columns are the vectors in the reduced basis, with norms roughly
4.7, 406.4, 443.6
Remember that Proposition 2 only guarantees that the first basis vector is within a factor
of
√
23−1 = 2 of the shortest nonzero vector. In this example however, it can be checked
that the first basis vector is actually a shortest vector. This does happen sometimes, an
instance of the fact that the worst-case performance and the average-case performance of an
algorithm can be quite different.
How does lattice reduction help with finding integer relations? (It may be good to follow
along with Example 4.) Given x1, . . . , xn ∈ R, let x′1, . . . , x′n ∈ R be close approximations
(such as, for instance, truncating after a certain number of digits). Let A be a multiplier.
EXPERIMENTAL MATHEMATICS 2022 29
Consider the Z-linear map Zn → Rn+1 given bym1...
mn
 7→

m1
...
mn
A
∑
imix
′
i

Let L denote the image of this map; it is a lattice (of rank n) inside an n-dimensional
subspace V of Rn+1.
In other words, we consider the matrix representation of the Z-linear map above
M =

1 0 . . . 0
0 1 . . . 0
0 0 . . . 0
...
...
. . .
...
0 0 . . . 1
Ax′1 Ax
′
2 . . . Ax
′
n

Letting b1,b2, . . . ,bn denote the columns of this matrix, we have
L = SpanZ{b1, . . . ,bn} ⊂ V = SpanR{b1, . . . ,bn}
The crucial observation is that
m1
...
mn
ε
 ∈ L if and only if m1x′1 + · · ·+mnx′n = εA.
If ε is relatively small then (since x′i ≈ xi)
m1x1 + · · ·+mnxn ≈ 0.
So we are interested in finding elements of L that have small components. But, as we have
seen above, the first vector in a reduced basis for L will be relatively small. So we apply
LLL to the basis b1, . . . ,bn of L, obtain a reduced basis c1, . . . , cn, and take the relation
determined by the first vector c1.
Tutorial
4.3. Lattices. Write a function lattice2d which takes as input (v1,v2,range1,range2)
(where v1 and v2 are vectors, typically 2d vectors) and returns a list of the points in the
lattice with basis {v1, v2} whose first coefficient (i.e., the coefficient of v1) is in the range
range1 and the second is in the range range2.
Note that Macaulay2 distinguishes between vectors and lists: one can go back and forth
between the two using vector and entries.
Repeat the previous exercise, this time with triplets of vectors (in three, or possibly higher,
dimensional space).
Solution. (both at once)
i5 : lattice2d = (v1,v2,r1,r2) -> toList splice apply(r1,
i -> apply (r2, j -> i*v1+j*v2))
o5 = lattice2d
o5 : FunctionClosure
30 P. ZINN-JUSTIN
i6 :
lattice3d = (v1,v2,v3,r1,r2,r3) -> toList deepSplice apply(r1,
i -> apply (r2, j -> apply(r3, k ->
i*v1+j*v2+k*v3)))
o6 = lattice3d
o6 : FunctionClosure
i7 :
Experiment with your function and the square lattice (basis (1, 0) and (0, 1)):
L = lattice2d(vector{1,0},vector{0,1},-2..2,-2..2)
o7 =
{( −2
−2
)
,
( −2
−1
)
,
( −2
0
)
,
( −2
1
)
,
( −2
2
)
,
( −1
−2
)
,
( −1
−1
)
,
( −1
0
)
,
( −1
1
)
,
( −1
2
)
,
(
0
−2
)
,
(
0
−1
)
,
(
0
0
)
,
(
0
1
)
,
(
0
2
)
,
(
1
−2
)
,
(
1
−1
)
,
(
1
0
)
,
(
1
1
)
,
(
1
2
)
,
(
2
−2
)
,
(
2
−1
)
,
(
2
0
)
,
(
2
1
)
,
(
2
2
)}
o7 : List
In Macaulay2 you can use
i8 : needsPackage "VectorGraphics";
i9 : listPlot(L,Join=>false)
o9 = x
y
o9 : GraphicsList
to visualize the list L of points.
Try a couple more lattices.
Compare the plots for lattice2d for the lattices given by
• (1, 0) and (0, 1), both ranges -2..2
• (2, 3) and (3, 5), both ranges -20..20
Solution.
EXPERIMENTAL MATHEMATICS 2022 31
i10 : L = lattice2d(vector{2,3},vector{3,5},-20..20,-20..20);
You may want to restrict the viewing window for the second plot; in Macaulay2,
i11 : listPlot(L,Join=>false,ViewPort=>{vector{-2,-2},vector{2,2}})
o11 = x
y
o11 : GraphicsList
4.4. LLL part 1. Reproduce the LLL calculation from Example 5 in the lecture notes.
In Example 5, check that the first column b1 is indeed a shortest nonzero vector in the
lattice.
Solution.
i12 : M = matrix { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 }, { 785398, 197395, 4184 } };
o12 : Matrix Z4 ←− Z3
i13 : M’ = LLL M
o13 =

1 −13 −52
−4 58 203
1 −296 184
2 272 345

o13 : Matrix Z4 ←− Z3
First a sanity check: the lattices generated by M and M ′ are the same. In Macaulay2, use
image to find the image (i.e., span of column vectors) of a linear map:
i14 : image M == image M’
o14 = true
(Macaulay2 automatically understands “span” as “Z-span” because our matrices are defined
over Z).
32 P. ZINN-JUSTIN
Next we’re going to use the previous exercise to compute a bunch of lattice points, and
check their norm. It’s better to use M ′ rather than M for this search:
i15 : L = lattice3d(M’_0,M’_1,M’_2,-10..10,-10..10,-10..10);
Note the final semi-colon since we probably don’t want to output to the screen a list of
∼ 8000 vectors.
Finally, we ask for the smallest norm of nonzero vectors:
i16 : nrm = v -> sum(entries v,x->x^2);
i17 : min apply(select(L,x->x != 0_(ZZ^4)),nrm)
o17 = 22
and we compare to the norm of the first column of M ′:
i18 : nrm M’_0
o18 = 22
Now check that the second column b2 is a shortest lattice element not in SpanZ(b1).
Solution. It’s easy to select vectors that are not proportional to b1:
i19 : min apply(select(L,x -> x != x_0*M’_0),nrm)
o19 = 165 133
and sure enough
i20 : nrm M’_1
o20 = 165 133
Is the third column b3 a shortest lattice element not in SpanZ(b1,b2)?
Solution. We’re not going to be too clever – there are better ways of testing sublattice
membership:
i21 : L2 = lattice2d(M’_0,M’_1,-10..10,-10..10);
i22 : min apply(select(L,x -> not member(x,L2)),nrm)
o22 = 196 794
i23 : nrm M’_2
o23 = 196 794
4.5. LLL part 2. To get a better sense of how the LLL algorithm works, we’re going to reimplement
it from scratch. Here’s the pseudo-code version:
Input: {b1,b2, . . . ,bn} (as a matrix)
Repeat two steps until LLL reduced basis found:
Step 1: Gram--Schmidt-like orthogonalization
for i = 1 to n do
for j = i− 1 to 1 do
m ← nearest integer of µij = bi·b
∗
j
b∗j ·b∗j
bi ← bi −mbj
Step 2: Check condition 2 and swap
for i = 2 to n do
EXPERIMENTAL MATHEMATICS 2022 33
if ‖b∗i + µi,i−1b∗i−1‖2 < 34‖b∗i−1‖2 then
swap bi−1 and bi
go back to step 1
Output {b1,b2, . . . ,bn} (as a matrix)
Note that Macaulay2 distinguishes between matrices and nested lists: one can go back
and forth between the two using matrix and entries.
A suggestion of how to implement Gramm-Schmidt orthogonalization, which is needed for
step 1, is given in tutorial3-GS.m2:
i24 : -- Gram-Schmidt orthogonalization
sc = (a,b) -> sum(a,b,times);
i25 : GS = M -> (
b := new MutableList from transpose entries M;
n := #b;
for i from 0 to n-1 do for j from 0 to i-1 do (
mu := sc(b#i,b#j)/sc(b#j,b#j);
b#i = b#i - mu * b#j;
);
matrix transpose toList b
)
o25 = GS
o25 : FunctionClosure
Check that your implementation produces the same result as the built-in function on the
example of the lectures.
Solution.
i26 : myLLL = M -> (
b := new MutableList from transpose entries M;
n := #b;
while true do (
G := GS M;
bs := transpose entries G;
for i from 0 to n-1 do
for jj from 0 to i-1 do (
j := i-1-jj;
mu := sc(b#i,bs#j)/sc(bs#j,bs#j);
b#i = b#i - round(mu) * b#j;
);
i := 1;
while iv := bs#i+mu*bs#(i-1);
sc(v,v)>=3/4*sc(bs#(i-1),bs#(i-1))) do i=i+1;
if itemp := b#i;
b#i = b#(i-1);
b#(i-1) = temp;
34 P. ZINN-JUSTIN
);
M = matrix transpose toList b;
if i==n then return M;
)
)
o26 = myLLL
o26 : FunctionClosure
i27 :
. . . and . . .
myLLL M
o27 =

1 −13 −52
−4 58 203
1 −296 184
2 272 345

o27 : Matrix Z4 ←− Z3
4.6. Recognizing algebraic numbers. (if time permits – it didn’t)
An algebraic integer is a root α of a monic polynomial with integer coefficients:
f(α) = 0, f(x) = xn + an−1xn−1 + · · ·+ a1x+ a0, ai ∈ Z.
Given an algebraic integer α, its minimal polynomial is the monic integer polynomial f of
smallest degree d such that f(α) = 0. We refer to d as the degree of α.
For instance, α =
√
2 has minimal polynomial x2 − 2.
Suppose you have a real number α and you think it might be an algebraic integer of degree
d. Reduce the problem of finding the minimal polynomial of α to an instance of the integer
relation problem.
Apply your strategy from the previous exercise to find the minimal polynomial for
α = 1 +
√
1 +
√
1 +
√
2,
which seems like it might have degree ≤ 8, doesn’t it?
Let H = {z ∈ C | =(z) > 0}. The Dedekind eta function η : H → C is defined by
η(z) = q1/24
∞∏
n=1
(1− qn), q = e2piiz
Consider the number
α =
η2
(√−5/2)
2η2
(
2
√−5)
Figure out whether α may be algebraic. If yes, what is its minimal polynomial (and
degree)?
Figure out whether pi (yes, the pi) is algebraic. If yes, what is its minimal polynomial
(and degree)?
EXPERIMENTAL MATHEMATICS 2022 35
5. Polynomial equations and ideals
From linear algebra, we know how to solve a system of linear equations by doing Gaussian
elimination (aka row reduction). But what can we do if we are faced with a system of
nonlinear polynomial equations such as
x2 + y2 − 1 = 0
x2 − y = 0
Okay, that’s not a very good example since we recognise the first equation as describing the
unit circle and the second as a parabola and it is easy to do a substitution and solve for the
two intersection points.
Exercise 2. Do this. (Find the real solutions (x, y) of this system.)
Does anything change if I ask for complex solutions?
Such ad hoc methods will not always be available, for instance given
6x2y − x+ 4y3 − 1 = 0
2xy + y3 = 0
There is a systematic approach to simplifying such systems with a view towards under-
standing their solutions. This approach is based on the notion of Gro¨bner basis, which we
will encounter in the next lecture. For now, we need to describe/review a bit of algebra.
5.1. Rings and ideals. Just as vector spaces are the natural environment for reasoning
about systems of linear equations, (multivariate) polynomial rings are the natural environ-
ment for systems of polynomial equations.
The first ingredient is a base field (or field of coefficients) k. For the purposes of our
discussion you can imagine this to be whatever field you prefer, e.g. Q or R or C. Then we
form polynomial rings like
k[x1, . . . , xn] =
{∑
ca1,...,anx
a1
1 . . . x
an
n | a1, . . . , an ∈ Z≥0, ca1,...,an ∈ k
}
,
where the sums have only finitely many terms.
We will refer to the elements of k[x1, . . . , xn] as polynomials and perform the arithmetic
operations of addition and multiplication in the usual way. We might sometimes call them
formal polynomials as opposed to the polynomial functions that you encounter in calculus
(broadly construed). The link between the two is given by the operation of evaluation.
For any b = (b1, . . . , bn) ∈ kn we have evb : k[x1, . . . , xn] → k given by the substitution
x1 7→ b1, . . . , xn 7→ nn. This allows us to think of f ∈ k[x1, . . . , xn] as (giving rise to) a
function from kn → k defined by (b1, . . . , bn) 7→ evb(f). While it is important to have an
awareness of the distinction between polynomial and polynomial function, we’ll now join the
rest of the planet in identifying them, thus writing
f(b1, . . . , bn) = evb(f).
(such a confusion is harmless when the base field is infinite, as distinct polynomials give rise
to distinct polynomial functions).
We’ve now encountered polynomial rings, which are an important class of commutative
rings (here commutative refers to multiplication; addition is always supposed to be com-
mutative!). What about ideals? They occur naturally when thinking about systems of
polynomial equations, as we shall see now.
36 P. ZINN-JUSTIN
Going back to the example above, we can write the set of solutions of the system of its
equations as
V(f1, f2) =
{
(b1, b2) ∈ R2 | f1(b1, b2) = 0, f2(b1, b2) = 0
}
where f1(x, y) = 6x
2y−x+4y3−1 ∈ R[x, y], f2(x, y) = 2xy+y3 = 0 ∈ R[x, y], and V stands
for vanishing set.
We can define similarly the vanishing set of any subset of polynomials A ⊂ k[x1, . . . , xn]:
V(A) = {(b1, . . . , bn) ∈ kn | f(b1, . . . , bn) = 0 for all f ∈ A}.
And we can turn this idea on its head and define, for any subset B ⊂ kn, the set of
polynomials that vanish identically on B:
I(B) = {f ∈ k[x1, . . . , xn] | f(b1, . . . , bn) = 0 for all (b1, . . . , bn) ∈ B}.
The interplay between these two functions V and I is the starting point of algebraic geometry
and a beautiful piece of mathematics.
But I digress. The subset I := I(B) ⊂ R := k[x1, . . . , xn] exhibits some interesting extra
structure:
• if f1 and f2 are in I, then f1 + f2 ∈ I
• if f is in I and r is in R, then rf ∈ I
A nonempty subset I of a (commutative) ring R with these properties is called an ideal. (In
particular, an ideal always contains 0, so {0} is the smallest ideal.)
A great way of getting lots of ideals is choosing some elements f1, . . . , fm ∈ R and com-
bining them as follows:
〈f1, . . . , fm〉 = {r1f1 + . . . rmfm | r1, . . . , rm ∈ R}.
(Compare this to the span of vectors f1, . . . , fm. It is similar in shape but definitely not the
same thing.)
Exercise 3. Prove that 〈f1, . . . , fm〉 is indeed an ideal of R. It is called the ideal generated
by f1, . . . , fm.
Going back to a system such as
6x2y − x+ 4y3 − 1 = 0
2xy + y3 = 0
our strategy will be to form the ideal I = 〈f1, f2〉 of R[x, y] with f1 = 6x2y−x+ 4y3− 1 and
f2 = 2xy + y
3 and try to gain a better understanding of this ideal. More precisely, we will
aim to find a simpler set of generators for this ideal, one that will illuminate the solutions
of the system.
In Macaulay2, here’s how to define polynomial rings and ideals: (we’ve chosen the base
field to be Q in order for the results to be exact, rather than a numerical approximation;
the methods that we discuss later do not strongly depend on the base field)
i1 : R = QQ[x,y]
o1 = R
o1 : PolynomialRing
i2 : f1 = 6*x^2*y-x+4*y^3-1
o2 = 6x2y + 4 y3 − x− 1
o2 : R
i3 : f2 = 2*x*y^3
EXPERIMENTAL MATHEMATICS 2022 37
o3 = 2x y3
o3 : R
i4 : I = ideal(f1,f2)
o4 = ideal (6x2y + 4 y3 − x− 1, 2x y3)
o4 : Ideal of R
To get a feel for this in more familiar settings, we consider two very special cases of the
problem.
5.2. Linear systems. We go back to the linear case momentarily, because we know how to
deal with it. Consider the ideal I = 〈x+ y − z, 2x+ 3y + 2z〉 in R[x, y, z]. Its vanishing set
is the set of solutions of the linear system
x+ y − z = 0
2x+ 3y + 2z = 0
In order to solve this system, we apply Gaussian elimination to the associated matrix[
1 1 −1
2 3 2
]
→
[
1 1 −1
0 1 4
]
which tells us that I = 〈x+y−z, y+4z〉, an arguably simpler set of generators than the one
we started with. We could even go all in with Gauss–Jordan elimination to get the reduced
row-echelon form [
1 0 −5
0 1 4
]
telling us that I = 〈x− 5z, y + 4z〉.
From these generators we see that we can express both x and y in terms of z, so the
solution set is one-dimensional and we can write down an explicit parametrisation in terms
of z.
So in this special case we used elementary row operations in order to simplify the generators
of I. Note that we implicitly ordered our variables so that x comes before y comes before z.
Had we decided to order them differently, say z before y before x, the same method
applies but we would get the differently-looking answer I = 〈−z + y + x, 2z + 3y + 2x〉 =
〈z − x/5, y + 4x/5〉. This is a feature of the whole multivariate polyonomial business, and
we’ll have to think about it more carefully soon.
Macaulay2 implements such simplifications automatically with trim:
i5 : R=QQ[x,y,z]
o5 = R
o5 : PolynomialRing
i6 : I=ideal(x+y-z,2*x+3*y+2*z)
o6 = ideal (x+ y − z, 2x+ 3 y + 2 z)
o6 : Ideal of R
i7 : trim I
o7 = ideal (y + 4 z, x− 5 z)
o7 : Ideal of R
38 P. ZINN-JUSTIN
5.3. Single variable higher degree case. The other special situation we consider is that
of higher-degree polynomials in a single variable x. In other words, we now work with
arbitrary elements of k[x].
Let’s take the ideal I = 〈x3− 2x2 + 2x+ 8, 2x2 + 3x+ 1〉. It’s not clear what the vanishing
set of I might be, so we try to find a simpler set of generators.
We can do this by polynomial long division:
x3 − 2x2 + 2x+ 8 =
(
1
2
x− 7
4
)(
2x2 + 3x+ 1
)
+
(
27
4
x+
39
4
)
,
from which we conclude that I =
〈
2x2 + 3x+ 1, 27
4
x+ 39
4
〉
. This is simpler than before (the
generators have degrees 2 and 1 instead of the original 3 and 2), but still not as simple as
could be3.
We repeat the process:
2x2 + 3x+ 1 =
(
8
27
x+
4
243
)(
27
4
x+
39
4
)
+
68
81
,
so I =
〈
27
4
x+ 39
4
, 68
81
〉
. We could do another long division of the first generator by the
second, but that’s overkill since we already have the condition 68
81
= 0, which indicates that
the system has no solutions.
Let’s check in Macaulay2:
i8 : R=QQ[x]
o8 = R
o8 : PolynomialRing
i9 : I=ideal(x^3-2*x^2+2*x+8,2*x^2+3*x+1)
o9 = ideal (x3 − 2x2 + 2x+ 8, 2x2 + 3x+ 1)
o9 : Ideal of R
i10 : trim I
o10 = ideal 1
o10 : Ideal of R
Note that there is no ambiguity of ordering of the variables here, as there is a single
variable.
The Euclidean algorithm for k[x] (which is what this process of repeated long division is
called) shows that every ideal in k[x] can be generated by a single element. This ceases to
be the case for the general situation k[x1, . . . , xn] (in the tutorial you’ll be asked to think of
an ideal in there that cannot be generated by a single element), but something useful is still
true:
Theorem 2 (Hilbert Basis Theorem). Every ideal in k[x1, . . . , xn] can be generated by a
finite set of elements.
Having considered the two special cases of this section and the previous one, we turn our
attention to the general case of an ideal (or system of equations) in k[x1, . . . , xn]. Before
we investigate the common generalisation of Gaussian elimination and of single-variable
polynomial long division, let’s note that many questions about the solution set of a system
can be expressed in terms of the associated ideal.
3You can however use the degree 1 generator to solve for x and then plug the solution into the other
generator to reach the same conclusion as in the next paragraph.
EXPERIMENTAL MATHEMATICS 2022 39
For instance, if k is an algebraically closed field (e.g. C) and I is the ideal associated
with a system of equations in k[x1, . . . , xn], then the system has no solutions if and only if
I contains the constant 1. (As we observed in the example of this section, one direction of
this claim is very easy and holds without any conditions on the field k.)
For the remainder of this chapter, we will focus on the ideals themselves rather than the
systems of equations that they may have originated from.
40 P. ZINN-JUSTIN
References
[Cox84] David A. Cox, The arithmetic-geometric mean of Gauss, Enseign. Math. (2) 30 (1984), no. 3-4,
275–330. MR767905.
[Cox85] , Gauss and the arithmetic-geometric mean, Notices Amer. Math. Soc. 32 (1985), no. 2,
147–151. MR779224.
[LLL82] A. K. Lenstra, H. W. Lenstra, Jr., and L. Lova´sz, Factoring polynomials with rational coefficients,
Math. Ann. 261 (1982), no. 4, 515–534, doi:10.1007/BF01457454. MR682664.
[LS10] Joel C. Langer and David A. Singer, Reflections on the lemniscate of Bernoulli: the forty-eight faces
of a mathematical gem, Milan J. Math. 78 (2010), no. 2, 643–682, https://case.edu/artsci/
math/langer/jlpreprints/fortyeight2010.pdf. MR2781856.
[Sal76] Eugene Salamin, Computation of pi using arithmetic-geometric mean, Math. Comp. 30 (1976),
no. 135, 565–570, doi:10.2307/2005327. MR404124.
[Sta12] Richard P. Stanley, Enumerative combinatorics. Volume 1, second ed., Cambridge Studies in Ad-
vanced Mathematics, vol. 49, Cambridge University Press, Cambridge, 2012. MR2868112.