Programming Praxis

We examined the bifid cipher, which uses a polybius square to convert letters to digits and back again, in a previous exercise. Today we will look at another tool for converting between letters and digits known as a straddling checkerboard. A sample straddling checkerboard, based on the key SHARPEN YOUR SAW with spaces at 2, 5, and 9, is shown below (the three underscores in the first row make the space character visible):

  0 1 2 3 4 5 6 7 8 9
  S H _ 8 A _ 1 R P _
2 E 5 N * Y O U W B 2
5 C 3 D 4 F 6 G 7 I 9
9 J 0 K L M Q T V X Z

The checkerboard has four rows and ten columns. Three of the positions in the first row are spaces, leaving twenty-six letters, ten digits, and an asterisk to represent the space character. The row numbers correspond to the positions of the three first-row spaces. Our method is traditional: the pass phrase followed by the alphabet, with duplicates removed and digits paired with the first ten letters of the alphabet (A=1, B=2, …, J=0). But any other method of filling out the key may be used.

To use the checkerboard, simply write one digit for letters that appear in the first row and two digits for letters in the subsequent rows, row digit first then column digit. For instance, the plain-text PROGRAMMING PRAXIS is converted to digits like this:

P R O   G   R A M   M   I   N   G   *   P R A X   I   S
8 7 2 5 5 6 7 4 9 4 9 4 5 8 2 2 5 6 2 3 8 7 4 9 8 5 8 0

Then the digits are processed in some way. A traditional method is double transposition, but we will use non-carrying (modulo 10) addition with key 2641:

8 7 2 5 5 6 7 4 9 4 9 4 5 8 2 2 5 6 2 3 8 7 4 9 8 5 8 0
2 6 4 1 2 6 4 1 2 6 4 1 2 6 4 1 2 6 4 1 2 6 4 1 2 6 4 1
0 3 6 6 7 2 1 5 1 0 3 5 7 4 6 3 7 2 6 4 0 3 8 0 0 1 2 1

Then the sums are converted back to letters and digits using the checkerboard in reverse. Although some letters are represented by one digit and other letters are represented by two digits, there is no ambiguity since the leading digits are known:

0 3 6 6 7 2 1 5 1 0 3 5 7 4 6 3 7 2 6 4 0 3 8 0 0 1 2 1
S 8 1 1 R 5   3   S 8 7   A 1 8 R U   A S 8 P S S H 5

Thus, the completed cipher-text is S811RS3S87A18RUAS8PSSH5. Decryption is the inverse operation.

Your task is to write functions that perform encryption and decryption using the straddling checkerboard as shown above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

[ There is a bug in the solution of this exercise. See the revised version of the exercise for a proper solution. ]

We examined the Miller-Rabin probabilistic primality checker in a previous exercise. Today, we examine a primality checker that combines the Miller-Rabin test with a test on Lucas pseudoprimes, devised by Robert Baillie and described by Baillie and Wagstaff in their article “Lucas Pseudoprimes” in Mathematics of Computation, Volume 35, Number 152, pages 1391-1417, October 1980; see also Thomas Nicely’s web page devoted to Baillie’s test. This is the same algorithm used in the PrimeQ function in Mathematica.

Lucas numbers are defined by a recurrence formula similar to Fibonacci numbers, where with and ; the Lucas numbers are 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, … (Sloane’s A000204). Lucas numbers have the rather startling property that, if n is prime, . But the converse is not true, and composite numbers n such that are known as Lucas pseudoprimes; the first few Lucas pseudoprimes are 705, 2465, 2737, 3745, 4171, … (Sloane’s A005845).

Lucas numbers are a special case of Lucas sequences. If P and Q are integers such that the discriminant , then the roots of are and . There are two Lucas sequences and for , which can be computed by the recurrence equations and . The Lucas numbers are given by the sequence and the Fibonacci numbers are given by the sequence .

We will want to compute the nth element of a Lucas sequence, mod n, for large n. Rather than computing the entire recurrence in time proportional to n, it is possible to use doubling and halving, in the same way as the exercise on Three Binary Algorithms, to compute the nth element in log n time. Such a computation is known as a Lucas chain.

Thus, to test whether an odd number n is a Lucas pseudoprime, we choose sequence parameters P and Q so that the the discriminant D is non-square and the Legendre number (otherwise the modular arithmetic would fail). Then we construct the Lucas chain; if the nth element is zero modulo n, then n is either prime or is a Lucas pseudoprime.

