Mathematics, read from a classifier's source

Co-authored by Claude · Version 1.0.0, changed in rev 133

Some time between 1800 and 1600 BC a student in southern Mesopotamia pressed a square and its two diagonals into a clay tablet, now catalogued as YBC 7289. Along one diagonal are two numbers in base sixty. The first, 1;24,51,10, is the ratio of a diagonal to its side, the square root of 2, to about 1.41421296, right to six decimal places. The second, 42;25,35, is the same ratio times 30, the length written on the side. The tablet is a calculation done and kept, and a good deal of mathematics has been that: a procedure fixed so that the next person gets the answer without redoing the reasoning.

The mathematics library for DKE, the Deterministic Knowledge Engine, writes some of that mathematics down in DKE Python, and this essay teaches the mathematics by reading it. You record three kinds of fact. A whole number is written with a row for every candidate divisor, and the library's rules conclude which divide it, the sum of its proper divisors, and whether it is prime or composite and perfect, abundant or deficient. A binary operation on a small finite set is written as its table, and the rules conclude its identity, whether it is closed, commutative and associative, which elements have inverses and whether it is a group. A triangle is written as its three sides in whole units, and the rules conclude whether the lengths make a triangle, and whether it is acute, right or obtuse and equilateral, isosceles or scalene. Each answer names the facts it rests on, and correcting a fact withdraws the answers it supported. The second half of the library calculates: Euclid's algorithm and what is built on it, primes and factors, counting, continued fractions and Pell's equation, roots, integrals and linear systems. Every program below runs as shown once the library's source file has been compiled into a store, and every output is copied from a run.

1. Classifying the whole numbers

A number is prime if nothing between 1 and itself divides it, and composite if something does. Number theorists sort by a second question too: add up a number's proper divisors, those that are not the number itself, and see whether the sum is equal to the number, more, or less. Equal, and the number is perfect; more, abundant; less, deficient. The library writes each number with a candidate for every possible divisor up to its half, since no proper divisor is larger, and lets the rules do the sorting:

@rule
def divisor_sum():
    for x in Number, d in Candidate:
        if d.of == x and d.divides == True and absent(x.disputed):
            x.divisor_sum = sum(d.k)

if x.n > 1 and none(d in Candidate: d.of == x and d.k > 1 and d.divides == True):
    x.prime = True

The first rule's last premise holds the sum back while a candidate is bound to its number twice, as a second source writing a candidate's k, even the same k, would bind it. A sum would add it twice, and 6 would read as abundant, which no source says; numbers_that_are("disputed") names such a number instead.

Here are the numbers from 1 to 60. The candidates grow as the square of the largest number, so this is for numbers in the hundreds and not the millions:

import mathematics

mathematics.numbers(1, 60)
print(mathematics.numbers_that_are("prime"))
print(mathematics.numbers_that_are("perfect"), mathematics.divisors(28))
print(mathematics.numbers_that_are("abundant"))
print(len(mathematics.numbers_that_are("deficient")), len(mathematics.numbers_that_are("composite")))

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59]
[6, 28] [1, 2, 4, 7, 14]
[12, 18, 20, 24, 30, 36, 40, 42, 48, 54, 56, 60]
46 42

There are seventeen primes. The perfect numbers are 6 and 28, and 28 is 1 + 2 + 4 + 7 + 14. The abundant numbers begin with 12, and the library has them to 60. Euclid showed, in Proposition 36 of Book IX of the Elements, written about 300 BC, that if 2^p − 1 is prime then 2^(p−1) times it is perfect: 6 is 2 times 3 and 28 is 4 times 7. Nicomachus noted 8128 about AD 100 as the fourth. Euler later proved that every even perfect number has Euclid's form. Nobody knows whether an odd one exists.

A conclusion answers why, and what failed answers why_not. 21 is composite because 3 and 7 divide it. 7 is not, and the answer names the rule's conditions and says that each is met on its own but never together by one candidate:

import mathematics

w = why(mathematics.Number.n21.composite)
for c in w:
    print(str(c.path) + " = " + str(c.value))
reason = why_not(mathematics.Number.n7.composite)
print(reason)

mathematics.Number.n21.composite = True
mathematics.Candidate.n21_3.of = n21
mathematics.Candidate.n21_3.k = 3
mathematics.Candidate.n21_3.divides = True
mathematics.Candidate.n21_3.remainder = 0
mathematics.Candidate.n21_3.n = 21
mathematics.Candidate.n21_7.of = n21
mathematics.Candidate.n21_7.k = 7
mathematics.Candidate.n21_7.divides = True
mathematics.Candidate.n21_7.remainder = 0
mathematics.Candidate.n21_7.n = 21
rule `mathematics.composite` needs mathematics.Candidate[`d`].of == "n7", mathematics.Candidate[`d`].k > 1 and mathematics.Candidate[`d`].divides == True together: each is met on its own, but no one `d` meets all of the rule's conditions at once

