The law of quadratic reciprocity is a beautiful result in elementary number theory which relates the Legendre symbols and . For and odd primes, it states that

One of the better-known proofs of quadratic reciprocity involves the quadratic Gauss sum, an exponential sum defined by

in which . In essence, this proof works by computing in two ways:

1. Euler’s criterion: . Note that this is well-defined because is an integer.
2. Fermat’s Little Theorem (for algebraic integers): Working in the cyclotomic field , we find that . Changing variables gives .

It follows that , which implies quadratic reciprocity.

In this note, we give a second proof of quadratic reciprocity using the analytic properties of the Jacobi theta function , which is a modular form of weight . In particular, we produce two formulas for the behavior of as , and derive quadratic reciprocity by comparing them. As in the previous proof, Gauss sums play a prominent role. Otherwise, this proof is distinct; in particular, it avoids Euler’s criterion and any semblance of algebraic number theory.

— THE JACOBI THETA FUNCTION —

The Jacobi theta function is defined on the upper half-plane by the series . It is a modular form of weight on the congruence subgroup , transforming via

in which for , for , and denotes Kronecker’s extension of the Legendre symbol. In addition, satisfies the functional equation , which follows from Poisson summation and analytic continuation. (In both functional equations, we use the principal square root on .)

We only need the second functional equation in our proof, but we will use it several times. For starters, it implies that as from above. And more generally, this asymptotic implies that

as , for any choice of . (Hint: Use the squeeze theorem to compare the asymptotics with different choices of .)

— THE BEGINNING OF A PROOF —

We now begin our proof in earnest. Let and be odd primes, and consider the limit of as . We have

In the limit as , each inner sum grows as . For the outer sum over , we let denote and observe that

It follows that as . We will use this asymptotic a few times in what follows.

From the functional equation , we see that

As , the right-hand side tends to . Meanwhile, the left-hand side equals . Ignoring the non-dominant term, this tends to

The sum over is a generalized Gauss sum. By the Chinese remainder theorem, we have

The first exponential sum is . The second depends only on and can be computed exactly; we see that it equals when and when . It follows that

in the limit as .

Of course, these two expressions for the limiting behavior must agree. After some simplification, we conclude that

Note that the right-hand side has absolute value because as a consequence of , and the same for . (Note that the cases in the previous line could be written .)

— THE SIGN OF THE GAUSS SUM —

Recall that the Gauss sum is one of the complex square roots of . At this point, we could prove quadratic reciprocity if we knew which square-root (depending on ) to choose. The correct choice of sign is of course well-known now, but it took Gauss a few years to determine the sign after first introducing the Gauss sums. Fortunately, we find ourselves in a position to determine the sign of immediately.

The key observation is that none of our formulas so far have required and to be odd primes. In fact, it is enough for and to be odd, coprime, and square-free. Taking prime and in the previous offset equation gives

Of course, as well. Returning to the same offset equation from before, it follows that

The right-hand side only depends on the values of and mod . And, as expected, it equals unless both and , when it equals . This completes our proof of quadratic reciprocity.

— EXERCISES —

Exercise: Prove the second supplement to quadratic reciprocity by considering as .

I’ve recently decided to learn Python, which means that I’ve been Googling a lot of simple search terms to learn command syntax. One of these searches led to an invitation to foobar, an online coding challenge run by Google. This post concerns one of the early problems in Google foobar, lovely-lucky-lambs.

I wouldn’t usually comment on a problem like this, but the accepted answer for the lovely-lucky-lambs problem is false. This doesn’t seem to be acknowledged anywhere, which is strange given the reach/duration of foobar (and Cunningham’s Law). In this post, I’ll introduce the lovely-lucky-lambs problem and discuss why the official test case answers are sometimes incorrect. This post contains major spoilers, so be warned.

Here is the lovely-lucky-lambs problem, in its entirety:

Lovely Lucky LAMBs

==================
Being a henchman isn’t all drudgery. Occasionally, when feeling generous, Commander Lambda hand out Lucky LAMBs (Lambda’s All-purpose Money Bucks). Henchmen can use Lucky LAMBs to buy things like a second pair of socks, a pillow for their bunks, or even a third daily meal!

However, actually passing out LAMBs isn’t easy. Each henchman squad has a strict seniority ranking which must be respected — or else the henchmen will revolt and you’ll all get demoted back to minions again!

There are 4 key rules which you must follow in order to avoid a revolt:

1. The most junior henchman (with the least seniority) gets exactly 1 LAMB. (There will always be at least 1 henchman on a team.)

2. A henchman will revolt if the person who ranks immediately above them gets more than double the number of LAMBs they do.

3. A henchman will revolt if the amount of LAMBs given to their next two subordinates combined is more than the number of LAMBs they get. (Note that the two most junior henchmen won’t have two subordinates, so this rule doesn’t apply to them. The 2nd most junior henchman would require at least as many LAMBs as the most junior henchman.)

4. You can always find more henchmen to pay – the Commander has plenty of employees. If there are enough LAMBs left over such that another henchman could be added as the most senior while obeying the other rules, you must always add and pay that henchman.

Note that you may not be able to hand out all the LAMBs. A single LAMB cannot be subdivided. That is, all henchmen must get a positive integer number of LAMBs.

Write a function called solution(total_lambs), where total_lambs is the integer number of LAMBs in the handout you are trying to divide. It should return an integer which represents the difference between the minimum and maximum number of henchmen who can share the LAMBs (that is, being as generous as possible to those you pay and as stingy as possible, respectively) while still obeying all of the above rules to avoid a revolt. For instance, if you had 10 LAMBs and were as generous as possible, you could only pay 3 henchmen (1, 2, and 4 LAMBs, in order of ascending seniority), whereas if you were as stingy as possible, you could pay 4 henchmen (1, 1, 2, and 3 LAMBs). Therefore, solution(10) should return 4-3 = 1.

To keep things interesting, Commander Lambda varies the sizes of the Lucky LAMB payouts. You can expect total_lambs to always be a positive integer less than 1 billion (10 ^ 9).

In summary, we have some number of LAMBs to distribute and rules which govern their distribution. We’d like to know the difference between the most number of recipients and the fewest, subject to those rules.

Rules (1) and (3) imply that the nth recipient (in order of increasing seniority) must receive at least LAMBs, where is the nth Fibonacci number (indexed starting at ). Let denote the sum of the first Fibonacci numbers. Then any distribution of LAMBs to henchmen requires at least LAMBs; in other words, if we have fewer than LAMBs, we must allocate LAMBs to fewer than henchmen.

Thus an upper bound for the number of recipients is the largest for which . This upper bound is clearly sharp. It’s also easy to compute:

def most_recipients(n):
    recipients = 0
    partial_sum = 0
    while partial_sum <= n:
        recipients += 1
        partial_sum = partial_sum + fibonacci(recipients)
    return recipients -1