Recall that the Miller-Rabin test used strong pseudoprime tests on fifty bases to check primality. Using Lucas pseudoprimes, we can reduce the number of tests substantially. It turns out that combining two strong pseudoprime tests, with bases 2 and 3, with a Lucas pseudoprime test, there are no known pseudoprimes; the test has been performed exhaustively on all numbers less than 10¹⁶, and no counter-examples have been found in over twenty-five years of use (in addition to fifteen minutes of fame, there is a monetary reward from the original authors for the finder of a counter-example).

Your task is to write a function that checks primality using the algorithm described above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

The weather forecast in the daily newspaper publishes the times of sunrise and sunset, a calculation that we examined in a previous exercise. The daily newspaper also publishes the phase of the moon, a calculation that we examine in today’s exercise.

The moon circles the earth every 29.530588853 days, so pick a starting point and count. A new moon occurred at julian date 2451550.1 (January 6, 2000). Then it is easy to count the number of days since the last new moon.

There are eight generally recognized phases of the moon: new, waxing crescent, first quarter, waxing gibbous, full, waning gibbous, last quarter, and waning crescent. To calculate the phase of the moon simply divide the days since the last new moon by eight and select the appropriate phase.

Your task is to write a function that calculates the phase of the moon for a given date. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

[ We studied methods of computing sines in a previous exercise. In today’s exercise, guest author Jos Koot, a retired physicist, amateur pilot, and Scheme enthusiast, gives us an example of the use of sines. Feel free to contact me if you have an idea for an exercise, or if you wish to contribute your own exercise. Jos writes this disclaimer: Do not use the example program for real life flight planning. It may not be fool proof!!! Use more appropriate equipment and procedures after studying aeronautics more profoundly. I do not accept any claims. There is a lot more to flying than presented in this exercise. ]

When preparing to fly, a pilot must know:

• The direction and distance from the airfield of departure to the airfield of destination. This direction is called the ground track (GT).
• The name (WN) and speed (WS) of the wind. The name is the direction from which the wind blows and is measured in degrees clockwise from north. The direction proper (WD) is opposite to that of the name.
• The air speed of the plane (AS) This is the speed of the airplane relative to the moving air.

Directions are measured in degrees clockwise from geographic north, 360 degrees on the compass-card. You may measure distances and times in any unit you like. Speeds should be in unit of distance per unit of time, of course. We take the air speed as a given datum and we ignore time needed for climbing and descending. Because the plane may be drifted by the wind, the pilot must calculate which direction to head for in order to compensate the drift. For this purpose the pilot may draw a triangle whose sides are speeds:

• Ground track. This side is to be drawn in the direction from takeoff to destination but with unknown speed. This still unknown speed is called the ground speed (GS).
• Drift caused by the wind. Both direction and speed are known thanks to meteorologists. The drift points in the direction the wind blows to and is opposite to that of the name of the wind. For example, wind from east causes drift to west.
• True heading (HT). This is the direction the pilot must choose in order to compensate for the drift. The direction is unknown, but the speed is: it is the air speed, also called the cruise speed.

It is instructive first to draw some navigation triangles by means of a pair of compasses and a ruler. The air speed is significantly greater than the wind speed, of course. First draw a line in the direction of the ground track. On this line choose a point B. Draw the drift ending in point B. Now we have point C. Draw a circle with point C for its centre and the air speed for its radius. The circle has two points of intersection with the ground track. One of the intersection points is A. Point A and the other point of intersection correspond to the two roots of the quadratic equation. Because the ground speed is directed from point A to point B it is not difficult to see which point of intersection is the right one. Now the vector equation AB = AC + CB is satisfied.

The known data are: ground track GT, distance (D) from departure to destination, wind name WN, wind velocity WS and the air speed AS. The air speed is significantly greater than the velocity of the wind, of course. From these data we want to compute: true heading TH, ground speed GS and time of the flight (FT=D/GS).

One solution uses the law of sines, which states that the ratio of the length of a side to the sine of the opposite angle is the same for all three sides of any triangle. We have , hence . The ground speed can be calculated as .

