A Pythagorean triple is a triple of integers (a,b,c) forming the side lengths of a right triangle and satisfying the Pythagorean equation a² + b² = c². There's some evidence that these triples were known to the ancient Babylonians, for whom the regular numbers (having only 2, 3, and 5 as prime factors) were especially important. So, what are the regular Pythagorean triples? We might as well limit our attention to primitive triples (those having no common divisor) because we can form any other triple by multiplying all the numbers in a primitive triple by the same scalar.

Theorem: The only regular primitive Pythagorean triple is (3,4,5).

( Proof )

Because every Pythagorean triple includes numbers divisible by 2, 3, and 5, this shows more generally that (3,4,5) and its multiples are the only Pythagorean triples with only three prime factors. This also shows that there are only four primitive Pythagorean triples (a,b,c) for which a and b are regular, ignoring the regularity of c: (3,4,5), (5,12,13), (8,15,17), and (9,40,41). Partial results of Stewart and Tijdeman on the abc conjecture can be used to show that, for any set P of primes, there are only finitely many Pythagorean triples (a,b,c) for which a and b are both P-smooth, but I don't know of an efficient method for listing all such triples in general.

I received my copy of Knuth's Art of Computer Programming, Volume 4, Fascicle 0, this week. It contains a long section on median graphs, and coincidentally I had just added a long article on median graphs to Wikipedia, so I've been reading that part with great interest.

One of the exercises, 7.1.1 ex.109, caught my attention. In it, Knuth describes the "binary majorization lattice" in which the elements are the n-bit binary numbers and in which x is greater than y if x has at least as many nonzeros as y with x's highest order nonzero bit larger than y's, x's second highest nonzero bit larger than y's, and so on. Here it is as drawn by my algorithm for drawing graphs that can be embedded on an integer grid — here the integer grid turns out to be three dimensional, and the awkward folds in the drawing are there because my algorithm doesn't know any better than to put them in.

( Read more... )

Find an algorithm for counting the number of isosceles triangles in an n × n grid.

( Read more... )

( Added much simpler analysis )

The binomial coefficients forming Pascal's triangle count, among other things, the number of lower sets in a two-dimensional rectangular grid. What about three-dimensional grids?

( Read more... )

Mark-Jason Dominus writes about taking a too-sophisticated approach to a puzzler about generating triangular numbers with many divisors, ending up with a program that produced the answer in a minute and a half when someone else's much shorter code took longer to compute but was much easier to write. But it took me only five minutes or so (less than it took to write up this blog entry, anyway) to write an even more sophisticated version of Dominus' idea, that produces the correct answer instantly even running in Python. The trick: I had a generator for already-factored numbers prewritten in PADS and could just plug it in. And of course taking advantage of the analysis Mark wrote up for his solution kept my own thinking time low...here's the code, slightly cleaned up from my original version:

( Cut for Python code )

All pairs of consecutive 47-smooth numbers, calculated by my Python code for Lehmer's improvement to Størmer's theorem. That is, this file lists pairs x, x+1 such that all prime factors of both x and x+1 are at most 47. There are 1502 such pairs, the largest of which is formed by 1109496723125 = 5⁴ × 7 × 17 × 19² × 31² × 43 and 1109496723126 = 2 × 3 × 11² × 23 × 29² × 41² × 47.

It took nearly two days on a dual-core Pentium, so the next prime 53 may be in reach but would likely take more computer time than I'm willing to devote to it. Judging by A002072, faster codes are out there.

I doubt I'll be awake in time, but Karl Rubin will be speaking about Fermat's last theorem, tomorrow (Tuesday) morning at 7:30, at the Beckman Center.

