charlesreid1.com blog

Computing Square Roots: Part 2: Using Continued Fractions

Posted in Mathematics

permalink

Table of Contents

Continued Fractions

Let's start part 2 of our discussion of computing square roots by talking about continued fractions. When we first learn mathematics, we learn to count in the base 10 system: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9. We can construct representations of all of the integers using these 10 digits, by arranging them in a different order. So, for example, saying 125 is equivalent to saying:

$$ 125 = 1 \times 10^2 + 2 \times 10^1 + 5 \times 10^0 $$

Later on in our mathematical lives, we learn that we can use other integers as our base, or radix, by expressing every integer as the sum of various powers of that integer. For example, 125 can be decomposed into powers of 2. In a base 2 system, we have only two symbols, 0 and 1, so 125 can be represented as 1111101, which is equivalent to saying:

$$ 125 = 1 \times 2^6 + 1 \times 2^5 + 1 \times 2^4 + 1 \times 2^3 + 1 \times 2^2 + 1 \times 2 $$

Note that 1111101 is close to 1111111, which is \(2^7 = 128\).

As Knuth points out in his Art of Computer Programming, Part II:

"If our ancestors had invented arithmetic by counting with their two fists or their eight fingers, instead of their ten "digits," we would never have to worry about writing binary-decimal conversion routines. (And we would perhaps never have learned as much about number systems.)

This idea of an alternative radix to the traditional base 10 leads to entirely new number systems with their own interesting properties. However, it goes even further - as a high school student in 1955, Donald Knuth invented a number system with an imaginary, 4-symbol radix in base \(2i\) (where \(i = \sqrt{-1}\)), called the quater-imaginary base.

There are also number systems that use an irrational radix, such as phinary, which uses the golden ratio as its radix. Thus, the Golden Ratio \(\phi = \dfrac{1 + \sqrt{5}}{2}\) becomes \(\phi = 1\), \(2 = \phi^1 + \phi^-2\) becomes \(2 = 10.01\), \(5 = \phi^3 + \phi^{-1} + \phi^{-4}\) becomes \(5 = 1000.1001\), and so on. In the case of the irrational number \(\sqrt{5}\), this number becomes a rational number in base \(\phi\): \(\sqrt{5} = \phi^2 + \phi^{-1}\) becomes \(\sqrt{5} = 10.1\).

The single central idea behind these number systems is that the abstract set of integers \(\mathbb{Z}\) are being represented on an algebraic ring, which is a fundamental mathematical object. Initially the idea of a ring may seem strange, but it creates the foundations of modern number theory.

Continued Fraction Representations

(Note that while the next two sections will have a lot of "magic numbers," we will step through the procedure in the following sections, and it will be more clear where these numbers come from. There is also an excellent continued fractions calculator available online here.)

We can use other kinds of mathematical objects to represent different numbers. Another technique is to use continued fractions to represent numbers. These are well-studied mathematical objects that date back to Euclid's Elements. The basic idea is to create a recursive expression that involves fractions nested beneath other fractions:

$$ a_0+\cfrac{1}{a_1 +\cfrac{1}{a_2 +\cfrac{1}{ \begin{array}{@{}c@{}c@{}c@{}} a_3 + {}\\ &\ddots\\ &&{}+ \cfrac{1}{a_n} \end{array} }}} $$

For ease of writing, this is often written as:

$$ a_0 + \dfrac{1}{a_1 +} \quad \dfrac{1}{a_2 +} \quad \dfrac{1}{a_3 + } \quad \dots \quad \dfrac{1}{a_{n-1} + } \quad \dfrac{1}{a_n} $$

or, even shorter,

$$ [a_0; a_1, a_2, a_3, \dots] $$

These variables \(a_i\) are called the terms of the continued fraction. Continued fractions can be used to represent rational numbers, in which case the continued fraction representation terminates after a finite number of terms. For example, \(\dfrac{125}{3} = [41; 1, 2]\),

$$ \dfrac{125}{3} = 41 +\cfrac{1}{1 +\cfrac{1}{2}} $$