Here, fibonacci(n) returns the nth Fibonacci number. I won’t include code for it so I don’t get yelled at for implementing the recursion inefficiently. This code overcounts and then corrects by 1. As an exercise, we can prove that and implement most_recipients() that way. This wouldn’t do much to speed up the code, though.

The lower bound for the number of recipients is more nuanced. Rules (1) and (2) imply that the nth henchman may receive no more than LAMBs. This implies that the total number of LAMBs awarded to henchman is at most . Conversely, any distribution of more than LAMBs requires at least recipients. This lower bound for recipients is also easy to code:

def lower_bound_recipients(n):
    import math
    return int(math.floor(math.log(n+1,2)))

This function returns the maximal for which the sum remains less than or equal to . In particular, with this distribution scheme, it is not possible to add another henchman at the top of seniority and pay them twice the rate of the previous henchman. Thus the lower bound is optimal, correct?

Yes and no. It’s true that defining

def solution(n):
    return most_recipients(n) - lower_bound_recipients(n)

passes the ten predetermined “random” test cases for the lovely-lucky-lambs challenge. (Note: I suspect that the test cases are determined once and fixed for all users, but I am not 100% certain. More on this later.)

When we get into the weeds of things, though, we see that lower_bound_recipients() need not be optimal. In fact, it breaks as early as . For , we compute a lower bound of 1 recipient, since is already too large. However, the solution which allocates one LAMB to a first henchman and pockets the rest is not a valid distribution scheme: it violates rule (4), which states that we must give the remaining LAMB to a second henchman:

4. You can always find more henchmen to pay – the Commander has plenty of employees. If there are enough LAMBs left over such that another henchman could be added as the most senior while obeying the other rules, you must always add and pay that henchman.

Thus any allocation of two LAMBs must go to exactly two henchmen, and the lower bound is incorrect.

There is an easy fix: we simply compute the maximal number of henchman who can be paid in geometric progression, and then check if the surplus LAMBs are enough to demand we hire another henchman. Here is one version of the corrected lower bound:

def fewest_recipients(n):
    recipients = 0
    partial_sum = 0
    # Overcount the geometric progression by 1
    while partial_sum <= n:
        partial_sum = partial_sum + 2 ** recipients
        recipients += 1
    # Correct for the overcount:
    recipients -= 1
    partial_sum = partial_sum - 2 ** recipients

    gap = n - partial_sum

    if recipients >= 2 and gap >= 2 ** (recipients -1) + 2 ** (recipients -2):
        return recipients + 1
    elif recipients == 1 and gap >= 2 ** (recipients -1):
        return recipients + 1
    else:
        return recipients

This disagrees with the previous lower bound for . For , they disagree on the union of the intervals .

To investigate whether or not this nuance could be responsible for my test case failures, I modified my code to probe the black box test cases. Since solution() defined above passed all the test cases, a version like

def solution(n):
    if n < M:
        return most_recipients(n) - lower_bound_recipients(n)
    else:
        return -1

would pass exactly those test cases with . If we don’t value our time, we can manually binary search up to 10^9 to determine the exact used in the test cases. For me, the test cases were

(Disclaimer: I don’t know if these test cases are universal, but I expect that they are. They appear to be chosen specifically to test edge cases. For example, 75024 is the sum of the first 23 Fibonacci numbers.)

The “corrected” code fails two test cases: and . For , we’ve already mentioned that any distribution of 2 LAMBs must go to two henchman, so the difference between the most and fewest recipients should be 0. The accepted answer is 1, which is incorrect.

The second failed test case is . Here, the accepted solution is 9 and my proposed correction is 8. The maximum number of recipients is 28 (since ). For the minimal number of recipients, we note that . This gives a lower bound of 19 recipients. However, even in this extremal case (in which each henchman is paid maximally), there remains enough leftover LAMB to pay a 20th henchman. We have 393216 remaining LAMB, which exactly equals . Thus a 20th henchman would need to be added, per rule (4). (Note that a 20th henchman is barely avoided in test case 8.)

Some concluding remarks:

1. The four rules which restrict LAMB allocation were carefully chosen to produce our upper and lower bounds. Rules (1) and (3) give a lower bound for LAMBs and thus an upper bound for recipients. Rules (1) and (2) give the upper bound for payouts. Some version of rule (4) was needed so that the lower bound for recipients wasn’t 0 (keeping all the LAMBs for oneself). The specific wording of rule (4) seems to have been misinterpreted by whoever coded the internal version of solution().
2. The test cases 917502 and 917503 seem specifically designed to test the precise interpretation of rule (4). I suspect that the test cases and internal solution were written by different people.
3. This issue with lovely-lucky-lambs is alluded to in a Stack Overflow post here. The original submission attempts to implement rule (4) correctly, but makes an error along the way (and ironically passes the test case by accident). A proposed solution passes all test cases by incorrectly implementing rule (4). The potential issue with is noted in a comment, but the veracity of the accepted answers is not discussed.

These days, Blaise de Vigenère is best remembered for the Vigenère cipher, a classical cipher which saw wide use between the sixteenth and nineteenth centuries. This is actually a misattribution, as the Vigenère cipher was invented by Giovan Battista Bellaso in 1553. Still, Vigenère was a capable cryptographer, who invented his own cipher, the autokey cipher, in 1568. The autokey cipher was never used as widely as the Vigenère cipher, perhaps because small encryption/decryption errors could render the whole message illegible.

The autokey and Vigenère ciphers differ only in their method for producing a keystream, a long string of characters (equal to the length of the plaintext) generated by a short key. In the Vigenère cipher, the keystream is formed by concatenating the key over and over again, as in lemonlemonlemon… In the autokey cipher, the keystream is formed from the plaintext itself, prepended by a short key. For example, to encrypt classicalcryptography with the key lemon in the autokey cipher, our keystream is lemonclassicalcryptog. In either case, the ciphertext is formed by “adding” corresponding letters in the plaintext and keystream, by converting a,b,…,z to 0,1,…,25 and performing addition mod 26:

The Vigenère cipher fell out of use in the second half of the nineteenth century, following successful cryptanalysis by Charles Babbage in 1854 and Friedrich Kasiski in 1863. In essence, one notes that the Vigenère cipher acts as a series of interlaced shift ciphers. If the key has length k, the ciphertext breaks into k parts, with each part transformed by a part-specific shift. Thus letter frequency analysis can determine the shifts and the corresponding key letters.

It is natural to ask whether or not frequency-based statistical analysis can be used in cryptanalysis of the autokey cipher as well. Certainly the methods which produce good results in the Vigenère cipher no longer apply, since the keystream is aperiodic. In any case, it appears that very little has been written on this topic. Wikipedia, for instance, lists statistical attacks on Vigenère but relies on cribbing to cryptanalyze the autokey cipher. The goal of this note is to show two ways in which statistical analysis can be used to analyze autokey ciphers: as a means to identify them “in the wild” and as a means to decrypt them and recover keys.