2. Operations with a table

An operation on a finite set can be written as a table whose row a and column b hold a · b, and the properties mathematicians care about are properties of the table. The operation is closed if every entry is in the set, commutative if the table is symmetric about its diagonal, and associative if (x · y) · z and x · (y · z) never differ. It has an identity if some element leaves every element alone from either side, and an element has an inverse if there is another that the operation takes to the identity, from both sides. An operation that is closed and associative, with an identity that every element has an inverse for, is a group:

if u.op == t and u.identity == True and t.closed == True and t.associative == True and none(e in Element: e.op == t and e.noninvertible == True):
    t.group = True

Évariste Galois was the first, in the 1830s, to use groups to settle which polynomial equations can be solved, and he gave them the name. Abstract groups took hold late in the nineteenth century. The tables here are the clock arithmetic of four hours, with addition, subtraction and multiplication, and each is closed:

import mathematics

els = ["0", "1", "2", "3"]
add = [[str((i + j) % 4) for j in range(4)] for i in range(4)]
sub = [[str((i - j) % 4) for j in range(4)] for i in range(4)]
mul = [[str((i * j) % 4) for j in range(4)] for i in range(4)]
mathematics.operation("add", els, add)
mathematics.operation("sub", els, sub)
mathematics.operation("mul", els, mul)
print(mathematics.properties("add"))
print(mathematics.properties("sub"))
print(mathematics.properties("mul"))
print(mathematics.nonassociative_triples("sub")[0])
print(mathematics.identity_of("mul"), mathematics.without_inverse("mul"))

['closed', 'commutative', 'associative', 'group', 'abelian']
['closed']
['closed', 'commutative', 'associative']
0 0 1
1 ['0', '2']

Addition is a group, and an abelian one, since its table is symmetric. Subtraction is closed and nothing else: the triple 0, 0, 1 shows it, for (0 − 0) − 1 is 3 and 0 − (0 − 1) is 1. Multiplication is closed, commutative and associative and has the identity 1, but is not a group, because 0 and 2 have no inverse: nothing times 2 is 1 on a clock of four.

3. A table that does not commute

The six symmetries of an equilateral triangle, the three rotations and the three reflections, are a group of order six, and the smallest group in which the order of two actions matters. Rotate it and then flip it, and you do not get what flipping first and rotating second gives. Here is its table, and the first pair the library finds that do not commute. Then one cell is rewritten, and the properties that cell held up go with it:

import mathematics

sym = ["e", "r", "r2", "s", "sr", "sr2"]
table = [["e", "r", "r2", "s", "sr", "sr2"], ["r", "r2", "e", "sr", "sr2", "s"], ["r2", "e", "r", "sr2", "s", "sr"], ["s", "sr2", "sr", "e", "r2", "r"], ["sr", "s", "sr2", "r", "e", "r2"], ["sr2", "sr", "s", "r2", "r", "e"]]
mathematics.operation("triangle", sym, table)
print(mathematics.properties("triangle"))
print(mathematics.noncommuting_pairs("triangle")[0])
mathematics.cell("triangle", "r", "s", "sr2")
print(mathematics.properties("triangle"))
print(mathematics.unsettled())

['closed', 'associative', 'group']
r s
['closed']
[]

The table is closed and associative and a group, and not commutative. The library names the pair r and s, a rotation and a reflection. Rewriting the cell r · s breaks the associativity, and the group goes with it, leaving a table that is only closed. No other source wrote the cell, so this is a correction and not a disagreement, and unsettled() has nothing to report.

4. Triangles

Euclid's Proposition 20 of Book I says that in any triangle the sum of any two sides is greater than the remaining one. The library's triangle is written by its sides in whole units, sorted so that a is at most b at most c, and it computes two numbers from them. slack is a + b − c, positive exactly when the lengths close up, and gap is *a*² + *b*² − *c*², whose sign is the cosine of the angle facing the longest side:

gap: int = a * a + b * b - c * c

slack: int = a + b - c

The rules read the sign: zero, a right angle; positive, acute; negative, obtuse. A triangle's area follows from Heron's formula, which Heron of Alexandria set down in his Metrica about AD 60 and which Archimedes may have known two centuries before:

import mathematics