Continued fractions can also be used to represent irrational numbers, in which case the continued fraction representation is a repeating pattern of variable length. For example, \(\sqrt{14} = [3; \overline{1,2,1,6}]\), where the line over the last digits indicates that the pattern repeats infinitely as \(1, 2, 1, 6, 1, 2, 1, 6, 1, 2, 1, 6, \dots\):

$$ 3 + \cfrac{1}{ 1 + \cfrac{1}{ 2+\cfrac{1}{ 1+\cfrac{1}{ 6+\cfrac{1}{ \begin{array}{@{}c@{}c@{}c@{}} 1 + {}\\ &\ddots\\ &&{} \end{array} } } } } } $$

A few useful properties of these patterns, for square roots, are:

  • First, the integer portion (3 in the case above) is the largest integer whose square is less than the number (14).

  • Second, the sequence of integers that repeats is always palindromic, and it always begins repeating once it reaches a value of \(2 a_0\).

Here are a few more square roots represented as continued fractions, to help illustrate the above properties:

$$ \begin{array} _ \sqrt{19} &=& [4; \overline{2, 1, 3, 1, 2, 8}] \\ \sqrt{115} &=& [10; \overline{1, 2, 1, 1, 1, 1, 1, 2, 1, 20}] \\ \sqrt{988} &=& [31; \overline{2, 3, 4, 1, 20, 6, 1, 14, 1, 6, 20, 1, 4, 3, 2, 62}] \end{array} $$

Next, we'll cover how to turn these terms \([a_0; a_1, a_2, \dots]\) into rational numbers \(\frac{P}{Q}\). These fractions are called convergents.

Convergents

We will show an example of how to compute a continued fraction in the next section, but we will cover one additional topic first. In order to be useful, we need some way to evaluate the continued fraction representation. One way to do this is to repeatedly compute common denominators, perform the fraction addition, and inert the result. The fraction that results is equivalent to the continued fraction expression, but is a rational number and therefore easier to evaluate.

It turns out that this rational number is a very important quantity called the convergent of \(\sqrt{n}\). However, it is cumbersome to perform the operation of de-rationalizing the continued fraction. There is a useful shortcut that takes the form of a recurrence relation.

The \(n^{th}\) convergent is defined as the rational number that results when the continued fraction representation is carried out to \(n\) terms, then expanded and simplified into a rational number. If a number like \(\sqrt{19}\) has an infinite continued fraction representation, it means we can compute progressively more accurate rational approximations.

For example, we know from the above that \(\sqrt{19} = [4; \overline{2,1,3,1,2,8}]\). Using this fact, we can use successive terms to write successive linear approximations:

$$ \begin{array} _ \sqrt{19} &\approx& 4 \\ \quad &\approx& 4 + \frac{1}{2+0} \approx \frac{9}{2} \\ \quad &\approx& 4 + \cfrac{1}{2+\cfrac{1}{1}} \approx \frac{13}{3} \\ \quad &\approx& 4 + \cfrac{1}{2+\cfrac{1}{\cfrac{1}{3}}} \approx \frac{48}{11} \\ \end{array} $$

If we continue this sequence, we get a slew of approximations (there are also additional approximations between each of these terms...):

$$ \begin{array} _ \sqrt{19} &\approx& \frac{61}{14} \\ \quad &\approx& \frac{170}{39} \\ \quad &\approx& \frac{1421}{326} \\ \quad &\approx& \frac{3012}{691} \\ \quad &\approx& \frac{4433}{1017} \\ \sqrt{19} &\approx& \frac{16311}{3742} \end{array} $$

The numerator and denominator of the \(k^{th}\) convergent are denoted \(P_k\) and \(Q_k\), respectively, and can be computed through the recurrence relation:

$$ \dfrac{P_k}{Q_k} = \dfrac{P_{k-2} + a_k P_{k-1}}{Q_{k-2} + a_k Q_{k-1}} $$

where the initial values are \(P_{-2} = 0, P_{-1} = 1\) for P and \(Q_{-2} = 1, Q_{-1} = 0\) for Q (making the first convergent equal to \(\frac{a_0}{1}\)). This can be used to compute successive approximations. Note that on a computational platform, you will quickly reach the end of your accuracy limit, so care must be taken for continued fraction sequences of longer than about 10 terms.