— 1. IDENTIFYING AUTOKEY CIPHERTEXT —

Our jumping-off point is the observation that the letter frequencies of English text are non-uniform. (This is the foundation of all frequency analysis.) Nevertheless, letters do follow a consistent distribution. In the autokey cipher, this implies that ciphertext follows the distribution of the sum (mod 26) of two English-distributed character sequences. The two streams are correlated, but this correlation grows quite weak once the key length exceeds a few letters, so we can ignore it in computations and assume that the streams behave as independent variables. Under this assumption, the following diagram shows the distribution of letter frequencies in an autokey ciphertext:

Just to be clear, this is a model for the letter frequencies of any sufficiently long and sufficiently standard autokey ciphertext. The distribution does not depend on the keyword, which is one of the ways in which the autokey cipher differs from the Vigenère cipher. More importantly, this letter distribution is fairly distinct, so a test could be used to test whether or not a given ciphertext is encrypted using the autokey cipher.

— 2. STATISTICAL CRYPTANALYSIS —

Now, given a ciphertext which is expected to have been encrypted using the autokey cipher, we discuss cryptanalysis. To motivate these ideas, consider the following example ciphertext:

If we expect that the autokey keyword is a common English word or phrase, we can apply a dictionary attack to save time over a brute-force keyspace attack. This requires an automated way to determine “English-ness” of the automated decryption, which we can accomplish using n-gram statistics. (I’ve found that trigrams are most useful for what follows.) A dictionary attack using 10,000 randomly generated keys of length between 3 and 10 drawn from a crossword dictionary identifies some keys as more probable than others; the top fifteen keys (as measured by trigram statistics) are listed below:

1. lysosome        6. bourgeon       11. wwii
2. beloveds        7. dreamgirls     12. lookpale
3. paragoge        8. regauged       13. altitude
4. runouton        9. brushpile      14. regorged
5. aquacade       10. damosels       15. dxii

The abundance of high-performing 8-letter words in our random sample suggests that the real keyword has length 8. The reason for this is not too hard to see: when our candidate keyword has the correct length k and at least one letter in the correct position, one out of every k letters in attempted decryption must be correct (and therefore English-like). Restricting to 8-letter keys produces a refined sample of near-keys. The fifteen highest-performing candidates are listed below:

1. hemostat        6. removals       11. lemonoil
2. lemaitre        7. lemonbar       12. letslide
3. lemonice        8. bemuzzle       13. lemurine
4. remoulds        9. prionace       14. remotely
5. demomode       10. remorses       15. demotist

The most common letters in each position form the string lemonide, which differs from lemonade in only one position. Decryption using the key lemonade recovers the plaintex:

In hindsight, we stress that our dictionary attack never found lemonade in our dictionary: the deduction was done merely by piecing together information from near-keys. In particular, we could reassemble the key by studying decryptions using randomized keys. Either way, high-performing keys should share letters with the actual key.

In our example, suppose that the key length is known to be 8. By testing random keys which begin with a and are otherwise random, we can determine a rough estimate of such keys. Replacing a with other characters allows us to rank the relative performance of different key starts. If we do this for each of the eight letters in our key, we produce the following table of rankings:

The leftmost column correctly identifies about half of the letters in the key, and the other letters can be derived from context or from partial decryption of the plaintext. We could also improve our results (at the cost of a lot of time) by looking at random keys with prescribed bigram pairs.

As with any statistical attack, our results improve quite a bit as the ciphertext increases in length. The example used here has 77 characters, which is probably on the low end of tractable given the key length. A good rule of thumb asks that the ciphertext be at least ten times as long as the key.

A positive square-free integer is said to be congruent if appears as the square-free part of the area of an integral right triangle. Thus for instance 5 is congruent, because the (9,40,41) right triangle has area 20. A nice proof-by-picture demonstrates that is congruent if and only if there exists an arithmetic progression of three (integer) squares with common difference for some :

An arithmetic progression of squares. The common difference 720 (= 5 * 144) appears geometrically as the union of the four red triangles.

If is congruent, Euclid’s classification of (primitive) Pythagorean triples gives the existence of coprime such that

Conversely, any integer expressible in this form is congruent. Taking prime and , we see that . Then divides with multiplicity , so we have found a congruent number with as a factor. In particular, congruent numbers may have arbitrarily large prime factors and there exist infinitely many congruent numbers.

In concrete terms, this establishes the existence of something like congruent numbers up to . This is laughably (and unsurprisingly) weak. Heegner sketched an argument in 1952 that the congruent numbers include each prime and thereby gave congruent numbers up to . (Apparently, the details for some of Heegner’s claims weren’t put into writing until Paul Monsky’s Mock Heegner Points and Congruent Numbers in 1990.)

The main purpose of this post is to give a third proof of the infinitude of congruent numbers. This third proof is analytic moreso than algebraic and gives only an ineffective result. Nevertheless I think it has merit, because its techniques may generalize (with a lot of blood, sweat, and tears) to prove the infinitude of congruent numbers in any reasonable congruence class.

The distribution of congruent numbers is an open problem. A heuristic from random matrix theory due to Keating conjectures the following:

Conjecture (Keating, 2005):
1. Any square-free integer congruent to 5, 6, or 7 modulo 8 is a congruent number.
2. Denote by the set of square-free integers less than or equal to that are congruent to 1, 2, or 3 modulo 8 and that are congruent numbers. Then has density 0; more precisely, there exists a strictly positive constant such that

For more discussion on this conjecture, see also Henri Cohen’s Number Theory (Volume I: Tools and Diophantine Equations), pg. 453.

— PART II: A DIRICHLET SERIES —

The bijection between Pythagorean triples and three-term arithmetic progressions of squares implies that

⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀

in which is the square-indicator function and is the set of hypotenuses of primitive integral right triangles with square-free part (sfp) of the area equal to . (In theory, this is counted with multiplicity; in practice, we’ve never found two primitive triangles with equal hypotenuse which determine the same congruent number.) It’s worth remarking that the series on is absolutely convergent in . This follows from dyadic summation and the bound

which itself follows from a bound one of my recent papers, A Shifted Sum in the Congruent Number Problem. (Here, is the rank of the elliptic curve .) The Dirichlet series is thus in . I’m not sure if it is possible to make the explicit in .

— PART III: CUT-OFF INTEGRALS —

To prove that there exist infinitely many congruent numbers, we apply an integral transform to the Dirichlet series from the previous section. Our function of choice is a smooth cutoff of compact support, with Mellin transform .

Proposition: We have

Proof: We shift contours in a Mellin transform to compute

Since is , the integral over is . On the other hand, the integral transform of splits into the terms

