Category Archives: Mathematics

Small Squares

There are many algorithms for finding the prime factors of a large integer n, but the most effective are based on the simple idea of finding integers x and y satisfying x2y2 (mod n). Then (x + y)(xy) ≡ 0 (mod n) and there is a good chance that gcd(x + y, n) will be a non-trivial factor of n. Different factorization algorithms use different methods to find x and y that have the same square modulo n, ranging from Fermat’s method which can be worse than trial division to the number field sieve which is the best known factoring algorithm. A simple method is to randomly choose an integer x in the range √n < x < n and check whether x2 mod n is a square integer y2.

Example: This article will use as a running example the 32-bit integer n = 2,768,439,589 which is the product of the two primes p = 57,641 and q = 48,029. Suppose x is randomly chosen to be 267,857,730. Then x2 mod n = 9, and therefore x2 ≡ 32 (mod n) and a non-trivial factor is recovered by gcd(x + 3, n) = q. □

The simple method can be greatly improved by using a sieve, which relaxes the requirement that x2 mod n is a square, and instead requires it to be a product of elements of a factor base FB = {-1} ∪ { p | p is a prime and pB }. Integers satisfying this condition are called B-smooth, and after discovering |FB| + 1 B-smooth integers of the form xi2 mod n linear algebra can be used to find a subset I ⊆ {1, …, |FB| + 1} such that ΠiI (xi2 mod n) is a square integer y2. Then (ΠiI xi)2y2 (mod n) as desired.

Example: Let F11 = {-1,2,3,5,7,11} and suppose the following 11-smooth integers have been discovered:

  • If x1 = 32,051,026 then x12 mod n = -20 = (-1)122305170110
  • If x3 = 135,988,713 then x32 mod n = -210 = (-1)121315171110
  • If x7 = 373,858,046 then x72 mod n = 42 = (-1)021315071110

Now modulo n

  • (x1x3x7)2
  • ≡  (x12 mod n)(x32 mod n)(x72 mod n)
  • ≡  ((-1)122305170110)((-1)121315171110)((-1)021315071110)
  • ≡  ((-1)224325272110)
  • ≡  ((-1)122315171110)2
  • ≡  (-420)2

and a non-trivial factor is again recovered by gcd(x1x3x7 – 420, n) = q. □

Since numbers with small absolute value are more likely to be B-smooth, most factorization algorithms using this approach are based on a particular method for generating integers x such that x2 mod n has small absolute value. The remainder of this article will describe a few different methods and a way to visualize the generated integers.

Visualizing Numbers

By the Chinese Remainder Theorem there is a one-to-one correspondence between numbers x mod n ∈ {0, …, n – 1} and pairs of numbers (x mod p, x mod q) ∈ {0, …, p – 1} × {0, …, q – 1}. Using this correspondence a number x mod n can be plotted as a circle inside a box of size p×q. The circle is colored blue for positive numbers (i.e., 2(x mod n) < n), and red for negative numbers (i.e., 2(x mod n) ≥ n). The size of the circle increases as the absolute value of x mod n decreases (where the absolute value of a number x mod n is itself if it is positive, and -x mod n if it is negative).

Example: Here is a plot of all numbers x mod n where the integer x satisfies -√n < x < √n. □

Numbers with absolute value at most sqrt(n)

Since the goal is to find integers x such that x2 mod n has small absolute value, the visualization is changed to communicate both pieces of information. The position of the circle plotted for the number x mod n is calculated as above, but the color and size of the circle is calculated from the sign and absolute value of x2 mod n.

Example: Here is a plot of all 1,352 numbers x mod n such that the absolute value of x2 mod n is at most 634. □

All square roots of numbers with small absolute value

The reflective symmetric pattern occurs because every square number x2 mod n is the square of four numbers with circles plotted at the following positions:

  • (x mod p, x mod q)
  • (x mod p, –x mod q)
  • (-x mod p, x mod q)
  • (-x mod p, –x mod q)