Another solution uses the cosine rule c² = a² + b² – 2ab cos(θ) to compute the third side of a triangle when two sides and the included angle are known, a quadratic equation is obtained. Solving this equation for the unknown ground speed, we have . Because we want a positive ground speed, take the root with the plus sign. Finally use the cosine rule again in order to compute the angle of correction (angle A), which is the angle between the ground track and the required true heading: , so the true heading is . This angle must be added or subtracted from the direction of the ground track depending on which side the wind drifts the plane.

Your task is to write two functions, as above, that take a distance, ground track, wind speed and direction, and air speed and computes the true heading, ground speed, correction and duration of flight. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Today’s exercise examines three binary algorithms, for basic arithmetic: multiplication, division, and greatest common divisor.

Four millenia ago, when the ancient Egyptians were building the pyramids, they invented a method for multiplying positive integers that works like this: They start by writing the two numbers to be multiplied at the head of two columns. Then they repeatedly divide the number in the left column by two (right-shift), dropping any remainder, and double the number in the right column (left-shift), writing the two new numbers immediately below their predecessors, until the number in the left column is one. Then they cross out all rows where the number in the left column is even, and add the remaining numbers in the right column, which is the desired product. For instance, the product fourteen times twelve is found like this:

14  12
 7  24
 3  48
 1  96
   168

It is easy to see why this method works if you use the grade-school method of multiplication, but with binary numbers instead of decimal numbers:

        1 1 0 0                     1 2
        1 1 1 0                     1 4
        - - - -                     - -
        0 0 0 0       0 x 1 2         0
      1 1 0 0         2 x 1 2       2 4
    1 1 0 0           4 x 1 2       4 8
  1 1 0 0             8 x 1 2       9 6
  - - - - - - -     - -   - -     - - -
1 0 1 0 1 0 0 0     1 4 x 1 2     1 6 8

In binary, 14 is 1 · 2³ + 1 · 2² + 1 · 2¹ + 0 · 2⁰. The odd numbers in the left column correspond to the 1 bits in the binary representation of the multiplicand.

Binary division is the inverse operation to binary multiplication, and is also done by a series of doubling and halving operations. We’ll begin with the grade-school algorithm, in binary:

                        1 0 0 1 1
            - - - - - - - - - - -
1 0 1 0 1 1 ) 1 1 0 1 0 0 0 1 0 1           8 3 7
              1 0 1 0 1 1                   6 8 8 / 4 3 = 1 6
              - - - - - -                   - - -
                  1 0 0 1 0 1 0             1 4 9
                    1 0 1 0 1 1               8 6 / 4 3 =   2
                  - - - - - - -             - - -
                      1 1 1 1 1 1             6 3
                      1 0 1 0 1 1             4 3 / 4 3 =   1
                      - - - - - -             - -
                        1 0 1 0 0             2 0

The quotient has a 1 wherever the divisor is less than the current portion of the dividend, and a 0 where the next digit is “brought down” from the dividend. The algorithmic version of n ÷ d is based on a quantity d · 2^k, which we call t, which is initialized with the largest k that leaves t less than n; this can be done by repeated left shifts. Then the quotient q is initialized to zero, the remainder r is initialized to n, and an iteration continues halving t until t is less than the original d; at each stage of the iteration, the quotient is doubled (left-shift) if the current remainder is less than t, otherwise the quotient is doubled (left-shift), then one is added, and the current remainder is reduced by the current value of t. An example is shown below:

t    q  r
43
86
172
344
688
1376
688  0  837
344  1  149
172  2  149
86   4  149
43   9  63
43/2 19 20

Josef Stein described a binary algorithm for greatest common divisor in 1967 that was devised by Chinese mathematicians in the first century. The greatest common divisor of two numbers is calculated by a recursive algorithm with four rules:

• If either number is zero, their greatest common divisor is the other number.
• If both numbers are even, the two numbers share a factor of two, so their greatest common divisor is double (left-shift) the greatest common divisor of the halves (right-shift) of the two numbers.
• If one number is even and the other is odd, they have no factor of two in common, so their greatest common divisor is the greatest common divisor of the odd number and half (right-shift) the even number.
• Otherwise, both numbers are odd, and a normal Euclidean step is performed, replacing the larger of the two numbers with their difference.

Binary algorithms such as these are quite common; see for instance the ipow and expm functions in the Standard Prelude. They are especially useful for hardware implementations in the microcode of the processor, since bit-shifting operations are extremely fast.