in which is the number of representation of as a sum of two integer squares. Pointwise estimates for bound this error term by . We produce our final answer by balancing errors.

— PART IV: THE FINAL PROOF —

We now begin our proof by contradiction in earnest. Suppose that there exist finitely many square-free which are congruent. Summing over these choices of has the effect of removing the at left in the above Proposition. In particular, we obtain

Note that the right-hand side is of size (by finiteness and the existence of at least one congruent number) and that we have removed the -dependence from by finiteness. To further understand the left-hand side, we introduce another Proposition.

Proposition: When is square, we have

It follows from the equation line preceding this Proposition that

We now complete our proof by proving that the associated Dirichlet series has a double pole at . Indeed, we may write

The Dirichlet series has a double pole, which contradicts the upper bound written above. This completes our proof.

A Lehmer Sieve made of LEGO.  The pins represent quadratic residues.

In his excellent and excellently-named expository piece A Tale of Two Sieves, Carl Pomerance tells the history of two factoring algorithms, the quadratic and number field sieves. This story begins with a factorization technique called factorization via congruence of squares, which attempts to factor by constructing a non-trivial congruence of the form .

If such a congruence is known, then , and the factors of distribute between and . If we’re lucky (or at least, not unlucky), each of and gets at least one prime factor of and we can detect such a factor by computing .

This is a great start but it doesn’t tell us how to produce these congruences of squares in the first place. Kraitchik’s idea back in the 1920s was to find them by multiplying together non-square quadratic residues with “simple” prime factorizations. We illustrate this by (stealing Pomerance’s) example:

Example 1. Let’s factor the integer . If we wanted to apply Fermat’s Factorization Method (factorization using differences of squares), we’d consider the values of , starting at , and look for squares. We produce the list

These are quadratic residues (by construction), but none of these are squares, so Fermat would keep on looking. But what Kraitchik observes is that some of these outputs factor entirely into small primes:

By looking at prime factorizations, it’s not hard to find a combination of these outputs that gives a perfect square. For example, . (There are a few others, too.) But, let’s not forget that these numbers came from values of . In this example, we used the -values and . Thus

Since , this implies . Then , and by computing

we get a non-trivial factor of . The full factorization is . //

The purpose of this article isn’t to rehash A Tale of Two Sieves. Rather, I’d like to explore the bits and pieces that make Pomerance’s example work so nicely. For contrast, here’s a bad example:

Example 2. Let’s factor . We define and compute , , and . We thus produce a prime, a number with the ungainly factor , and a perfect square. In this case, Kraitchik’s Method reduces to the method of Fermat.  It would have been easier to compute .  //

But that’s not the only way things can go wrong. Here’s a failure of a different sort:

Example 3. Let’s factor . We define and compute values of starting at . In line with Pomerance’s example, we want to focus on -values for which admits no prime factors larger than a given bound ; such integers are called -smooth. To stay specific, let’s choose here. The smallest value we test for smoothness is ; after this, we compute , , and so forth. Each has a prime factor larger than , so we choose to discard them all.

Around now, the numbers start getting large. For the sake of efficiency, we might stop factoring completely. Instead, since we’re looking for -smooth outputs, we can perform trial division by small primes and give up if seems to misbehave. But herein lies the problem: if we cut off trial division at , say, we’d miss out on and keep on going long after Fermat’s Factorization Method terminates.

This spirals. Searching even two million values turns up none which are -smooth. Are there tricks we could use to salvage this problem? Of course. We could increase the smoothness bound. We could switch to looking at when the values of become too large. Still, these are things we’d hate to have to do when presenting an early example. //

Given all this, just how special is Pomerance’s example ? Is it misleadingly (or in some sense even optimally) nice, or is it just one of many reasonable to consider? To answer this, let’s fix some properties that a “nice” introductory example should have:

1. should have no factors less than . I’ll be the first to admit that the bound is somewhat arbitrary, but passes this test (barely) and we certainly want to exclude integers which fall too easily to trial division. (At the bare minimum, things get strange if is even.)
2. should be a semiprime. This goes hand in hand with (1), since any number with three distinct prime factors of size at least already exceeds . Moreover, it seems natural to look at semiprimes because one of the main “industrial applications” of factoring is the breaking of RSA keys.
3. The example shouldn’t devolve into Fermat’s Method. Factorization based on congruences of squares is often pitched as an improvement to Fermat’s Method. It’s therefore silly to work through an example that reverts back to this simpler technique. (See Example 2.)
4. The example shouldn’t require an overly large smoothness bound. Here, large is a bit subjective. Pomerance notes that the optimal smoothness bound for large is

   

   Being charitable, we’ll call an example “bad” if the values of appear to escape to infinity before a square-producing subset of -smooth numbers is found, where is as above.

These conditions create a fairly sparse of potential . Pomerance’s example remains in the list, as do smaller examples:

Example 4. Let’s factor using Kraitchik’s Method. We form and compute and . Then . The computation

gives a non-trivial factor of . //

I suppose the funny thing about this example is that it has succeeded a bit too well. To better highlight Kraitchik’s Method, let’s also restrict to the cases where first two outputs of don’t multiply to a perfect square. This eliminates and ; overall, only a slight change.

Example 5. Let’s factor using Kraitchik’s Method. As before, we form and compute:

Then , and a GCD computation gives the non-trivial factor of . //

This is the first of our examples that feels like a solid introduction to Kraitchik’s Method. We have some smooth outputs, some not-so-smooth outputs, and enough terms to get you thinking about how you might find a square-producing subset in a more general case.

Pomerance does a few other computations with which may explain his particular choice. It demonstrates, for example, that the first subset of quadratic residues which produces a perfect square might require more than two terms. The case requires four terms: , , , and , although relations with fewer terms appear if one continues the computation.

If one wants to arrange for the first square-producing subset to use more than two terms, the list of admissible shrinks further. The full list of the such is given below:

Example 6. Let’s factor using Kraitchik’s Method. As always, we form and compute some values of . In this case, quite a few are insufficiently smooth. The useful ones are as follows:

In this case, the first congruence of squares we find is not useful: it reads . We’d need to keep going. Unfortunately, the next smooth value we find is , so this example (secretly) devolves into Fermat’s Method.

Example 7. Let’s factor using Kraitchik’s Method. As always, we form and compute some values of . Omitting the insufficiently smooth ones, we consider

Then , which is enough to find the non-trivial factor of . At long last, we’ve produced a truly beautiful example. //

In the previous example, we defined “sufficiently smooth” as -smooth. We are, after all, being five times more lenient in our definition of smooth than what is asymptotically optimal. Tightening this bound from to disqualifies and makes the new minimal example:

Example 8. Let’s factor using Kraitchik’s Method. As always, we form and compute some values of . Omitting the insufficiently smooth ones, we consider

Then and we produce the non-trivial factor of . This examples matches in complexity but requires less trial division to perform. //