This provides a visual interpretation of the core factorization approach: given two of these circles that are not diagonally opposite each other, adding them will produce a point that is on the edge of the box, and the edges consist of all the integers that have a nontrivial common divisor with n.

The blue corners are the result of the integers {1, 2, …} having small positive squares. This makes the bottom left corner blue, and the reflective symmetry means that what happens in one corner happens in them all.

Beyond these features, it is striking how random the visualization appears: finding numbers with small squares is a number theory version of finding a needle in a haystack. The following three sections will each describe an algorithm for searching these kind of haystacks.

Square Roots

The first algorithm for finding integers x such that x2 mod n has small absolute value proceeds by checking integers around √(kn) for small k ∈ {1, 2, …}, since for i ∈ {0, 1, …}

  • |(⌈√(kn)⌉ + i)2 mod n|
  • =  |(√(kn) + (i + (⌈√(kn)⌉ – √(kn))))2 mod n|
  • =  |(kn + 2(i + (⌈√(kn)⌉ – √(kn)))√(kn) + (i + (⌈√(kn)⌉ – √(kn)))2) mod n|
  • =  |(2(i + (⌈√(kn)⌉ – √(kn)))√(kn) + (i + (⌈√(kn)⌉ – √(kn)))2) mod n|
  • ≤  2(i + (⌈√(kn)⌉ – √(kn)))√(kn) + (i + (⌈√(kn)⌉ – √(kn)))2
  • ≤  2(i + 1)√(kn) + (i + 1)2

and similarly

  • |(⌊√(kn)⌋ – i)2 mod n|  ≤  2(i + 1)√(kn) + (i + 1)2

Minimizing the absolute value of x2 mod n is thus achieved by choosing integers x from a range of the form √(kn) – m < x < √(kn) + m.

Example: Here is a plot of all 20,000 integers x in the ranges √(kn) – m < x < √(kn) + m, where k ∈ {1, …, 20} and m = 500.

Numbers of the form floor(sqrt(k * n)) + i

The plot of each range around √(kn) looks like a battleship that is painted blue for the numbers greater than √(kn) and red for the numbers less than √(kn). Unfortunately, none of the numbers on the battleships are marked as larger “hits”, meaning that the absolute value of all of the squares is greater than 634. □

Checking integers close to √(kn) may seem like a simple idea, but the quadratic sieve does exactly this and is the second fastest factorization algorithm after the number field sieve. The quadratic sieve works so well because it defines another sieve algorithm that can search whole ranges at a time to identify all integers x such that x2 mod n are smooth numbers.

Continued Fractions

The second algorithm for finding integers x such that x2 mod n has small absolute value is based on continued fractions. For any real number r, its continued fraction approximation is a sequence of rationals ai/bi each having the property that

  • |ai/bir|  ≤  1/bi2

If a/b is a continued fraction approximation to √n, then

  • |a/b – √n|  ≤  1/b2
  • ⇔  √n – 1/b2  <  a/b  <  √n + 1/b2
  • ⇔  bn – 1/b  <  a  <  bn + 1/b
  • ⇔  b2n – 2√n + 1/b2  <  a2  <  b2n + 2√n + 1/b2
  • ⇔  -2√n + 1/b2  <  a2b2n  <  2√n + 1/b2
  • ⇒  |a2b2n|  <  2√n + 1/b2
  • ⇒  |a2 mod n|  <  2√n + 1/b2

The continued fraction factorization algorithm computes the sequence of continued fraction approximations ai/bi to √n and searches for smooth numbers in ai2 mod n.

Example: Here is a plot of the first 20,000 numerators ai in the sequence of continued fraction approximations to √n.

Numerators of continued fraction approximations to sqrt(n)

The visualization clearly shows two differences between the continued fraction factorization algorithm and the quadratic sieve. Firstly, the generated numbers ai mod n are not clustered in groups, but rather are scattered throughout the whole range {0, …, n – 1}. Secondly, the effect of the tighter bound on the absolute value of ai2 mod n can be observed by the presence of larger circles (each marking an absolute value of at most 634). □