Your task is to write functions that perform the three operations given above, simulating how they are performed in hardware. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

[ Today’s exercise was written by guest author Bill Cruise, who blogs at Bill the Lizard, where he is currently describing his adventures studying SICP. Feel free to contact me if you have an idea for an exercise, or if you wish to contribute your own exercise. ]

We calculated the value of pi, and logarithms to seven digits, in two previous exercises. We continue that thread in the current exercise with a function that calculates sines.

Sines were discovered by the Indian astronomer Aryabhata in the sixth century, further developed by the Persian mathematician Muhammad ibn Mūsā al-Khwārizmī (from whose name derives our modern word algorithm) in the ninth century. Sines were studied by European mathematicians Leibniz and Euler in the seventeenth and eighteenth centuries. It was Euler who coined the word “sine”, based on an earlier mis-translation (to the Latin “sinus”) of the word “jya” used by Aryabhata. The sine of an angle is the ratio of the length of the opposite side to the length of the hypotenuse in a right triangle. (You may remember the mnemonic SOHCAHTOA if you’ve ever taken a course in trigonometry.)

One way to calculate the sine of an angle expressed in radians is by summing terms of the Taylor series:

Another method of computing the sine comes from the triple-angle formula . Since the limit , a recursion that drives x to zero can calculate the sine of x.

Your task is to write two functions to calculate the sine of an angle, one based on the Taylor series and the other based on the recursive formula. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

This exercise marks the one-hundredth exercise since the beginning of Programming Praxis. Thanks to all my readers, who give me the energy to continue writing these exercises. To celebrate, we ought to have a party, and parties need games, so in this exercise we implement a program to play the game of nim.

Nim is played by placing stones in piles; there may be any number of piles, and any number of stones in each pile. Two players alternate turns in which each player removes as many stones as desired from any pile desired; all the stones removed must come from the same pile, it is not permitted to take stones from multiple piles. The player who removes the last stone from the last non-empty pile is the winner. (In some variants of the game the player who removes the last stone is the loser, but we will not consider those variants.)

The winning strategy was devised by Charles Bouton in 1901 and revolves around a mathematical calculation called the nim-sum, which is the exclusive-or of the number of stones in each pile; for instance, with piles of 25, 32, 19, 4, and 17 stones, and representing the exclusive-or operation as ⊕, the nim-sum is calculated as 25 ⊕ 32 ⊕ 19 ⊕ 4 ⊕ 17 = 63. If a player can remove stones from a pile in such a way that he reduces the nim-sum to zero, he can always force a win. Such a move can be computed by exclusive-or’ing the nim-sum with the number of stones in each pile in turn, choosing a pile for which the result of that exclusive-or is less than the number of stones in the pile, and removing from that pile sufficient stones to reduce the pile to the result of that exclusive-or; for instance, removing one stone from the second pile in the example above reduces the nim-sum of the piles to zero, forcing a win. However, if the nim-sum is already zero at the beginning of a player’s turn, he will lose against proper play, so any move may be chosen at random.

Your task is to write a program that plays nim against a human player. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

If you like to program math puzzles, you probably know about Project Euler; in fact, it was that site that inspired Programming Praxis. We don’t usually do math puzzles here, because Project Euler does them better, but the occasional math puzzle can be entertaining, and lead to some fun programming. Hence, today’s exercise.

Your task is to write a function that, given a number n, finds all pairs of numbers x and y, with x ≥ y ≥ 0, such that x² + y² = n; for instance, 50 = 7² + 1² = 5² + 5², 48612265 = 6972² + 59² = 6971² + 132² = …, and 999 has no solutions. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Happy New Year! And best wishes for a healthy and prosperous year for all my readers and their families.

On the first day of a new year it seems appropriate to write an exercise based on calendars. Unix provides a cal command to print calendars. There are several forms:

• cal — prints a calendar for the current month
• cal year — prints a twelve-month calendar for the specified year; note that year 10 occurred during the time of Christ, so you must specify 2010 for the current year
• cal month year — prints a calendar for the specified month and year; the month is given as a number from 1 to 12
• cal -3 — prints a three-month calendar for the prior month, current month and next month

The current date is highlighted wherever it appears.

Your task is to implement cal. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2