As a parting gift, here’s a full list of which are admissible under our many conditions and the stricter smoothness bound :

This list contains Pomerance’s example and others. I imagine that most, if not all, would be fine “first examples” in the study of Kraitchik’s Method. (Although I have not verified that they avoid the pitfalls of Example 6.)

A magnitude-argument plot of the gamma function.

Stirling’s approximation is a useful approximation for large factorials which states that the th factorial is well-approximated by the formula

as . Stirling’s approximation was first proven within correspondence between Abraham de Moivre and James Stirling in the 1720s; de Moivre derived everything but the leading constant, which Stirling eventually supplied (without proof; it’s not known how Stirling guessed it).

Given the rapid growth of the factorial, it is often more convenient to take logarithms and work with those instead. Doing so, Stirling’s approximation assumes the form

Most proofs of Stirling’s approximation work with the formulation above. One typical (modern-ish) proof begins by studying the integral approximation

This version (as written) is too imprecise to recover Stirling’s Approximation, but the error term can be improved massively by introducing more terms using the Euler–Maclaurin formula (ca. 1735). This is the proof sketched on Wikipedia, for example.

The precise statement of the Euler–Maclaurin formula makes reference to the Bernoulli numbers . These show up all over combinatorics, but as a multiplicative number theorist I will always think about Bernoulli numbers in the context of the Riemann zeta function. Specifically, one has , so the Bernoulli numbers relate to special values of the zeta function.

This point of view inspired me to derive Stirling’s approximation (and the additional terms making up Stirling’s series) in a way which makes the role of the zeta function obvious. For convenience, we’ll phrase everything in terms of the gamma function; this affects the shape of our formula in a small and readily-understandable way. Without further ado, here’s the proof:

Proof: We begin with Weierstrass’ infinite product for the gamma function (ca. 1854),

in which is the Euler–Mascheroni constant. (Euler gave an infinite product for the gamma function back in 1729 which might also work for this proof technique.) Taking a principal branch of the logarithm, we produce

Our next ingredient is a contour integral representation of the logarithm; namely,

in which the line of integration is the vertical line with real part . This formula can be proven by Cauchy’s residue formula (and analytic continuation in ). We note as well an interpretation via Mellin transforms: this formula states that the Mellin transform of is . (See here for a table listing this and other Mellin pairs.)

In any case, by combining our formulas and shifting the contour to the right, we produce

Absolute convergence of the contour integral justifies an interchange of sum and integral and allows us to recognize a zeta function:

We can now shift the contour of integration back to the left, extracting residues as we go. The residue from the double pole at is , hence

(From here on out, we identify the Riemann zeta function with its meromorphic continuation.) Shifting farther left passes a second double pole at , which has residue , hence

A bound in absolute value shows that the integral at right in the line above is and therefore small enough to recover Stirling’s approximation.

With the proof complete, I offer two final remarks:

1. There’s nothing stopping us from shifting the line of integration in our formula even farther left. Doing so passes simple poles at the negative odd integers, which we can extract as residue terms in our formula. (The poles at negative even integers are cancelled by trivial zeros of the zeta function.) This yields

which is equivalent to Stirling’s series.

2. In one of Keith Conrad’s online notes, he subdivides proofs of Stirling’s approximation according to the origin of the mysterious factor . Conrad lists references in which is obtained through (a) infinite products (such as the Wallis product), (b) the Gaussian integral, and (c) reflection formulas for the gamma function. In the proof we give here, we obtain from the series expansion

in the limit as . Because the functional equation of zeta invokes gamma functions, you could argue that this proofs falls in camp (c). However, the proof here is quite different from that given in the Conrad’s example for (c) (which was taken from section 5.2.5 in Ahlfors’ book Complex Analysis). Both proofs begin with the Weierstrass product, but their use of contour integration is quite different. The proof given here is also shorter.

Members of the 4-free set G_4, in red.

— 1. The Erdős–Turán Conjecture —

In combinatorics, a set of positive integers is called large if the sum of the reciprocals of the elements of diverges. In other words, large sets are exactly those with divergent harmonic sums. Some examples of large sets include:

1. The set of positive integers in a given congruence class, by comparison to the divergent harmonic series..
2. Any set with positive natural density, by definition of density.
3. The set of primes, by a result of Euler.

The Erdős–Turán Conjecture claims that any large set of positive integers must contain arithmetic progressions of arbitrary long (finite) length. This conjecture is known to hold for the three sets above; trivially for (a), non-trivially for (b) by Szemeredi’s Theorem (1975), and non-trivially for (c) by the Green–Tao Theorem (2004).

— 2. Constructing Large Harmonic Sums —

A set of integers is called k-free if it contains no arithmetic progressions of length k. Thus the Erdős–Turán Conjecture may also be phrased in the following way: the harmonic sum of any k-free set is finite. Moreover, if the Erdős–Turán Conjecture holds, a remark of Gelver (1977) implies that these harmonic sums can be bounded, uniformly, by a function of k.

It is a natural question to ask what these maximal harmonic sums look like. Nothing concrete is known about upper bounds (which is unsurprising, since their existence implies Erdős–Turán), but lower bounds may be found through experimentation. For example, let be the k-free set constructed using a greedy algorithm. (That is, construct inductively by repeatedly adding the smallest integer which maintains k-freedness.) This gives:

These sequences (and links to many more related ones) can be found on the OEIS at the following links: A003278, A005837, and A020655.

Note that no element of is divisible by five. In fact, it’s not hard to show that is exactly the set of integers for which the base-5 representation of omits the digit 4. Likewise, is the set of integers for which the base-3 representation of omits the digit 2. This pattern falls apart for ; one can show that this fractal structure develops in exactly when k is prime.

Here, having structure is a good thing, because structure in limits the ability of to form arithmetic progressions “by accident.” For this reason, and are of slow-growth and therefore have competitively-large harmonic sums (among 3-free and 5-free sets, respectively). The harmonic sum of is 3.00794 and the harmonic sum of is 7.8723.

I should mention that convergence to these harmonic sums is actually quite slow. While it is not hard to show that and are small sets (i.e. have convergent harmonic sums), the problem of effectively estimating their harmonic sums remained open until the work of Baillie–Schmelzer (2008). And, for sets like without fractal self-symmetry, no algorithm is known to efficiently estimate the harmonic sum. Brute-force computation using the first 10,000 terms of shows that this harmonic sum exceeds 4.191, but does not give an error envelope.

All this leads me to suspect that we can do better than . By finding a 4-free set with more structure, we can hope to (1) avoid arithmetic progressions in a more systematic (and efficient?) way and (2) better estimate the harmonic sum of what we produce.

— 3. Fractal Sets Without Long Arithmetic Progressions —