The Pell Equation

The third algorithm for finding integers x such that x2 mod n has small absolute value is based on the Pell equation. Given a positive integer n the Pell equation is a2 = nb2 + 1, and the goal is to find all integers a and b that solve the equation. If n is a square, then every solution must correspond to two integer squares that are distance 1 apart, which can only occur in the trivial solution a = ±1 and b = 0. If n is not a square, then there are infinitely many integer solutions which can all be derived from the smallest nontrivial solution computed by the Chakravala method. The Chakravala method generates a sequence of triples (ai, bi, ki) that satisfy the equation ai2 = nbi2 + ki. Each triple is used to generate the next triple according to a fixed rule which sets up a condition and then finds the minimal |ki| that satisfies the condition. The initial triple is (a0, b0, k0) = (m, 1, m2n) where m is chosen to be whichever of ⌊√n⌋ and ⌈√n⌉ minimizes |k0|. The sequence is guaranteed to converge to the smallest nontrivial integer solution of the Pell equation, which is the values of ai and bi when ki = 1.

Example: Given the Pell equation a2 = 61b2 + 1, here is the sequence computed by the Chakravala method:

  • [(8,1,3), (39,5,-4), (164,21,-5), (453,58,5), (1523,195,4), (5639,722,-3), (29718,3805,-1), (469849,60158,-3), (2319527,296985,4), (9747957,1248098,5), (26924344,3447309,-5), (90520989,11590025,-4), (335159612,42912791,3), (1766319049,226153980,1)]

The smallest nontrivial integer solution a = 1,766,319,049 and b = 226,153,980 can be read from the last triple in the sequence. □

Note that each triple (ai, bi, ki) generated by the Chakravala method satisfies ai2ki (mod n), so the Chakravala method can be turned into a factorization algorithm by generating triples (ai, bi, ki) and searching for smooth numbers in ai2 mod n.

Example: Here is a plot of the first 20,000 integers ai from triples (ai, bi, ki) generated by the Chakravala method.

Integers a_i from triples (a_i,b_i,k_i) generated by the Chakravala method

The visualization is similar to the one for the numerators of the continued fraction approximations, covering the whole range {0, …, n – 1} and containing larger circles marking that the absolute value of the square is at most 634. However, closer examination shows that the Pell equation visualization is overall darker than the continued fraction visualization, indicating that the Chakravala method is finding more small squares than the continued fraction approximations. Perhaps the idea of turning the Chakravala method into a factorization algorithm is worth further study. □

References

All of the computations and images in this article were generated using Haskell.

Advertisements

A Natural GCD Algorithm

The greatest common divisor (gcd) of two non-negative integers is the largest integer that divides both of them, so for example

gcd 15 6 == 3
gcd 21 10 == 1
gcd 10 0 == 10
gcd 0 0 == 0  -- by convention

Two numbers that have a gcd of 1 (like 21 and 10) are said to be coprime. Although in general it’s a hard problem to factor an integer into divisors

91 == 7 * 13
69 == 3 * 13

it’s an easy problem to compute the gcd of two integers

gcd 91 69 == 13

and if they are not coprime then a non-trivial divisor of both is revealed.

Here is a simple recursive definition that efficiently computes gcd:

gcd :: Integer -> Integer -> Integer
gcd a 0 = a
gcd a b = gcd b (a `mod` b)

How would we check that this definition always computes the gcd of two non-negative integers? If we define a divides relation

divides :: Integer -> Integer -> Bool
divides 0 b = b == 0
divides a b = b `mod` a == 0

then we can easily test that gcd returns some common divisor:

commonDivisorTest :: Integer -> Integer -> Bool
commonDivisorTest a b = let g = gcd a b in divides g a && divides g b

But how to check that gcd returns the greatest common divisor? Consider the linear combination g of a and b