mathematics.triangle("egypt", 3, 4, 5)
mathematics.triangle("flat", 1, 2, 3)
mathematics.triangle("none", 1, 2, 5)
mathematics.triangle("equi", 5, 5, 5)
mathematics.triangle("wide", 2, 3, 4)
mathematics.triangle("tall", 5, 5, 8)
for n in ["egypt", "flat", "none", "equi", "wide", "tall"]:
    print(n, mathematics.kinds(n), mathematics.area(n))

egypt ['right', 'scalene'] 6.0
flat [] 0.0
none [] 0.0
equi ['acute', 'equilateral', 'isosceles'] 10.825317547305483
wide ['obtuse', 'scalene'] 2.9047375096555625
tall ['obtuse', 'isosceles'] 12.0

The 3, 4, 5 triangle is right-angled and scalene, with an area of 6. 1, 2, 3 has no slack and is a flat line, and 1, 2, 5 would never close; neither is a triangle, and both have an empty list of kinds and an area of 0. The equilateral triangle is also isosceles and acute. 2, 3, 4 is obtuse, and 5, 5, 8 is obtuse and isosceles.

The Babylonian tablet Plimpton 322, about 1800 BC, has fifteen rows each relating a right triangle with whole sides. The rule usually called Euclid's formula makes such triples from two numbers m and n with m greater than n, coprime and of opposite parity: the sides are *m*² − *n*², 2*mn* and *m*² + *n*². The next program makes the triangles that rule gives and asks the library what they are:

import mathematics

for m in range(2, 6):
    for n in range(1, m):
        if mathematics.gcd(m, n) == 1 and (m - n) % 2 == 1:
            name = "euclid_" + str(m) + "_" + str(n)
            mathematics.triangle(name, m * m - n * n, 2 * m * n, m * m + n * n)
            print(m, n, m * m - n * n, 2 * m * n, m * m + n * n, mathematics.kinds(name))

2 1 3 4 5 ['right', 'scalene']
3 2 5 12 13 ['right', 'scalene']
4 1 15 8 17 ['right', 'scalene']
4 3 7 24 25 ['right', 'scalene']
5 2 21 20 29 ['right', 'scalene']
5 4 9 40 41 ['right', 'scalene']

Each is right-angled, which is a test of Euclid's formula against the rule that decides it from the sides and says nothing of why the formula works.

5. Euclid's algorithm

To find the greatest common divisor of two numbers, divide the larger by the smaller and replace the larger by the remainder, until the remainder is nothing. Euclid sets it out in Book VII of the Elements, about 300 BC, and it was probably known earlier. Donald Knuth called it the granddaddy of all algorithms, because it is the oldest nontrivial one that survives. The heart of the library's gcd() is its remainder step:

r = x % y

The extended form keeps track of how each remainder is made from the two numbers, and ends with the pair that satisfy a·x + b·y = the gcd, which makes inverses in modular arithmetic. In 1977 Rivest, Shamir and Adleman published a public-key system built on them; Clifford Cocks at GCHQ had described a similar one in 1973, which was declassified in 1997. The small example with the primes 61 and 53 that the textbooks use goes through the library from end to end. The modulus is their product, 3233, and the number of residues coprime to it is the product of the primes less one, 3120. The key e = 17 has an inverse modulo 3120, 2753, and a message of 65 is raised to the 17th power modulo 3233 and recovered with the 2753rd:

import mathematics

print(mathematics.gcd(1071, 462), mathematics.lcm(21, 6))
print(mathematics.extended_gcd(240, 46))
print(mathematics.modinv(17, 3120), mathematics.modinv(6, 9))
n = 61 * 53
phi = 60 * 52
print(n, mathematics.totient(n), phi)
c = mathematics.modpow(65, 17, n)
print(c, mathematics.modpow(c, 2753, n))

21 42
[2, -9, 47]
2753 0
3233 3120 3120
2790 65

gcd(1071, 462) is 21, the example found in most accounts of the algorithm. The extended form of 240 and 46 gives 2, −9 and 47, for 240 × (−9) + 46 × 47 = 2. modinv() finds the inverse of 17 modulo 3120 as 2753 and gives 0 for 6 modulo 9, because 6 and 9 share a factor and no inverse exists. totient() counts the residues coprime to 3233 and agrees with the 3120 computed by hand. The message 65 goes to 2790 and comes back as 65. The same example is often written with 413 for the secret exponent, which works with the lcm of 60 and 52, 780, in place of 3120; 2753 is 413 plus three times 780, so the two are the same key.

6. Primes and their factors