To generalize the nice properties of and , we study arithmetic progressions in the residue field , with p a prime. Formally, call a set of residues an arithmetic progression modulo p if the have some representatives in in an arithmetic progression (reordered, if necessary). For example, is arithmetic progression mod 7: up to reordering, it’s equivalent to . By analogy, a set is called k-free modulo p if it contains no k-term arithmetic progressions modulo p.

The following theorem gives a reason to care about k-free sets modulo p:

Theorem: Let p be prime. If is k-free modulo p and , then the set of integers whose base-p digits are restricted to forms a k-free set.

Proof: Suppose that is k-free modulo p and that the set of integers with restricted base-p digits contains the arithmetic progression . If , then has k distinct residues in modular arithmetic progression, a contradiction. Otherwise, let be the base-p units digit of . The progression has digits restricted to and a smaller common difference. We conclude by infinite descent.

This theorem gives a quick proof of the fact that is p-free when p is prime, since the base-p digit set is clearly p-free modulo p. (Technically, is obtained by then shifting all elements in this digit-based set up by 1, but this has no effect on the existence of arithmetic progressions.)

This theorem also provides a new avenue for constructing k-free sets of fractal type. From any digit set which is k-free mod p, we form an auxiliary set based on base-p digit restriction which is k-free. Checking k-freedness mod p is a finite problem and can be done algorithmically (and efficiently), which gives us a reasonable means for producing many examples of k-free sets. Among many examples, we hope to find a few with large harmonic sums.

— 4. Numerical Results —

In our quest to best and improve the numerical results for 4-free sets, we begin with a (weakly-pruned) brute-force search for 4-free sets modulo p within the subsets of . This is of course computationally taxing (the run-time is exponential in p). In practice, I’m limited to primes , though more effective pruning (or completely novel techniques) could push our search a lot farther.

Fortunately, we don’t have to go very far at all before we find some sets which produce larger harmonic sums than . The simplest is the set of integers for which the base-11 digits of restrict to the set . This set has a harmonic sum of 4.4217, which seems to beat that of by a considerable margin. (You might argue that this sum only beats because Baillie–Schmelzer gives better sum estimation. It’s possible that outperforms, although this seems unlikely since a comparison of the first 10,000 terms in each set shows that is already behind in harmonic sum.)

In the course of my search, I found approximately thirty other sets with harmonic sums larger than 4, of which three had harmonic sums over 4.1. These three employ the 4-free digit sets

and produce the harmonic sums 4.1968, 4.1477, and 4.1267, respectively.

The lead of this base-11 set over all others in my data-set is quite striking. I doubt that that set has a maximal harmonic sum (among sets of this type), but good solutions seem sparse and it appears that new techniques are needed to advance this search.

That’s it for today. I’m open to ideas advancing this problem, both computationally and theoretically. The application of these generalized fractal sets in the Erdős–Turán Conjecture seems new and I wonder if their properties will be helpful in that context.

Substitution ciphertext from Edgar Allan Poe’s “The Gold Bug.”

A (monoalphabetic) substitution cipher is an encryption scheme in which characters in the plaintext and ciphertext correspond in a one-to-one fashion. Substitution ciphers include, as special cases, well-known and classic ciphers like ATBASH and the Caesar shift. More generally, the underlying permutation (known as the key) will have no special structure: it will be a random permutation drawn from the space of 26! potential keys. The size of this keyspace is fairly large, even by modern-ish (ca. 1990s) standards, but a large keyspace does not in itself guarantee security.

It was known as early as ~850 CE (in the work of al-Kindi) that frequency analysis could be used to crack substitution ciphers. Despite this, substitution ciphers (in various and increasingly-complex incarnations) saw continued use until the end of the late medieval period (ca. 1500 CE), when they were supplanted by more robust creations, like the Vigenère cipher.

In general, the number of characters required to carry out a brute-force keyspace attack on a cipher is called the unicity distance of the cipher. The substitution cipher on an English lowercase alphabet has a keyspace of size . All keys being equally likely (which may or may not be true in practice!), we’d expect to need 88.4 bits of information to decrypt a ciphertext encrypted in this way. Theoretically, an English letter could represent bits; in practice, redundancies limit this to around 3.2 bits of information. Thus one expects to need at least 28 characters to decrypt a substitution cipher of this particular form.

More realistically, about 50 characters are required to crack a substitution cipher, although this number increases if the text is unusual (for example, if the letter distribution is unusually flat or if the writing is otherwise constrained). The “e”-less phrase,

of 57 characters, is one such example. Given a “normal” string of text (such as an un-parsed lower-case character string from a random article on Wikipedia), 50 characters is a pretty good stretch goal for computer-assisted solves.

In what I’ve seen, there are two main algorithms employed for the automated solution of substitution ciphers. The first is combinatorial in nature: word lengths are guessed, then letters are assigned to pattern-match words in a large dictionary. (This technique is particularly efficient in the variant problem where spacing is preserved. See one implementation here.) The second algorithm is more probabilistic: a statistical model is imposed on the k-grams of the text, and letters are swapped so as to maximize the “fitness” of the decrypted text. (See one example here.)

There is plenty of variation in each approach. In the second, for example, one can select letter swaps based on hill-climbing or genetic algorithms. In my experience, genetic algorithms are more useful because the fitness landscape in this problem is rife with local maxima.

I’ve linked a Mathematica file (and some auxiliary data files) at the end of this post which implements a genetic algorithm to solve un-parsed substitution ciphers. The algorithm makes an initial guess at the key based on letter frequencies (i.e. 1-grams), then mutates and appraises candidate keys based on progressively longer k-grams. It can decrypt most 70-character strings in fewer than 100 generations. For shorter strings, the main challenge seems to be in the choice of fitness functions: it’s hard to produce fitness functions for which the desired plaintext is a global (or even local) maximum.

Click here to access a Dropbox folder containing the aforementioned files.

At this moment, the fastest general-purpose factoring algorithm for extremely large integers is the general number field sieve (GNFS), invented by Pomerance in 1989 (as a generalization of Pollard’s special number field sieve from 1988). The GNFS has a heuristic runtime which is sub-exponential in the size (in bits) of the input. In particular, the GNFS is much, much faster than trial division at finding factors of large integers. (The worst-case runtime for trial division is exponential in the bit size.)

This is not to say that the GNFS is the exclusive tool used in finding factorizations today. The algorithmic overhead required to implement the GNFS is immense, so much so that it only starts to eclipse other algorithms when used to factor integers of 400 bits or more. (This translates to around 120 decimal digits.) For smaller inputs, different methods will out-perform the GNFS:

• Trial Division / Look-up Tables
•   SQUFOF / Pollard’s Method
•   Elliptic Curve Method (ECM)
•   Quadratic Sieve
•  General Number Field Sieve (GNFS)

For one concrete example, the computer algebra system Mathematica uses a combination of trial division, Pollard’s Method, Pollard’s Method, ECM, and the quadratic sieve to implement its built-in function FactorInteger[]. Note that Mathematica doesn’t even bother to implement the GNFS. In practice, the GNFS is only used in dedicated, distributed computing projects.