Example: Continued Fraction Coefficients of \(\sqrt{14}\)

Let's walk through an example of how to compute the continued fraction representation \([a_0; a_1, a_2, \dots, a_n]\) for a square root, and how to compute the \(k^{th}\) convergent.

Begin with the square root of a given number. For variety, we will use \(\sqrt{14}\). We begin by computing the nearest integer to 14's square root. If we don't have a square root routine, we can try integers by squaring them, and find the last integer that, when squared, is less than 14. This is \(a_0 = 3\). To interpret, \(a_0\) is the largest integer portion of our square root, with the residual portion \(\sqrt{14}-3\) the portion that will be represented by the continued fraction.

Here's what we have:

$$ \sqrt{14} = 3 + (\sqrt{14} - 3) \\ $$

We can drop the residual for an initial rational approximation of \(\sqrt{14} \approx 13\), but we can do better by going another step.

We recognize that the residual, which we want write in the form \(\dfrac{1}{\mbox{something}}\), can be written as \(\dfrac{1}{\dfrac{1}{\sqrt{14}-3}}\):

$$ \sqrt{14} = 3 + \dfrac{1}{ \frac{1}{\sqrt{14}-3} } $$

Now examine the inverse residual term \(\dfrac{1}{\sqrt{14}-3}\) for step 1, which we will call \(r_1\). Repeat the procedure we performed above: find the nearest integer portion to this quantity. In this case,

$$ \dfrac{1}{r_1} = \dfrac{1}{\sqrt{14}-3} = 1.34833... $$

Now split this into its integer portion, which becomes \(a_1\), and its remaining fractional portion, .34833...:

$$ a_1 = floor( \dfrac{1}{r_1} ) $$

and the new residual \(r\) is written in terms of the old residual and the coefficent \(a_i\):

$$ r_2 = \dfrac{1}{r_1} - a_1 $$

Generalizing, we get an iterative procedure to determine the coefficients \(a_i\):

$$ a_i = floor( \dfrac{1}{r_i} ) $$

followed by:

$$ r_{i+1} = \dfrac{1}{r_i} - a_i $$

where the initial values \(a_0, r_0\) are computed as mentioned above, and the rest of the values in the sequence follow.

Continuing for \(\sqrt{14}\), we get:

$$ \begin{array} _ a_0 &=& 3 \\ r_1 &=& \dfrac{1}{\sqrt{14}-3} = 1.348331 \\ a_1 &=& 1 \\ r_2 &=& \dfrac{1}{.348331} = 2.870829 \\ a_2 &=& 2 \\ r_3 &=& \dfrac{1}{0.870829} = 1.1483311 \\ a_3 &=& 1 \\ r_4 &=& \dfrac{1}{0.1483311} = 6.741676 \\ a_4 &=& 6 \end{array} $$

It is at this point that we see \(2 a_0\) and know that our (palindromic) sequence will repeat. (When we evaluate the convergents, we will utilize the palindromic nature of this sequence.)

Collecting these terms gives us the expected result: \(\sqrt{14} = [3; 1, 2, 1, 6]\).

This gives us an algorithmic procedure for computing the continued fraction representation \([a_0; a_1, a_2, a_n]\) of a number - with the important caveat, as mentioned above, that some integers have sequences that are extremely long before they repeat, making it impossible to find the full continued fraction representation without arbitrary precision libraries.

That being said, if you are only interested in the first 10 or so terms of the continued fraction representation, here is a static Java method to compute them:

    /** Find the (shorter than 10) continued fraction sequence for sqrt(n). 
     * This returns a list of integers, [a0, a1, a2, ...]
     *
     * @params n Number to compute the square root of.
     */
    public static List<Integer> continuedFractionSqrt(int n) {
        if(isPerfectSquare(n)) {
            throw new IllegalArgumentException("Error: n cannot be a perfect square.");
        }
        int niters = 10; // handbrake

        int ai = 0;
        double val = 0;
        double remainder = 0;
        List<Integer> coeffs = new ArrayList<Integer>();

        // Fencepost
        remainder = 1.0/Math.sqrt(n);

        for(int i=0; i<niters; i++) {
            val = 1.0/remainder;
            ai = floor(val);
            remainder = val - ai;
            coeffs.add(ai);
            if(coeffs.get(i) == 2*coeffs.get(0)) {
                break;
            }
        }
        return coeffs;
    }

    /** Check if x is a perfect square. */
    public static boolean isPerfectSquare(int x) { 
        int s = (int)(Math.round(Math.sqrt(x)));
        return x==(s*s);
    }

    /** Find the floor of a double. */
    public static int floor(double j) {
        return (int)(Math.floor(j));
    }