Størmer's theorem is really not so much a theorem but an algorithm, for finding pairs of consecutive smooth numbers. I implemented Lehmer's improved algorithm, in Python, as I needed an implementation to generate some examples while expanding the WP article (and maybe more importantly, to check my understanding of Lehmer's paper). It includes some careful bignum implementation of rounding quadratic irrationals, finding continued fraction expansions of them, and solving Pell's equation, which might be useful for other purposes.

Some years before Dijkstra attributed it to Hamming, Gingerich considered an interesting variant of Hamming's problem: list, in order, only those regular numbers in the interval [60^k,60^k+1).

The motivation comes from the study of Babylonian mathematics. In the sexagesimal system the Babylonians used, which we still use a modified version of for times, latitudes and longitudes, and angles, the power of the initial base-60 digit was not specified. That is, if I tell you that a duration of time is 1:12, it's not obvious from that notation whether I mean one hour and twelve minutes (that is, 72 minutes), or one minute and twelve seconds (that is, 72 seconds); both are written the same way. So, in a Babylonian table of the reciprocals of regular numbers, one would not have separate entries for 72, 72×60, and 72×60×60; all three would just be written as 1:12. Grouping equivalent numbers in this way allowed the Babylonians a considerable space savings in their arithmetic tables: there are O(k²) different k-digit sexagesimal representations of regular numbers to be listed, while there are Θ(k³) actual regular numbers that can be represented with k digits. The sorted order of the regular numbers in Gingerich's range is the order in which the Babylonians would have listed these numbers.

Gingerich listed these numbers out of order, and then sorted them. Can we do better, and list them in sorted order in a constant number of arithmetic operations per output value? Yes!

( Read more... )

I just merged today two articles in Wikipedia on the same problem, regular numbers, that is, numbers of the form 2ⁱ3^j5^k. The ancient Babylonians liked these numbers because their reciprocals had terminating sexagesimal representations; more recently, "Hamming's problem" is the task of generating these numbers efficiently. Dijkstra described a very simple but inefficient functional solution: if H is the sequence of numbers to be generated, then H consists of the number 1 concatenated with the merger of 2H, 3H, and 5H. But there are two flaws with this. First, it generates the same numbers many times: e.g. 6 is generated both as 2*3 and as 3*2. This can be alleviated by letting A be the powers of 2, B be A merged with 3B, and H be B merged with 5H. But even if that improvement is made, and one does the merges in the most obvious and natural lazy-functional way, the sequence is generated inefficiently: the number of merges a generated number is involved in is O(j+k), and all these merges cause generating n numbers in this way to take Θ(n^4/3) arithmetic operations.

Below is an approach that instead generates n numbers in a linear number of operations. The basic principle is, when solving equations like H = merge(B, 5*H), instead of generating 5*H using a recursive instance of the same algorithm, to cache a single list of the whole sequence of H and reuse that stored sequence when computing 5*H.

( Cut for Python code )

However, the pure functional approach uses less space to generate the same sequence, O(n^2/3) instead of O(n) for the caching approach. I can find algorithms that generate the sequence in O(n^2/3) space and O(n log n) time (essentially, by turning it into a graph shortest path problem in an infinite grid graph, applying Dijkstra's algorithm for shortest paths, and only keeping track of nodes on the search frontier), but is it possible to solve optimally both in time and space?

ETA: Yes, it is, by a very simple change to the above code. Simply throw away the parts of the cached sequence that the back pointer has moved past. ( Read more... )

ETA2, March 14: Writing in a more imperative style actually gives tighter code, though perhaps it is not as readable: ( Read more... )

The interesting number paradox results from attempting to classify numbers as interesting or dull according to whether there is some simply-stated property that describes the number; for instance 255 is interesting (some might say) because it's the smallest perfect totient number with three different prime factors. The paradox comes from trying to find the value of the smallest number that is not interesting.

But: Wikipedia also has a description of notable numbers that seems essentially the same as the interesting numbers: "numbers with some remarkable mathematical property". There is at any point in time a smallest integer not deemed notable by the WP editors (as of today, that number is 202). And as "not deemed notable by the WP editors" is a property that is itself neither mathematical nor deemed remarkable by the WP editors, there is no paradox!

So is the disappearance of this paradox itself a paradox?

Today's submission to OEIS:

5, 13, 37, 109, 2917, 19131877, 57395629, 16210220612075905069

They are the primes of the form 4×3ⁿ+1, for n up to 39 (after which I got lazy and stopped, but it may be possible to find larger ones without specialized algorithms).

To my knowledge T. Venkataraman was the first to consider these primes; he showed in 1975 that multiplying them by 3 produces perfect totient numbers.

Venkataraman, T. (1975). "Perfect totient number". The Mathematics Student 43: 178. MR0447089.

ETA: The 1972 Math. Comp. reference in A005537 is a little older, but less specific, so I think the title of this post is still appropriate.

My recent subgreedy conjecture on Egyptian fractions implies, not just that the odd greedy algorithm terminates, but also that the corresponding even greedy algorithm (in which at each step we subtract the smallest even unit fraction from the remaining value to be expanded) produces finite expansions. But there's a much simpler proof for this case:

( Cut for spoiler )

The greedy algorithm for Egyptian fractions takes as input a fraction x/y, and repeatedly expands x/y = 1/d + x'/y', where d=ceiling(y/x) is the smallest possible denominator such that 1/d is at most x/y. Once one reaches a remaining fraction x/y such that (in its simplest form) x = 1, the algorithm terminates.

Analogously, define a subgreedy expansion of x/y to be the result of repeatedly expanding either x/y = 1/d + x'/y' as above, or x/y = 1/(d+1) + x'/y', where d=ceiling(y/x) as before. Also as before, the expansion terminates as soon as the remaining fraction has a unit numerator in its simplest form.

( Examples )

Conjecture. Any fraction x/y has a finite number of subgreedy expansions. More precisely for x/y < 1 the number of expansions is O(y).

The significance of this conjecture is that the odd greedy expansion is a subgreedy expansion. It's a well-studied open problem whether the odd greedy algorithm always terminates, but its termination would follow directly from this conjecture, because if the odd greedy algorithm didn't terminate one could find infinitely many subgreedy expansions by running odd greedy for some number of steps then switching to greedy. However, the conjecture abstracts away the complications of oddness, leaving a perhaps simpler problem. In addition it allows the focus to be on the number of solutions rather than the length of an individual solution, which seems a more promising line of attack.

( Evidence )

2/7 = 1/5 + 1/(5*3) + 1/(5*3*5) + 1/(5*3*5*3) + 1/(5*3*5*3*5) + ... and
3/7 = 1/3 + 1/(3*5) + 1/(3*5*3) + 1/(3*5*3*5) + 1/(3*5*3*5*3) + ...

More generally,
2/(12n+7) = 1/(6n+5) + 1/(6n+5)(4n+3) + 1/(6n+5)(4n+3)(6n+5) + ... and
3/(12n+7) = 1/(4n+3) + 1/(4n+3)(6n+5) + 1/(4n+3)(6n+5)(4n+3) + ...

So the analogue to the odd greedy Egyptian fraction problem for Engel expansion has an immediate negative answer.

The same infinite 2/7 and 3/7 expansions also occur as the Costé prime expansion of these numbers, which is infinite for the same reasons.

The sieve of Eratosthenes is normally used to generate prime numbers. But it's almost as easy and almost as efficient to make it generate the sequence of factorizations of all natural numbers: instead of eliminating all multiples of each prime p, add p to the list of factors of all its multiples. An implementation of this (with some more care about handling prime powers) is now in Eratosthenes.py in my PADS library.

As an application, I included in the same file a generator for practical numbers, based on a formula described in the link for testing whether a number is practical given its factorization. Despite being highly composite, practical numbers have a lot of similarities as a sequence to prime numbers; one of those similarities is that this formula is a lot like a sieve. But rather than implement it as a sieve, it seemed easier to use the simpler Eratosthenes sieve to generate the factorizations and then test each factorization independently.

For a prime p, define f(p) to be the minimum i such that all residues mod p can be expressed as sums of subsets of the values {1/1, 1/2, 1/3, ...1/i} mod p. We need only consider residues other than 0 and 1, since 0 is always the sum of the empty subset and 1 = 1/1. E.g.,

f(2) = 1, because 1 = 1/1.

f(3) = 2, because 2 = 1/2 mod 3.

f(5) = 3, because 2 = 1/3, 3 = 1/2, and 4 = 1/1 + 1/2 mod 5.

f(7) = 3, because 2 = 1/2 + 1/3, 3 = 1/1 + 1/2 + 1/3, 4 = 1/2, 5 = 1/3, and 6 = 1/1 + 1/3 mod 7.

f(11) = 4, because 2 = 1/2 + 1/3 + 1/4, 3 = 1/4, 4 = 1/3, 5 = 1/1 + 1/3, 6 = 1/2, 7 = 1/1 + 1/2, 8 = 1/1 + 1/3 + 1/4, 9 = 1/2 + 1/4, and 10 = 1/2 + 1/3 mod 11.

f(13) = 5. 2 = 1/2 + 1/5 (mod 13) and 12 = 1/1 + 1/2 + 1/3 + 1/5 (mod 13); all other residues can be expressed as sums of subsets of {1/1, 1/2, 1/3, 1/4}.

How quickly does f(p) grow with p? It must be at least log₂p else there would not be enough different subsets to cover all the residues. It seems plausible to me (e.g. using an argument in which we replace each unit fraction by a random variable) that f(p) = Θ(log p).

The relevance of this is for bounding the maximum denominator of an Egyptian fraction as a function of its denominator. For each p, some fraction x/p requires a maximum denominator of at least f(p)p in any Egyptian fraction expansion, as can be seen by considering only the terms in the expansion in which the denominator is a multiple of p. The simple lower bound for f(p) above thus leads to a lower bound of Ω(y log y) on the maximum denominator in an expansion of x/y, for infinitely many pairs x, y. On the other hand, it seems plausible (although not completely obvious due to issues with repeated fractions and overlarge expansions) that something like the prime-multiples method of my Egyptian fraction code can expand any fraction x/y using max denominator f(p)y, where p is the prime factor of y with the largest value of f(p). This approach might have some hope of improving the Tenenbaum-Yokota O(y log^2 y / log log y) upper bound on the max denominator.

The Egyptians expanded 2/101 = 1/101 + 1/202 + 1/303 + 1/606. That expansion is not the shortest (1/51 + 1/5151) but it has the smallest possible maximum denominator (606) of any expansion as a sum of distinct unit fractions. The expansion 2/n = 1/n + 1/2n + 1/3n + 1/6n works for all n; for how many other values of n does it produce the optimal result?

Answer: infinitely many, by the following theorem.

Theorem. Let D(q) denote the smallest d such that q has an Egyptian fraction with maximum denominator d, and let k be a fixed integer. Then for all but finitely many prime numbers p, D(k/p) = p D(k).

( Proof )

Testing all expressions of the form 2 - sum 1/a_i for a_i ≤ 5 shows that the only possible odd primes with 2/n expansions having denominators better than 6n are 2, 3, 5, 7, 11, 13, 17, 29, 31, 43, and 73. Here are some better expansions for these numbers:

( Read more... )

The Prime Game, Jeff Shallit. There are only finitely many primes that do not have another smaller prime hidden in them as a subsequence of their decimal representations. True more generally for any set instead of "primes" and any base instead of decimal, but for decimal primes Jeff can actually tell you what the minimal ones are.

An Egyptian fraction is a sum of distinct unit fractions, but it's easy to convert sums of non-distinct fractions into Egyptian fractions, by repeatedly finding pairs of equal fractions 1/k + 1/k and replacing them either by 2/k (if k is even) or 2/(k+1) + 2/k(k+1) as I described here. For example,

5/18 => 1/18 + 1/18 + 1/18 + 1/18 + 1/18 => 1/18 + 1/9 + 1/9 => 1/18 + 1/5 + 1/45.

It also works to replace 1/k + 1/k by 1/k + 1/(k+1) + 1/k(k+1), as Graham and Jewett proved:

3/7 => 1/7 + 1/7 + 1/7 => 1/7 + 1/8 + 1/8 + 1/56 + 1/56 => 1/7 + 1/8 + 1/9 + 1/72 + 1/56 + 1/57 + 1/3192.

Both of these methods produce fractions that may have significantly larger denominators than the starting one, because the replacement formulas involve quadratic polynomials. But maybe that quadratic blowup is unnecessary? What if we try a much simpler replacement: 1/k + 1/k => 1/k + 1/2k + 1/3k + 1/6k?

3/4 => 1/4 + 1/4 + 1/4 => 1/4 + 1/8 + 1/8 + 1/12 + 1/12 + 1/24 + 1/24 = 1/4 + 1/8 + 1/12 + 1/12 + 1/16 + 1/24 + 1/24 + 1/24 + 1/48

and now we have a problem: the three copies of 1/4 led to three copies of 1/24, which will lead to three copies of 1/144, which...won't ever terminate.

Can we fix this by a more complicated replacement rule? What about, say, 1/k + 1/k => 1/2k + 1/3k + 1/7k + 1/42k? No such expansion rule, in which all the denominators of the expansion are linear in the duplicated input denominator, can work! For, in a rule with s terms in the expansion, starting with a copies of the same fraction and expanding one of these copies once, one copy twice, etc., results in an expansion with ≥ s^a-1 terms, all of which belong to a set of O(a^s) different multiples of the starting denominator. So for any sufficiently large a, we can find after these forced expansions a larger denominator d with at least a copies of the fraction 1/d, setting off the same kind of infinite loop we saw for the 1/2 + 1/3 + 1/6 expansion.