The quadratic sieve was the major inspiration for the GNFS and operates by searching for congruences of the form . Pollard’s Method is a sort of Monte Carlo method, which obtains factorizations from cycles in the iterates of a pseudo-random map.

The other three algorithms Mathematica uses (trial division, Pollard’s Method, and the elliptic curve method) are all examples of algebraic group law factorization algorithms. Broadly speaking, these algorithms perform computations in an algebraic group defined mod for which the group structure decomposes (via Chinese Remainder Theorem) into components corresponding to the unknown prime powers dividing . If we can find an element of which projects to the identity in at least one but not all of the components of (and can detect this!) we can factor .

This note discusses some of the more successful algebraic group law factorization algorithms, which include trial division, Pollard’s Method, Williams’ Method, and the ECM. We start simple, and work our way up in complexity.

— 1. TRIAL DIVISION —

To find a non-trivial factor of using trial division, we simply check if for each . Unless is prime (which is detectable in polynomial time), we’ll find a factor before . Now, this isn’t what I’d call an algebraic group law factorization algorithm (AGLFA). To appreciate trial division as an AGLFA, we need to re-interpret some of our steps.

The underlying group in this case is the integers under addition. Viewed mod , we are looking at the additive group of . For simplicity, let’s assume that has prime factorization (as in the case of an RSA modulus). By the Chinese Remainder Theorem, we can write

Recall, as a general goal, that we’re looking for an element which equals the identity in exactly one of the components or . One way to do this is to fix an element (at random, say) and look at the iterates of under the group action. If we can find an integer such that

we will have found our partial identity. But how can we find without knowing or ? The idea is to search for in the sequence . To identify a successful choice of , we need a way to determine if is a partial identity. We can do this with a GCD computation, since if and only if is at least a partial identity. So, without further ado, here’s “trial division” as an AGLFA:

1. Choose a convenient value for .
2. Starting at , replace with .
3. Compute . If , return , a non-trivial factor of (and stop).
4. Increment and return to Step 2.

It’s easy to look at the algorithm above and see the ways in which it might be slower than traditional trial division. The trial division you’re used to takes mod operations, in which is the smallest prime factor of . This takes time when and are both large. On the other hand, this algorithm involves modular multiplications of size and GCD computations. The runtime is therefore . The runtimes are thus comparable when is large, although the implicit constants certainly favor the traditional version.

Remark — I’m aware just how silly this algorithm looks. The point here is just to create an analogy which we’ll develop in the next few sections. Some of the core ideas, like taking or finding factors by computing GCD against , will be seen time and time again.

— 2. POLLARD’S METHOD —

Pollard’s Method originates in a paper by Pollard from 1974 and is typically viewed as the archetypal AGLFA. In this case, the underlying group is , the multiplicative group of units in . Since the Chinese Remainder Theorem is a ring homomorphism, we have .

As in the additive case, we’ll search for a factor by taking a random and looking at the iterates of under the group action. If we can find an integer such that

we can recover as . As before, we search for in the sequence . This leads to Pollard’s Method, which — when written in pseudo-code — looks almost exactly like our version of trial division:

1. Choose a convenient value for .
2. Starting at , replace with .
3. Compute .
   • If , return , a non-trivial factor of (and stop).
   • If , return to Step 1 and choose a new .
4. Increment and return to Step 2.

Just how fast is Pollard’s Method? We note that if and only if , and that by Lagrange’s Theorem. In a worst-case (but not atypical) scenario, we’d have to increase until . As a rough heuristic, Pollard’s Method frequently terminates when is the largest prime factor of .

The runtime of Pollard’s algorithm depends dramatically on the shape of the hidden factors and . For example, Pollard’s algorithm factors

in just twelve iterations, while the seemingly-related

takes over a million iterations to complete. When Pollard’s method is at it’s worst — for example when and are both safe primes — the algorithm is no faster than trial division. For this reason, Pollard’s Method is typically classified as a special-purpose factoring algorithm; it won’t be efficient for most inputs, but there are cases when it is the only tractable approach. This is the role it plays in Mathematica.

The runtime for Method is typically , in which is the largest prime factor of (or , if is found first). This can be contrasted with our algorithm based on trial division, which ran in time , in which is the largest prime factor of (i.e. ). Obviously, will be smoother than , but their smoothnesses might be of equal magnitude.

The difference between (in Pollard’s Method) and (in trial division) can be explained in terms of Lagrange’s Theorem. Moving from to , we move from a group of order to one of order . Since a lot of the variety in AGLFAs comes from choosing groups with different orders, expect more of this in the following sections.

— 3. WILLIAMS’  METHOD —

Williams’ Method was first presented in a paper by Williams in 1982. This special-purpose factoring algorithm is directly inspired by Pollard’s Method (the prototypical AGLFA) and, as the name suggests, works best when has a prime factor for which is smooth.

The underlying group in Williams’ Method is the multiplicative group of the quadratic ring extension

in which is a quadratic non-residue mod . If , the Chinese Remainder Theorem implies that , which also extends to their unit groups.

When is prime, every element of has a multiplicative inverse. To be specific,

in which the denominators represent multiplication by . (If this last inverse fails to exist, , which implies since was chosen to be a quadratic non-residue.) Note that this implies that by Lagrange’s Theorem.

If we adopted Pollard’s Method wholesale, we’d now attempt to find a unit and an integer for which

Unfortunately, it’s not so easy, and we’re not quite done. Since the multiplicative group of a field is cyclic, it’s certainly possible (and not at all unlikely, even) that . In this case, searching for a candidate in the sequence would require to be at least as large as the largest prime factor of . Since , we’re impeded by large factors of and of , which makes our current approach no better than Pollard’s Method.

In theory, fixing this is simple: we just consider instead. But, lest we forget, we don’t know the value of ahead of time. And yet, somehow, this doesn’t stop us. To see why, consider

in which we’ve used the binomial theorem, the fact that unless , Euler’s Criterion, Fermat’s Little Theorem, and the fact that is a quadratic non-residue. It follows that

This formula still refers to , but we can make the visible -dependence go away by computing the right-hand side as an element . Now, we attempt to find an for which

.

By construction, is a divisor of . Moreover, we expect once gets as large as the largest prime factor of . But once this occurs, we can detect it by computing a certain GCD. This comes together as Williams’ algorithm:

1. Choose some which is a quadratic non-residue mod all primes .
2. Choose a random/convenient value and compute
   
3. Starting at , replace with .
4. Compute or .
   •  If , return , a non-trivial factor of (and stop).
   • If , return to Step 1 and choose a new and .
5. Increment and return to Step 3.

As you might guess, the typical runtime of Williams’ Method is , in which is the largest prime factor of . There are plenty of cases, for example