Here is a short program that uses this routine to compute the continued fraction representation of \(\sqrt{14}\):

public class SquareRootCF {
    public static void main(String[] args) {
        System.out.println("sqrt(14) = "+Convergents.continuedFractionSqrt(14));
        System.out.println("sqrt(19) = "+Convergents.continuedFractionSqrt(19));
    }
}

and the corresponding output:

$ javac SquareRootCF.java && java SquareRootCF
sqrt(14) = [3, 1, 2, 1, 6]
sqrt(19) = [4, 2, 1, 3, 1, 2, 8]

Example: Convergents of \(\sqrt{14}\)

We now turn to the task of computing the convergents of the continued fraction, which will yield successive rational numbers that are progressively better approximations to \(\sqrt{14}\).

We start with the expression given above for the \(k^{th}\) convergent:

$$ \dfrac{P_k}{Q_k} = \dfrac{P_{k-2} + a_k P_{k-1}}{Q_{k-2} + a_k Q_{k-1}} $$

with initial values \(P_{-2} = 0, P_{-1} = 1\) for P and \(Q_{-2} = 1, Q_{-1} = 0\) for Q. This yields the first "real" convergent:

$$ \dfrac{P_1}{Q_1} = \dfrac{1 + 1 \times a_1}{0 + a_1 \times 1} = \dfrac{1 + 3}{0 + 1} = \dfrac{4}{1} $$

Successive approximations will use the values \(P_1\) and \(Q_1\) to compute the next convergents.

$$ \dfrac{P_2}{Q_2} = \dfrac{P_0 + P_1 a_1}{Q_0 + Q_1 a_1} = \frac{11}{3} $$

Continuing in this fashion gives:

$$ \begin{array} \quad \dfrac{P_2}{Q_2} &=& \dfrac{15}{4} \\ \dfrac{P_3}{Q_3} &=& \dfrac{101}{27} \end{array} $$

and so on. This recurrence relation is easy to code up. It starts with the continued fraction coefficients for the given square root, and computes successive values of P and Q. The number of terms computed is specified by the user. Once it reaches the end of the sequence of continued fraction coefficients, it can start at the beginning again (the sequence is palindromic).

Finally, it returns the values of \(P_k\) and \(Q_k\), and of the successive convergents:

$$ \begin{array} _ \sqrt{14} &\approx& 3 \\ &\approx& 4 \\ &\approx& 11/3 \\ &\approx& 15/4 \\ &\approx& 101/27 \\ &\approx& 116/31 \\ &\approx& 333/89 \\ &\approx& 449/120 \\ &\approx& 3027/809 \\ &\approx& 3476/929 \\ &\approx& 9979/2667 \end{array} $$

Here is a static method in Java that will compute the convergents of a square root:

    /** 
     * Compute the convergents (rational representation of
     * continued fraction).
     *
     * This uses the recurrence relation:
     *
     * P_n     P_n-2 + a_n P_n-1
     * ---- = -----------------
     *  Q_n    Q_n-2 + a_n Q_n-1
     */
    public static long[] convergents(int n, int nterms) {
        long[] convergents = new long[2];

        List<Integer> cfrepr = continuedFractionSqrt(n);

        // Initial values for convergent recurrence relation
        long Pnm2 = 0; // P_{n-2}
        long Pnm1 = 1;
        long Qnm2 = 1;
        long Qnm1 = 0;
        long P = 0;
        long Q = 0;

        // Term 0 is the constant value a0.
        int accessindex = 0;
        for(int i=0; i<=nterms; i++) { 
            int an = cfrepr.get(accessindex);

            P = Pnm2 + an * Pnm1;
            Q = Qnm2 + an * Qnm1;

            Pnm2 = Pnm1;
            Pnm1 = P;

            Qnm2 = Qnm1;
            Qnm1 = Q;

            if(accessindex+1>=cfrepr.size()) { 
                // Ensure we keep repeating the sequence
                // if the sequence has fewer terms than
                // the user asks for.
                // This allows us to get really good
                // approximations for numbers.
                // This only works because the sequence
                // is palindromic.
                accessindex = 1;
            } else {
                accessindex++;
            }
        }

        convergents[0] = P;
        convergents[1] = Q;

        return convergents;
    }