g = s * a + t * b

where the coefficients s and t are arbitrary integers. Any common divisor of a and b must also be a divisor of g, and so if g is itself a common divisor of a and b then it must be the greatest common divisor. The extended gcd algorithm computes the coefficients s and t together with the gcd:

egcd :: Integer -> Integer -> (Integer,(Integer,Integer))
egcd a 0 = (a,(1,0))
egcd a b =
    (g, (t, s - (a `div` b) * t))
  where
    (g,(s,t)) = egcd b (a `mod` b)

which allows us to test that egcd returns the greatest common divisor:

greatestCommonDivisorTest :: Integer -> Integer -> Bool
greatestCommonDivisorTest a b =
    let (g,(s,t)) = egcd a b in
    divides g a && divides g b && s * a + t * b == g

The coefficients s and t are useful for much more than checking the gcd result: for any coprime a and b we have

s * a + t * b == 1
(s * a) `mod` b == 1  -- assuming b > 1
(t * b) `mod` a == 1  -- assuming a > 1

which shows that s is the multiplicative inverse of a modulo b, and t is the multiplicative inverse of b modulo a.

Computing Over the Natural Numbers

The extended gcd algorithm is a classic integer algorithm, but what if we wanted a version that computed over the type of natural numbers? This would be useful in environments in which unsigned arithmetic is more natural and/or efficient than signed arithmetic (one such environment being the natural number theory of an interactive theorem prover, which is where this need arose).

The immediate problem is that one of the s and t coefficients returned by the extended gcd algorithm is usually a negative number, which is to be expected since they must satisfy

s * a + t * b == g

where g is generally less than both a and b. The solution is to modify the specification and ask for natural number coefficients s and t that satisfy

s * a == t * b + g

Fortunately for us, it turns out that such coefficients always exist except when a is zero (and b is non-zero).

Unfortunately for us, our modified specification is no longer symmetric in a and b. Given coprime a and b as input, s is the multiplicative inverse of a modulo b as before, but t is now the negation of the multiplicative inverse of b modulo a (the additive inverse of the multiplicative inverse—funny). This asymmetry breaks the recursion step of the extended gcd algorithm, because there is no longer a way to derive coefficients s and t satisfying

s * a == t * b + g

from coefficients s’ and t’ satisfying

s' * b == t' * (a `mod` b) + g

[If you try you’ll discover that a and b are the wrong way around.] This problem is solved by unwinding the recursion so that the algorithm makes two reductions before calling itself recursively, giving coefficients s’ and t’ satisfying

s' * (a `mod` b) == t' * (b `mod` (a `mod` b)) + g

Although this is expression is more complicated, it is now possible to define coefficients

s = s' + (b `div` (a `mod` b)) * t'
t = t' + (a `div` b) * s

which satisfy the desired equation

s * a == t * b + g

Here is the finished implementation:

egcd :: Natural -> Natural -> (Natural,(Natural,Natural))
egcd a 0 = (a,(1,0))
egcd a b =
    if c == 0
    then (b, (1, a `div` b - 1))
    else (g, (u, t + (a `div` b) * u))
  where
    c = a `mod` b
    (g,(s,t)) = egcd c (b `mod` c)
    u = s + (b `div` c) * t

This satisfies a modified test confirming that it produces a greatest common divisor (for non-zero a):

greatestCommonDivisorTest :: Natural -> Natural -> Bool
greatestCommonDivisorTest a b =
    let (g,(s,t)) = egcd a b in
    divides g a && divides g b && s * a == t * b + g

And we can also check the upper bounds on the s and t coefficients (for non-zero a and b at least 2):

coefficientBoundTest :: Natural -> Natural -> Bool
coefficientBoundTest a b =
    let (_,(s,t)) = egcd a b in
    s < b && t < a

Application: The Chinese Remainder Theorem

The Chinese Remainder Theorem states that for coprime a and b, given an x and y satisfying