in which the Williams’ Method outperforms both trial division and Pollard’s Method. Still, we mustn’t forget that Williams’ Method is a special-purpose factoring algorithm, which is typically no faster than trial division.

Remark — There is no known algorithm to deterministically perform Step 1 without essentially factoring along the way. In practice, one chooses at random and hopes that it was a quadratic non-residue mod the prime for which is smoothest. This will happen about 50% of the time, so it’s a good idea to run Williams’ algorithm more than once if it fails the first time. It’s also worth noting that Williams’ Method degenerates into a less efficient version of Pollard’s Method when is a quadratic residue, which is not ideal but also not the worst outcome.

There’s nothing stopping a program like Mathematica from incorporating Williams’ Method into their factoring suite. That being said, I’ve never seen it put to the test in a serious context. It may be that the Williams’ Method fails too often to bother implementing.  (See here for further comparison to Pollard’s Method.)

— 4. FACTORING WITH CYCLOTOMIC POLYNOMIALS —

The which appears in Williams’ Method arises from the computation . If we replace the quadratic extension with an extension of degree , we’d produce a factoring method which relied on the smoothness of . With a bit of care, we can eliminate some of the factors of and require only that be smooth, where is the -th cyclotomic polynomial. This completely subsumes the algorithms due to Pollard and Williams, since and .

This technique was first published in 1989 by Bach and Shallit. Unfortunately, it remains mostly a curiosity, since is far less likely to be smooth as increases in size. The degree of is , which is at least two for . Since the likelihood that is smooth decreases as , the odds of success with a larger cyclotomic polynomial are negligible. This is a common problem with AGLFAs, one which really restricts the viability of algorithms in this class.

— 5. LENSTRA’S ELLIPTIC CURVE METHOD —

The elliptic curve method (ECM) was first presented in a paper by Lenstra in 1987. Like Williams’ Method, the ECM was directly inspired by Pollard’s Method. Here, the group of integers is replaced by the group of points on an elliptic curve.

More specifically, one picks a random elliptic curve over and a non-identity point on . Using the addition law on elliptic curves, we compute for in the sequence . By the Chinese Remainder Theorem, this addition can be viewed mod , for each prime . If ever , the -coordinate of will be divisible by . But this can be detected with a GCD computation, which gives us a method for factoring .

If we restrict to a single curve , this algorithm functions exactly like Pollard’s Method, except with replaced by the order of the group . By Hasse’s Theorem, , in which depends on and satisfies . If we’re lucky and is particularly smooth, we’ll factor in short order.

More likely, will not cooperate. But here lies the advantage of the ECM over Pollard’s Method: we can throw out and repeat our process with a new elliptic curve. Under the Sato–Tate Conjecture, varies in the interval according to a well-defined distribution. In particular, by trying many elliptic curves, we expect to find one for which is sufficiently smooth in a predictable amount of time.

By optimizing our definition of sufficiently smooth (and applying some heuristics for the number of smooth numbers in an interval), we obtain an algorithm which is sub-exponential in the smallest factor of . For this reason, the ECM is most commonly used to identify factors of moderate size (i.e. fifty-ish digits) in enormous composites. The largest factor found to date with ECM is an 83-digit divisor of , found by Propper in 2013. (The leftovers form a 202-digit composite.)

The generalization from elliptic curves to curves of higher genus is easy to do, but in the end, seems to be unhelpful. The reason for this is the same reason that factoring with cyclotomic polynomials is impractical: the odds of finding group with sufficiently smooth order decreases dramatically as the order increases in size.

Image Credit: Angela Yuriko Smith

Around the fifteenth century, technological advances such as the magnetic compass, lateen sail, and stern-post rudder enabled Europe to enter what is now known as the Age of Discovery. Oceanic routes, at first restricted to coastal areas, spread into open ocean with the adoption and subsequent improvements to the astrolabe and quadrant.

Navigation became increasingly computational without visual clues from the coastline. In the sixteenth century, navigators relied on ephemerides (journals of positions of astronomical bodies at known times) to determine position and heading. It’s no stretch to say that navigators became astronomers, who used spherical trigonometry to create models of the earth and the sky above.

The fundamental identity in spherical trigonometry is the spherical law of cosines, which – like its Euclidean cousin – relates the angles and side lengths inside a (spherical) triangle. With reference to the figure below, this formula may be written

If all but one quantity is known in the formula above, the last can be obtained through multiplication, division, and application of a trigonometric table. In this process, hand multiplication became an algorithmic bottleneck, since it was infeasible to create large multiplication tables. (Unlike trigonometric functions, multiplication is a binary operation.)

Today, we know that the algorithmic problem of multiplication can be solved using the logarithm (and tables of its values), since the identity

relates a product (inside a logarithm) to a sum (of logarithms). To multiply , one simply finds and , adds them, and looks in the table to find an for which . The process of constructing logarithmic tables was still time-consuming, but it was a process that only had to be done once (and could be done on dry land). Once made, application of these tables was easy. In the (perhaps apocryphal) words of Pierre-Simon Laplace, logarithms were

…an admirable artifice which, by reducing to a few days the labor of many months, doubles the life of the astronomer, and spares him the errors and disgust inseparable from long calculations.

While logarithms are, in hindsight, the simplest way to perform multiplication with a single variable look-up table, they were not the first. That honor goes to the Babylonian Quarter Square Algorithm, which applies the identity

to reduce the multiplication of to two squaring operations. The algorithm could be made efficient through the use of a table of squares, but it appears that this was not implemented until Antoine Voisin’s table of 1000 integral squares in 1817.

There was another technique, invented in the late sixteenth century, that saw widespread use in navigation for about twenty-five years before the introduction of Napier’s logarithm in 1614 and the first logarithmic table in 1617. This technique was called prosthaphaeresis, after the Greek prosthesis (addition) and aphaeresis (subtraction), and used trigonometric angle addition formulas to reduce multiplication to addition/subtraction with the assistance of trigonometric tables. For example, one has

To use prosthaphaeresis to multiply , one scales down and to lie in , then represents each as a value of sine (through application of a sine table) to find and . The values are added (prosthesis) and subtracted (aphaeresis), then and are computed using another table look-up. Their average gives , which is rescaled to give .

The angle addition laws, as relayed within European mathematicians, are believed to originate with Jost Bürgi (1552-1632). Globally, they predate Bürgi by several centuries, having been derived by the Persian mathematician and astronomer Abū al-Wafā Būzhjānī sometime in the late tenth century.

As a historical anecdote, the method of prosthaphaeresis highlights the role that technological advancement (in particular, the newfound ability to traverse open ocean) plays in the development of mathematics. Prosthaphaeresis also highlights the way in which trigonometric functions paved the way for other transcendental functions like the logarithm (1614), exponential function (1748, although the idea of exponentiation predates this), and the solutions to differential equations.