Here is a simple driver program that prints out several convergents for \(\sqrt{14}\):

public class SquareRootCF {
    public static void main(String[] args) {
        for(int i=1; i<=10; i++) {
            long[] conv = Convergents.convergents(14,i);
            System.out.println("Convergent "+i+": "+conv[0]+"/"+conv[1]);
        }
    }
}

and the corresponding console output:

$ javac SquareRootCF.java && java SquareRootCF
Convergent 1: 4/1
Convergent 2: 11/3
Convergent 3: 15/4
Convergent 4: 101/27
Convergent 5: 116/31
Convergent 6: 333/89
Convergent 7: 449/120
Convergent 8: 3027/809
Convergent 9: 3476/929
Convergent 10: 9979/2667

References

  1. Hardy, G. H. A Course of Pure Mathematics. Cambridge University Press, Tenth Edition (1908-1967).

  2. Knuth, Donald. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. Addison-Wesley Publishing Company, Second Edition (1975).

Computing Square Roots: Part 1: Using Newton's Method

Posted in Mathematics

permalink

Table of Contents

Newton's Method for Finding Roots of Equations

Suppose we have a function \(f(x)\) and we want to compute values of \(x\) for which \(f(x)=0\). These values of \(x\) are called the roots of \(f(x)\).

We can compute the roots using Newton's Method, which utilizes the derivative of the function to iteratively compute the roots of the function.

Newton's method begins with an initial guess. It evaluates the derivative of the fnction at the initial guess, which gives the slope of the tangent line, and computes the root of the tangent line as the next approximation of the root of the function.

This is based on the point-slope formula,

$$ y - y_0 = m(x - x_0) $$

Now the left side becomes

$$ f(x) - f(x_0) = f'(x_0)(x - x_0) $$

and at the root, we know \(f(x)=0\), so rearranging this equation into an expression for the root \(x\) gives:

$$ x = x_0 - \dfrac{f(x_0)}{f'(x_0)} $$

Now, this is the equation for the next approximation for the root. To turn this into an iterative procedure, we write this as

$$ x_{i+1} = x_{i} - \dfrac{f(x_i)}{f'(x_i)} $$

Newton's Method then allows us to evaluate the above expression as many times as we would like to achive the desired accuracy.

Two caveats with Newton's Method:

  • The function must be relatively well-behaved; Newton's Method does not converge for functions with large, high-order derivatives.
  • The convergence of the method depends on the accuracy of the initial guess. If you can make a good guess, do it!

Newton's Method for Finding Square Roots

Note that we can use the procedure and equation from above to compute the numerical value of a given function to an arbitrary degree of accuracy, so long as we have the derivative (and a program that can keep track of arbitrary digits).

Suppose we want to use Newton's Method to compute \(\sqrt{2}\). Then we can solve for the roots of the following function:

$$ f(x) = x^2 - a $$

for \(a = 2\).

Note that the derivative of this function is computed using the power rule, which is trivial to implement, so we can also use this method to compute general \(n^{th}\) roots of \(a\), by solving for roots of:

$$ f(x) = x^n - a $$

For a monomial with nonzero power, the derivative is always l\(f'(x) = n x^{n-1}\), so in the square root case we have \(f'(x) = 2x\). Now we plug these two functions into the iterative formula for Newton's Method to get:

$$ \begin{array} a x_{i+1} &=& x_{i} - \dfrac{f(x_n)}{f'(x_n)} \\ x_{i+1} &=& x_{i} - \dfrac{x_i^2 - a}{2 x_i } \end{array} $$

Newton's Method for Finding Square Roots: Program

This can be implemented in a programming language to yield an iterative method for computing square roots. We can either specify a tolerance (better), or a number of iterations. Here is a static method for computing a square root using Newton's method. The user specifies the value of \(a\), the number they want to compute the square root of; \(x_0\), the initial guess; and a tolerance, which controls the number of decimal places of accuracy of the final answer.

    /** Compute a square root using Newton's Method, to a specified tolerance.
     *
     * @param a Compute the square root of a.
     * @param x0 Initial guess.
     * @param tol Tolerance (stopping criteria).
     */
    public static double nmsqrttol(double a, double x0, double tol) { 
        double xi, xip1, err;
        xi = xip1 = x0;
        err = x0;
        while(err > tol) { 
            xip1 = xi - (xi*xi-a)/(2*xi);
            err = Math.abs( xip1 - xi );
            xi = xip1;
        }
        return xip1;
    }

(There is no check for infinite loops because our function is a smooth polynomial and Newton's Method will always converge.)

Here is sample output from a program that varies the tolerance and prints the corresponding value of the square root that was computed:

$ javac SquareRoot.java && java SquareRoot
Actual value sqrt(2) = 1.4142135623730951
Testing Newton's Method, Specifying Tolerance:
Tol         sqrt(2)
0.1         1.4166666666666667
0.01        1.4142156862745099
0.001       1.4142135623746899
0.0001      1.4142135623746899
1e-05       1.4142135623746899
1e-06       1.4142135623730951
1e-07       1.4142135623730951
1e-08       1.4142135623730951

Accuracy

Now that we've coded up Newton's Method, we can determine the number of accurate digits. There are a couple of ways to do this, but I went with a string comparison method. I start with a text file containing thousands of digits of the square root of 2, then I compute the square root of 2 using Newton's Method. The longest common substring gives me the number of accurate digits, plus "1.", so if I subtract 2 I get the total number of accurate digits after the decimal place.

javac SquareRoot.java && java SquareRoot
Tolerance = 0.1     Number of accurate digits = 2
Tolerance = 0.01        Number of accurate digits = 5
Tolerance = 0.001       Number of accurate digits = 11
Tolerance = 0.0001      Number of accurate digits = 11
Tolerance = 1e-05       Number of accurate digits = 11
Tolerance = 1e-06       Number of accurate digits = 15
Tolerance = 1e-07       Number of accurate digits = 15
Tolerance = 1e-08       Number of accurate digits = 15
Tolerance = 1e-09       Number of accurate digits = 15
Tolerance = 1e-10       Number of accurate digits = 15
Tolerance = 1e-11       Number of accurate digits = 15
Tolerance = 1e-12       Number of accurate digits = 15
Tolerance = 1e-13       Number of accurate digits = 15
Tolerance = 1e-14       Number of accurate digits = 15

Speed

It is also important to measure performance, in the form of speed. How fast is Newton's Method relative to the built-in square root function in the math library?

Turns out the performance of Newton's Method is much worse than the built-in math library's square root function. The Newton's Method defined above is around 100 times slower. Here are the results of a simple timing test, in which we compute the square root 10 million times, timing how long it takes, and divide by the number of operations to yield the time per operation (or rather, as it is slightly easier to grasp, the time per 1k operations);:

    /** Time Newton's Method.
     *
     * How long does it take to achieve 10 digits of accuracy? */
    public static void testTime() { 

        int Nops;
        double a;
        double initialGuess;
        double tol;
        double time;

        Tim tim;

        Nops = 10000000;
        a = 2;
        initialGuess = 1;
        tol = 1E-8;
        tim = new Tim();
        tim.tic();
        for(int i=0; i<Nops; i++) { 
            nmsqrttol(a, initialGuess, tol);
        }
        tim.toc();
        time = 1000*tim.elapsedms()/Nops;
        System.out.println("Newton's Method Time (ms) per 1k operations: "+time);


        Nops = 10000000;
        a = 2;
        tim = new Tim();
        tim.tic();
        for(int i=0; i<Nops; i++) { 
            Math.sqrt(a);
        }
        tim.toc();
        time = 1000*tim.elapsedms()/Nops;
        System.out.println("Math Library Time (ms) per 1k operations: "+time);

    }

and the results:

javac SquareRoot.java && java SquareRoot
Newton's Method Time (ms) per 1k operations: 0.016
Math Library Time (ms) per 1k operations: 3.0E-4

While the accuracy of Newton's Method for square roots may not be that great, it ain't bad, for 9 lines of code.

Next Steps

In a blog post to follow, we'll compare the speed and accuracy of square root computations using Newton's Method to an alternative approach involving the continued fraction representation of square roots. This particularly interesting technique can also be used to solve the Pell equation, a quadratic Diophantine equation of the form:

$$ x^2 - D y^2 = 1 $$

But more on that in a future post...

References

  1. Press, William et al. "Numerical Recipes in C." Cambridge Unviersity Press (2007).

  2. "SquareRoot.java". Charles Reid. git.charlesreid1.com. <https://git.charlesreid1.com/cs/java/src/master/numerics/newtons-method>

CSE 143 Final Project: Hilbert Sort: 3. The Code

Posted in Computer Science

permalink

Table of Contents

This is the third in a series of three posts detailing the Hilbert Sort problem, its solution, and its implementation. This post deals with the code to solve the Hilbert Sort problem.

Hilbert Sort: Pseudocode

From our prior post, here is the psudocode for our Hilbert Sort function:

define hilbert_sort( unsorted queue, square dimension ):
    create southwest queue
    create northwest queue
    create northeast queue
    create southeast queue
    for each point:
        if in southwest:
            create new point using X -> Y, Y -> X
            add to southwest queue
        if in northwest:
            create new point using X -> 2X, Y -> 2Y - S
            add to northwest queue
        if in northeast:
            create new point using X -> 2X - S, Y -> 2Y - S
            add to northeast queue
        if in southeast:
            create new point using X -> S - 2Y, Y -> 2S - 2X
            add to southeast queue

        hilbertsort(southwest queue, square dimension)
        hilbertsort(northwest queue, square dimension)
        hilbertsort(northeast queue, square dimension)
        hilbertsort(southeast queue, square dimension)

        create new results queue
        add points from southwest into results queue
        add points from northwest into results queue
        add points from northeast into results queue
        add points from southeast into results queue
        return results queue

Because we are manually sorting, and we want order to be preserved, we should be using a queue to organize points as we sort them. That way, we add them in sorted order, and we are then able to remove them in sorted order.

Hilbert Sort: Code

We begin by covering a utility class used by the Hilbert Sort method to store \((X,Y)\) points. This is a simple example of a composition design pattern. Next, we cover the bulk of the problem solution: the recursive sort method that partiions points into quadrants. Finally, we cover the main method, which demonstrates reading data from an input file and passing it to the sort method.

Hilbert Sort: Utility Classes

To organize \((X,Y)\) point data, we use a simple class using composition. This is defined next to the HilbertSort class.

/**
 * An (x,y) Point class. 
 */
class Point {
    int x, y; // (x,y) point.
    String name; // Each (x,y) point has a name in the file. Used for output.
    /** Constructor. */
    public Point(int x, int y, String name) { 
        this.x = x; this.y = y; this.name = name;
    }
    /** String representation (x,y). */
    public String toString() { 
        return "("+this.x+","+this.y+","+this.name+")";
    }
}

Recursive Sort Function

Following is the recursive sort method, which (like merge sort) consists of a split step, which partitions an \(S \times S\) square into quadrants and distributes points in the square into their corresponding quadrants, and a merge step, which stitches together each quadrant in the correct order.

    /** Recursive implementation of a Hilbert sort. */
    public static Queue<Point> hilbertSort(Queue<Point> inputP, int S) {
        // Recursive method:
        // Apply the Hilbert geometrical quadrant division 
        // to sort points by when they are visited by a Hilbert curve.
        //
        // Base case: 
        // There are 1 or fewer points in each quadrant.
        // Keep splitting into quadrants until we reach the base case. 
        if(inputP.size()<1) {
            return new LinkedList<Point>();
        } else if(inputP.size()==1) {
            return inputP;
        }

        // split by quadrant
        Queue<Point> qSW = new LinkedList<Point>();
        Queue<Point> qNW = new LinkedList<Point>();
        Queue<Point> qNE = new LinkedList<Point>();
        Queue<Point> qSE = new LinkedList<Point>();

        // Sort points by dividing into quadrants
        for(Point p : inputP) { 

            // Prepare for the tricky part.
            //
            // Rotate the quadrant, and points in it,
            // so that everything is now translated to fit
            // how the template of the Hilbert curve is being drawn.
            // (SW->NW->NE-SE)

            boolean inSWquadrant = (2*p.x <= S) && (2*p.y <= S);
            boolean inNWquadrant = (2*p.x <= S) && (2*p.y >= S);

            boolean inNEquadrant = (2*p.x >= S) && (2*p.y >= S);
            boolean inSEquadrant = (2*p.x >= S) && (2*p.y <= S);

            // Each time we sort (x,y) points into quadrants,
            // we also transform each coordinate point 
            // in such a way that it rescales to an S x S square,
            // but does not modify the order of the points. 
            //
            // Note that we can keep everything as integers by
            // continuing to look at an S x S square,
            // and double the x and y values to shift them over/up.
            //
            // Two easy cases:
            if(inNWquadrant) {
                // Northwest quadrant: 
                // - shift y down by S/2
                // - keep x and y in same order
                qNW.add( new Point(2*p.x, 2*p.y-S, p.name) );

            } else if(inNEquadrant) {
                // Northeast quadrant:
                // - shift x and y down by S/2
                // - keep x and y in same order
                qNE.add( new Point(2*p.x - S, 2*p.y - S, p.name) );

            } else if(inSWquadrant) { 
                // Southwest quadrant:
                // - x and y need to swap places 
                // - that's it.
                qSW.add( new Point(2*p.y, 2*p.x, p.name) );

            } else if(inSEquadrant) { 
                // Southeast quadrant:
                // - trickiest quadrant.
                // - We want to preserve S - x, distance from right side
                // - we want to use it as the new y coordinate
                qSE.add( new Point(S - 2*p.y, 2*(S - p.x), p.name) );

            }

        }
        // Sort til you reach the base case.
        qSW = hilbertSort(qSW, S); 
        qNW = hilbertSort(qNW, S); 
        qNE = hilbertSort(qNE, S); 
        qSE = hilbertSort(qSE, S);

        Queue<Point> result = new LinkedList<Point>();
        for(Point q : qSW) result.add(q);
        for(Point q : qNW) result.add(q);
        for(Point q : qNE) result.add(q); 
        for(Point q : qSE) result.add(q);

        return result;
    }

Main Method

The last part of the code is the portion that loads the points and their labels from a file, and populates a Queue of Point objects from it. This queue of points is then sorted and returned in order.

    /** Main driver. */
    public static void main(String[] args) { 

        Scanner stdin = new Scanner(new BufferedReader(new InputStreamReader(System.in)));

        int n = stdin.nextInt();
        int S = stdin.nextInt();

        // n lines of 3 tokens each
        Queue<Point> inputPoints = new LinkedList<Point>();
        for(int i=0; i<n; i++) { 
            int x0 = stdin.nextInt();
            int y0 = stdin.nextInt();
            String label = stdin.next();
            inputPoints.add(new Point(x0,y0,label));
        }
        Queue<Point> sortedPoints = hilbertSort(inputPoints, S);
        for(Point p : sortedPoints) { 
            System.out.println(p.name);
        }
    }

References

  1. "ACM Pacific Region Programming Competition." Association of Computing Machinery. Accessed 19 June 2017. <http://acmicpc-pacnw.org/>

  2. "Hilbert Sort." Git repository, git.charlesreid1.com. Charles Reid. Updated 16 June 2017. <https://git.charlesreid1.com/cs/finalproject-143/src/master/hilbert/HilbertSort.java>

March 2022

How to Read Ulysses

July 2020

Applied Gitflow

September 2019

Mocking AWS in Unit Tests

May 2018

Current Projects