The sieve of Eratosthenes, which Nicomachus attributes to a mathematician of the third century BC, finds every prime up to a limit by crossing off the multiples of each prime as it is found. Testing one number is another problem. The library's is_prime() uses the test Gary Miller devised in 1976 and Michael Rabin made unconditional in 1980, with the fixed bases 2, 3, 5 and 7, which are enough to decide every number below 3,215,031,751. factorize() divides by trial and isqrt() finds the largest whole number whose square is at most n:

import mathematics

print(mathematics.primes(100))
print(mathematics.is_prime(561), mathematics.factorize(561))
print(mathematics.is_prime(2147483647), mathematics.is_prime(2147483649))
print(mathematics.factorize(4294967297))
print(mathematics.factorize(360), mathematics.totient(360))
print(mathematics.isqrt(1000000007), mathematics.isqrt(9223372036854775807))

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
False [3, 11, 17]
True False
[641, 6700417]
[2, 2, 2, 3, 3, 5] 96
31622 3037000499

561 is the smallest Carmichael number: it is composite, 3 × 11 × 17, and it passes the older Fermat test to every base coprime to it, as every such number does. Václav Šimerka found the first seven in 1885, Robert Carmichael published 561 in 1910, and Alford, Granville and Pomerance proved in 1994 that there are infinitely many. No composite number passes Miller and Rabin's test to every base, so the library's test says False for 561. 2147483647 is the Mersenne prime 2^31 − 1, which Euler proved prime in 1772, when it was the largest prime known. 4294967297 is 2^32 + 1: Fermat conjectured that numbers of that form are all prime, and Euler showed in 1732 that this one is 641 × 6,700,417. The last line is isqrt() of 1,000,000,007 and of the largest 64-bit integer, 3,037,000,499, which is the bound of modpow() and is_prime().

7. Counting, and where an integer ends

An int is 64 bits, and the library refuses a result past 9,223,372,036,854,775,807 rather than wrapping. 20! is the largest factorial that fits, and the library refuses 21!. The Fibonacci numbers, which Fibonacci brought into Western mathematics in the Liber Abaci of 1202 from sources in India, stop at the 92nd, and it refuses the 93rd. binomial() divides out common factors at each step, so it reaches any coefficient that fits and not only those whose factorials do. And collatz() counts the steps for a number to reach 1 when it is halved if even and tripled and increased by one if odd:

import mathematics

print(mathematics.factorial(20), mathematics.fibonacci(92))
print(mathematics.binomial(10, 3), mathematics.binomial(60, 30))
print(mathematics.collatz(27), mathematics.collatz(97), mathematics.collatz(6171))

2432902008176640000 7540113804746346429
120 118264581564861424
111 118 261

The Collatz question was posed by Lothar Collatz in 1937, and the conjecture that every start reaches 1 has been checked up to 2^71 and not proved. 27 takes 111 steps, climbing to 9,232 on the way, and 97 takes 118 and 6171 takes 261. The library gives up after 1,024 steps and returns −1.

8. Continued fractions and Pell's equation

Euclid's algorithm applied to a fraction produces its continued fraction, the quotients it finds. Applied to pi to fifteen places, it gives 3, 7, 15, 1, 292, and the fractions made from the first few terms, the convergents, are the best rational approximations there are: 22/7, 333/106, 355/113 and 103993/33102. Zu Chongzhi in the fifth century gave 22/7 and 355/113 as approximations. The 292 is why 355/113 is so good: it is right to seven places, with a numerator and denominator that take six digits to write. The square roots of whole numbers have continued fractions that repeat, and Pell's equation, *x*² − n·*y*² = 1, is solved by one of them:

import mathematics

print(mathematics.continued_fraction(3141592653589793, 1000000000000000))
terms = mathematics.continued_fraction(3141592653589793, 1000000000000000)
c = mathematics.convergents(terms)
print(c[0], c[1], c[2], c[3], c[4])
print(mathematics.sqrt_period(2), mathematics.sqrt_period(61))
print(mathematics.pell(2), mathematics.pell(61))
p = mathematics.pell(61)
print(p[0] * p[0] - 61 * p[1] * p[1])

[3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 4, 2, 3, 1, 12, 5, 1, 5, 20, 1, 11, 1, 1, 1, 2]
[3, 1] [22, 7] [333, 106] [355, 113] [103993, 33102]
[1, 2] [7, 1, 4, 3, 1, 2, 2, 1, 3, 4, 1, 14]
[3, 2] [1766319049, 226153980]
1