x < a
y < b

there is a unique n satisfying

n < a * b
n `mod` a == x
n `mod` b == y

But how to compute this n given natural number inputs a, b, x and y? We build on the natural number extended gcd algorithm, of course:

chineseRemainder :: Natural -> Natural -> Natural -> Natural -> Natural
chineseRemainder a b =
    \x y -> (x * tb + y * sa) `mod` ab
  where
    (_,(s,t)) = egcd a b
    ab = a * b
    sa = s * a
    tb = (a - t) * b

The λ expression in the body maximizes the precomputation on a and b, allowing efficient Chinese remaindering of different x and y inputs once the a and b inputs have been fixed. As with all computation over the natural numbers we must be careful with the subtraction operation (a - t) to ensure that the first argument is always larger than the second, but here we know from the coefficientBoundTest that t is always less than a.

References

All of the natural number computations in this blog post have been formally verified and are available as the Haskell package opentheory-divides.

Large Knight’s Tours

I was at middle school when I first heard about the knight’s tour problem: place a knight on any square on a chess board, and move it to every square on the board without visiting any square more than once. A chess knight moves by jumping two squares forward and one square sideways, so a big part of the challenge is getting into and out of all the corners. I remember filling up pages and pages of squared notebooks with attempted solutions, although my record from that time was a tour that covered only 61 of the 64 squares.

Enter the Computer

Well, it’s one thing to look for knight’s tours by hand, but quite another to search for them with a computer. It is actually a fun programming exercise to search for knight’s tours: it is possible to use a standard backtracking algorithm, but it needs a couple of tricks to make it work. Firstly, if you have a choice of squares to jump to, then you should consider them in ascending order of accessibility (i.e., from how many unvisited squares could you jump to this square). This helps to ensure that you don’t waste last chances to visit squares. However, even with this trick, it is possible to make mistakes early on in the search that you can never recover from, and end up lost forever in a maze with 1051 turns and no exits. An ugly way to fix this is to assign a move budget—I chose 100 times the number of squares—and if you don’t find a solution within this many jumps then give up and start again from an empty board.

Armed with a fast computer, a search algorithm and a couple of tricks, I was finally able to realize a childhood dream and find a knight’s tour:

8x8 knight's tour

Note that we always begin our tours in the bottom right corner, so that the search program only has to work out how to get into and out of 3 corners: that time spent filling squared notebooks taught me something!

Building Up

Once the initial problem is solved, I naturally start wondering whether we can find knight’s tours on larger boards. It is easy to modify the search program to change the board dimensions from an 8×8 square to an m×n rectangle, but as the board size increases the time it takes to find a knight’s tour grows rapidly:

Board size (n×n) 8 10 12 14 16 18 20 22 24 26 28 30
Search time (ms) 1 3 7 11 26 51 102 225 401 701 1,074 1,738

Extrapolating these results to a 100×100 board, it seems that my search program would take over a thousand years to find a knight’s tour! Clearly I need a new approach to find knight’s tours for large boards.

Our 8×8 knight’s tour above looks like a little circuit, and drawn the same way a 100×100 knight’s tour would look like a huge circuit. But let’s take some inspiration from the way that huge circuits are designed in practice: not all at once, but rather by wiring together smaller circuits. Consider this special 6×6 knight’s tour:

Special 6x6 knight's tour

One reason this knight’s tour is special is because it has no start and end: it is a loop! There are probably clever ways to find knight’s tours that are loops, but I took a simple-minded approach. I modified my program to keep searching for knight’s tours until it found one where the end was a knight’s jump away from the start, and then I connected the ends of the tour to make a loop.

There is another reason that this knight’s tour is special, which is not so easy to see. You can put two of them together:

Special 6x6 knight's tourSpecial 6x6 knight's tour

then switch two of the jumps around to connect the tours together:

Special 6x6 knight's tourSpecial 6x6 knight's tour