The continued fraction of the square root of 2 is 1 followed by 2 for ever, and that of 61 is 7 followed by a repeating part of eleven terms, ending in 14, which is twice 7. The smallest solution of *x*² − 61*y*² = 1 is 1,766,319,049 and 226,153,980; the last line checks it, and gets 1. Bhaskara II found it in 1150 with a general method, and Fermat in 1657 posed such equations as a challenge to English mathematicians, of whom Brouncker is credited with the first European solution. Euler misattributed it to John Pell, and the name stuck. Lagrange developed the theory of its relation to continued fractions in 1766-1769.

9. Roots, integrals and pi

How the Babylonian student found the square root of 2 is not recorded. A rule that converges to it, called the Babylonian method and now seen as a case of Newton's method, takes a guess g and replaces it by the mean of g and 2/g. root() does that for any whole root, from above and stopping when a step no longer brings it down. trapezoid() and simpson() integrate from samples at equal steps. Thomas Simpson's rule of 1743 weights the samples one, four, two, four, and so on, and a similar rule was derived by Kepler in 1615 for measuring wine barrels. The integral of the sine from 0 to pi is 2. machin() is John Machin's formula of 1706, with which he calculated 100 digits of pi:

import mathematics
import math

print(mathematics.root(2.0, 2), math.sqrt(2.0))
print(mathematics.root(2.0, 3), mathematics.root(1000000.0, 6))
h = math.pi / 10
ys = [math.sin(i * h) for i in range(11)]
print(mathematics.trapezoid(ys, h), mathematics.simpson(ys, h))
print(mathematics.machin(), math.pi)

1.414213562373095 1.4142135623730951
1.2599210498948732 10.0
1.9835235375094546 2.0001095173150043
3.1415926535897936 3.141592653589793

The root agrees with the machine's square root to within a unit of the last place. From eleven samples a tenth of a pi apart, Simpson's rule gets 2.0001 and the trapezoid 1.9835. Machin's formula gives pi to within a unit of the last place.

10. Linear systems and the matrix that defeats them

Eliminating variables from a set of linear equations is in The Nine Chapters on the Mathematical Art, which is known by 179 AD and may be as old as 150 BC; the algorithm took Gauss's name only in the 1950s. determinant() and solve() use it with partial pivoting, choosing the row with the largest entry in each column, which is what keeps the arithmetic stable. solve() returns an empty list for a singular matrix. David Hilbert introduced his matrix in 1894, 1/(i + j + 1) in row i and column j from zero, and it is the standard example of an ill-conditioned one. The program below builds the Hilbert matrix of three sizes and a right-hand side that makes the exact answer all ones, solves it, and reports the determinant and the largest error:

import mathematics

print(mathematics.determinant([[2.0, 1.0], [1.0, 3.0]]), mathematics.solve([[2.0, 1.0], [1.0, 3.0]], [3.0, 5.0]))
print(mathematics.solve([[1.0, 2.0], [2.0, 4.0]], [1.0, 2.0]))
for n in [4, 8, 12]:
    h = mathematics.hilbert(n)
    b = [sum(h[i]) for i in range(n)]
    x = mathematics.solve(h, b)
    print(n, mathematics.determinant(h), max([abs(x[i] - 1.0) for i in range(n)]))

5.0 [0.8, 1.4]
[]
4 1.6534391534393e-07 4.558575739110893e-13
8 2.7370501431517162e-33 1.3720025349162057e-07
12 2.8332547472051984e-78 0.005954221523187808

The 2 by 2 system is solved exactly, and the matrix with a second row twice the first is singular. For the Hilbert matrix of 4 the answer is right to within a trillionth, for 8 the error is 1 × 10⁻⁷, and for 12 it is 0.006, with a determinant of 10⁻⁷⁸. The determinant falls toward zero faster than elimination can follow.

11. What the library leaves out

A number is a candidate list, so classification is for hundreds. A table is n² cells and n³ triples, so an operation is for a few elements, and nothing in the library says a group is cyclic, simple or isomorphic to another. A triangle has whole sides, and the library finds its kind and area and nothing of its angles. The calculations are of 64-bit integers and 64-bit reals, with the ceilings each script's comment names, and they are not proofs: is_prime() is the answer of a test that is right below its bound, and collatz() is a count, not a theorem. Numerical methods are the plain ones, with fixed steps and no error estimates. These are the textbook cases, and each has an answer known some other way.

Closing

Mathematics keeps its procedures, and the tablet of 1700 BC was one. The library keeps a few of them. It sorts the numbers by what divides them, reads the properties of an operation off its table and says which cell breaks one, tells a triangle from three lengths that cannot close, and runs Euclid's algorithm, which is older than any of the rest. And because each conclusion names what it rests on, a reader who doubts one can ask why, and follow the answer down to the divisor, the cell or the side it came from.