and the result is a knight’s tour on a 12×6 board. To see that the result is indeed a tour, let’s visualize the two separate knight’s tours as loops using a schematic diagram:

Knight's tour loopKnight's tour loop

And now when we perform the connection step the schematic diagram makes it clear that we have created one big loop from the two separate loops:

Knight's tour loopKnight's tour loop

The same thing also works by vertically stacking the special 6×6 knight’s tour:

Special 6x6 knight's tour
Special 6x6 knight's tour

and connecting them like so:

Special 6x6 knight's tour
Special 6x6 knight's tour

Assembling Large Tours

I’ve shown how to connect special 6×6 knight’s tours together, and the next step is to connect them in an m×n grid to make a knight’s tour of a 6m×6n board. Connecting them together in a grid requires a little bit of care to avoid fragmenting the knight’s tour into disconnected paths. Here’s a schematic diagram for a 2×2 grid illustrating how tours can be fragmented into disconnected paths:

Knight's tour loopKnight's tour loop
Knight's tour loopKnight's tour loop

One way to avoid tour fragmentation is to connect along some path through the grid that visits every 6×6 knight’s tour. I chose a spiral, starting at the bottom right of the grid, travelling anti-clockwise and ending in the middle. Here is the schematic diagram for a spiral path on a 3×3 grid:

Knight's tour loopKnight's tour loopKnight's tour loop
Knight's tour loopKnight's tour loopKnight's tour loop
Knight's tour loopKnight's tour loopKnight's tour loop

And here’s the resulting knight’s tour on an 18×18 board:

Special 6x6 knight's tourSpecial 6x6 knight's tourSpecial 6x6 knight's tour
Special 6x6 knight's tourSpecial 6x6 knight's tourSpecial 6x6 knight's tour
Special 6x6 knight's tourSpecial 6x6 knight's tourSpecial 6x6 knight's tour

Filling in the Gaps

At this stage we have the machinery to construct knight’s tours for arbitrarily large boards of dimension 6m×6n. But what about boards of other dimensions, such as the 100×100 board that we set as our target? These boards can be handled by searching for more special tours of dimensions m×n, where 6 ≤ m,n ≤ 11, and using these in place of 6×6 knight’s tours on the bottom and right edges of the grid. For example, here is a knight’s tour of a 23×23 board:

Special 6x6 knight's tourSpecial 6x6 knight's tourSpecial 11x6 knight's tour
Special 6x6 knight's tourSpecial 6x6 knight's tourSpecial 11x6 knight's tour
Special 6x11 knight's tourSpecial 6x11 knight's tourSpecial 11x11 knight's tour

This has the same 3×3 grid pattern as the 18×18 knight’s tour, but with larger pieces around the right and bottom edges. Also, the knight’s tour in the bottom right corner is not a loop—you can see one of its ends in the bottom right corner—so the schematic diagram for the assembled tour ends up looking a little different:

Knight's tour loopKnight's tour loopKnight's tour loop
Knight's tour loopKnight's tour loopKnight's tour loop
Knight's tour loopKnight's tour loopKnight's tour

In fact, it is impossible to find a knight’s tour that is a loop on this board, because it has an odd number of squares. To see why, imagine the squares of the board painted light and dark colors in a checkerboard pattern:

As a knight jumps around the board it alternates between light and dark squares, and a tour of a board with an odd number of squares will contain an even number of jumps. This means the start and end squares of the tour will have the same color, and so there’s no way to jump from the end square of the tour back to the start square to close the loop.

Summary

Using these techniques you can easily assemble knight’s tours for large boards: you can view a tour of a 100×100 board or even assemble your own knight’s tour.

But these techniques only work for boards where the width and height are at least 6: there are also an infinite number of thin boards that might contain knight’s tours:

Can you think of techniques for assembling knight’s tours on thin boards? ♘

Appendix

The performance results were gathered using a MacBook Air with a 1.8GHz Intel Core i7 processor and 4Gb of RAM, running OS X 10.7.5.