Calculating Babylonian square roots (ii)

How is the Babylonian square roots algorithm implemented? Let’s try it in Fortran:

program squareroot
   implicit none

   real :: epsilon = 0.0000001
   real :: n, y, r

   write(*,*) "SQRT of what number? "
   read(*,*) n

   ! Make an initial guess
   r = 1.0
   ! Divide the original number by the guess
   y = n / r
   ! Iterate until accuracy is achieved
   do while (abs(r-y) > epsilon)
      r = (r + n/r) / 2.0
      y = n / r
   end do

   write(*,10) "sqrt(",n,") = ", r
   10 format(a5,f5.2,a4,f10.8)

end program squareroot

Here is the code running:

SQRT of what number?
2
sqrt( 2.00) = 1.41421354

We could also print out the value of r in the loop:

sqrt( 2.00) = 1.50000000
sqrt( 2.00) = 1.41666675
sqrt( 2.00) = 1.41421568
sqrt( 2.00) = 1.41421354

Here is the same code in Ada:

with ada.Text_IO; use Ada.Text_IO;
with ada.Float_Text_IO; use Ada.Float_Text_IO;
with ada.Numerics.Elementary_Functions;
use ada.Numerics.Elementary_Functions;

procedure babylonianSQRT is
   epsilon, n, y, r : float;

begin
   put_line("SQRT of what number? ");
   get(n);
   epsilon := 0.0000001;

   -- Make an initial guess
   r := 1.0;
   -- Divide the original number by the guess
   y := n / r;

   -- Iterate until accuracy is achieved
   while abs(r-y) > epsilon loop
      r := (r + n/r) / 2.0;
      y := n / r;
   end loop;

   put("sqrt(");
   put(n, 2, 3, 0);
   put(") = ");
   put(r, 2, 8, 0);

end babylonianSQRT;

 

 

 

 

Calculating Babylonian square roots (i)

There is an ancient algorithm which can be used to approximate the square root of a number. It is found written in cuneiform (wedge-shaped script) on YBC 7289, one of the best known of the Babylonian mathematical clay tablets. It originates from the region of southern Mesopotamia (present day Iraq), and is likely from the first third of the second millennium BC [1]. The Babylonians worked in sexagesimal, or base 60 (the vestiges of this system are seen in our method of system of time). They also worked in floating-point sexagesimal numbers. On the tablet, the square root of two is calculated to three sexagesimal places. The square-root is written as:

babylonianSQRT

In modern characters, this would be written approximately as 1 ; 24, 51, 10, where a semicolon is used to separate the integer and fractional parts, and the comma is used as a separatrix for the sexagesimal positions. This Babylonian calculation of √2 is performed as 1 + 24/60 + 51/60² + 10/60³, which is approximately 1.414212963.

It is defined by C.B. Boyer in his “History of Mathematics” [2] as:6

“Let x = √a be the root desired, and let a₁ be a first approximation to this root; let a second approximation b₁ be found from the equation b₁=a/a₁. If a₁ is too small, then b₁ is too large , and vice versa. Hence the arithmetic mean a₂=½ (a₁ + b₁) is a plausible next approximation. Inasmuch as a₂ always is too large, the next approximation b₂=a/a₂ will be too small, and one takes the arithmetic mean a₃=½ (a₂ + b₂) to obtain a still better result; the procedure can be continued indefinitely.”

The technique is equivalent to Heron’s method for extracting a square root, found in Metrica I 8 [4]. Here is an example for calculating √2, with N=2, and r₀=1 as the initial guess (blue digits represent those accurately calculated).

r₁ = (r₀ + N/r₀)/2 = (1 + 2/1)/2 = 1.5
r₂ = (r₁ + N/r₁)/2 = (1.5 + 2/1.5)/2 = 1.416666
r₃ = (r₂ + N/r₂)/2 = (1.416666 + 2/1.416666)/2 = 1.4142156
r₄ = (r₃ + N/r₃)/2 = (1.4142156 + 2/1.4142156)/2 = 1.4142135

 

[1] Fowler, D., Robson, E., “Square root approximations in old Babylonian Mathematics: YBC 7289 in context”, Historia Mathematica, 25, pp.366-378 (1998)
[2] Boyer, C.B., History of Mathematics, John Wiley & Sons (1968).
[3] Knuth, D.E., “Ancient Babylonian algorithms”, Communications of the ACM, 15(7), pp.671-677 (1972)
[4] Smyly, J.G., “Square roots in Heron of Alexandria”, Hermathena, 63, pp.18-26 (1944)

Loading language syntax into the brain is expensive

The basic cognitive architecture of the brain is composed of sensory memory, working-memory and long-term memory. Information enters the brain through the sensory memory is consciously processed through the working memory, and is stored long-term in the long-term memory, or at least that is the ideal process: long-term storage. The sensory memory has a latency of several hundred milliseconds, the short term memory a decay period anywhere from 15-30 seconds. How we store store things in the long-term memory is different for everyone, and truthfully not fully understood. The entire process of loading content into the brain is expensive, because it takes time, and often the experience of loading that information may not be complete.

The process of loading programming content into the brain is one of loading syntactic structures, and how they can be used to implement a particular algorithm. This can be problematic if the structure in question is somewhat ambiguous. Take for example, the following if statement in C:

if (x < 100)
   y = calculate_Vortex();

It seems pretty straightforward, and easy to remember. Basically it says if x is less than 100, perform calculate_Vortex(). Nothing untoward there, in fact the concept of a generic if statement is easy for a novice programmer to learn, largely because they encounter decisions in everyday life. It is the syntax that becomes problematic. If the novice programmer were now to add a second task when x<100, they might create a piece of code like this:

if (x < 100)
   y = calculate_Vortex();
   x = y * 4.87;

But this code in C is incorrect because it is missing the block encapsulators, and effectively does not associate the statement “x = y *4.87;” with the condition when x<100. A { } pair is added, which transform the code into this:

if (x < 100) {
   y = calculate_Vortex();
   x = y * 4.87;
}

Effectively there are a number of differing syntactic variations when using an if statement in C. The novice programmer then has to modify the learned behaviour of an if statement. This is different to the syntax of if statements in other languages, where there is no syntactic variation, and the retention of memories is therefore less expensive. Here is the same statement in Julia.

if x < 100
   y = calculate_Vortex()
   x = y * 4.87
end

In Julia,  the same if structure exists, regardless of how many statements it encapsulates. The same cognitive issue occurs with for loops in C. Any programming language whose structures force a novice programmer to doubt how they are used causes an increase in cognitive load..

Removing goto’s that jump backwards

One of the issues when re-engineering an old Fortran program is the global goto’s that jump backwards. They kind-of emulate a loop, but sometimes cross over. Here is an example:

jumpBackwards

In this situation, if the conditional, C is true, then the program goto’s to label 20, otherwise it goto’s to label 10. In both situations, the code in B is executed, but only if C is false is the code in A executed. This code could look something like this:

   program gotoWITHloop
   integer :: turn
   turn = 1
10 write(*,*) "code A"
   ! This part of the code only happens when turn is 1
20 write(*,*) "code B"
   ! This part of the code only happens when turn is both 1 and 2

   if (turn .eq. 1) then
      turn = 2
      goto 20
   else
      turn = 1
      goto 10
   end if
   end

The conditional is related to whether the variable turn is 1 or 2. First time into the loop both code A and B are activated, and then the conditional is applied. If turn is equal to 1, the program jumps to label 20, then back to the conditional etc. If turn is equal to 2, then the program jumps to label 10 at the top, and flows through etc.

How can this be fixed, to remove the goto statements, and their associated labels?

One way of re-engineering this involves the use of a loop, and a subroutine.

program gotoWITHloop
   integer :: turn
   turn = 1
   do
      if (turn .eq. 1) then
         call codaA()
      end if
      write(*,*) "code B"
      ! This part of the code only happens when turn is both 1 and 2
      if (turn .eq. 1) then
         turn = 2
      else
         turn = 1
      end if
   end do
end
! This part of the code only happens when turn is 1
subroutine codeA()
   implicit none
   write(*,*) "code A"
end subroutine codeA

As the stuff in code A is only activated with turn equal to 1, it can be situated in an if statement. It also might be nice, to encapsulate the code itself into a subroutine. The labels can then be removed, and the entire code situated inside a loop.

This of course is just one way to deal with this sort of a situation. In many cases removing these goto’s involves a combination of loops, if statements, and subroutines.

 

Random integers in a range

Amazingly, sometimes the internet has very little proper information on a topic. In the context of generating random numbers, an algorithm should distribute numbers evenly across a range of n numbers such that the probability of generating a number is 1/n. This is called a uniform distribution. Now this is fine, but there may be times when we want to generate a number within a certain specified range, say between 0 and p. For example, to get a number between 0 and 9, we could use the following pseudocode (assuming randomNum is an integer of some size > 0):

d = mod(randomNum,10)

To make the range between 1 and 10:

d = mod(randomNum,10) + 1

Any positive number divided by 10 will give a remainder of between 0 and 9, adding 1 gives a range of between 1 and 10. If you want a range [p..q], try:

d = p + mod(randomNum,(q-p+1))

Here is a C program that generates 1000 random numbers in the range 12 to 73:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

int main(void)
{
   int i, p=12, q=73, r;
   // Seed the random number generator
   srand((unsigned int)time(0));

   // Generate 1000 random numbers between p and q
   for (i=1; i<=1000; i=i+1){
      r = p + (rand() % (q-p+1));
      printf("%d ", r);
   }
   return 0;
}

 

 

 

What is throw-away code?

A piece of code created to prototype something should be discarded when the concept has been proven – and rewritten properly. This code has several names: throwaway-code, a quick hack, kleenex code, or disposable code. Why do programmers write throwaway code? Firstly as a proof of concept. A good example is writing simple Python scripts to process images – it’s easy to test the outcome of ideas by writing nice scripts. Secondly, because it doesn’t take a lot of time – there may be no time to write an elegant, well designed and tested program. Thirdly, it may be written as an alternative to using someone else’s code. Throwaway code is often messy – it may lack structure, lack readability, have little or no documentation, or be highly inefficient.

Unfortunately it often becomes permanent. Throwaway code it meant for prototyping – it is not meant to become the implementation.

 

[1] Foote, B., Yoder, J., “Big ball of mud”, Conf. on Patterns Languages of Programs (1997).

7 Degrees of Correctness

Novice programmers often believe that a program is either right or wrong. This comes from a false sense of security provided by the compiler. If it compiles it must be correct. Not so. Scowen and Ciechaniwicz [1]  pose seven degrees of correctness for a program. We will look at these seven degrees, using a program to calculate factorials as an example.

1. Nonsense

These programs are not syntactically correct. The compiler is usually adept enough to find these problems.

1  #include <stdio.h>
2
3  int main(void)
4  {
5   
6     double fact
7     for (i=2; i<100; i=i+1)
8        fact = fact * i;
9     printf(“Factorial = %d”, fact);
10    return 0;
11 }

This program contains a whole series of syntax errors:

  • Line 5: No declaration for the variable i, used in line 7.
  • Line 6: missing semicolon
  • Line 9: %d used to print a double.

The program should look like:

1  #include <stdio.h>
2
3  int main(void)
4  {
5     int i;
6     double fact;
7     for (i=2; i<100; i=i+1)
8        fact = fact * i;
9     printf(“Factorial = %f”, fact);
10    return 0;
11 }

2. Illegal

These programs are syntactically correct but their actions are undefined at runtime. For example:

1  #include <stdio.h>
2
3  int main(void)
4  {
5     int i;
6     double fact;
7     for (i=2; i<100; i=i+1)
8        fact = fact * i;
9     printf(“Factorial = %f”, fact);
10    return 0;
11 }

When this program runs, we’ll get an answer, but it won’t be correct, because the value of fact was never initialized. Line 6 in the program should read:

6     double fact=1.0;

Another good example is the “array out of bounds” error that occurs in languages such as C.

3. Almost Working

The program gives the correct answer sometimes, not always. Such programs are sometimes described as 90% complete. For example:

1  #include <stdio.h>
2
3  int main(void)
4  {
5     int i, n=5;
6     double sum=0.0, a[5]={1,2,3,4,5};
7     for (i=0; i<n; i=i+1)
8        sum = sum + a[0];
9     printf(“Sum = %f”, sum);
10    return 0;
11 }

This program really only works when the value of n=0, such that the value of a[0] is summed. If n has any other value, the result is always a multiple of a[0], as there is no variable index used.

4. Fragile

These programs are syntactically correct and output the correct answer with valid input data. But illogical data will result in the program going nut-so. Consider the following program which divides one number by another.

1  #include <stdio.h>
2
3  int main(void)
4  {
5     double num, den, div;
6     scanf(“%lf %lf”, &num, &den);
7     div = num / den;
8     printf(“div = %f”, div);
9     return 0;
10 }

This program works when the denominator (den) does not equal zero. If it does equal zero, then a divide-by-zero error occurs.

5. Suspicious

Programs that are syntactically and semantically correct, but contain operations which have no logical effect. The program contains operations for no apparent reason. Many modern compilers check for used variables, but don’t detect variables which are assigned a value, but never used beyond that. Consider the following example of a function which calculates the area of a circle, given the radius:

double area(double radius)
{
   double circ;
   circ = 2.0 * pi * radius;
   return pi * radius * radius;
}

The calculation of the circumference (circ) in the function is superfluous – it is never used beyond the assignment.

6. Wrong

A program may be correct, robust and never fail, yet it might be wrong because the client asked for a different function to be calculated. For example calculating a mean instead of a median, the measure in the specifications. Consider the case of calculating arithmetic versus geometric means in investment returns. If the returns for five years were 12%, 9%, 27%, 5% and -16%, the arithmetic mean would be:

(12+9+27+5-16) / 5 = 47/5 = 9.4%

However annual investment returns are not independent of each other. If there is a major drop in the market one year, there is less capital to invest the next year. The geometric mean will provide a more accurate measurement of the actual average annual return. The geometric mean would be:

(1.12 × 1.09 × 1.27 × 1.05 × 0.84)1/3 = 1.065 = 6.5%

An average of 6.5% per annum increase over five years. Note that 1 is added to each investment return to avoid problems with negative returns. The use of arithmetic mean in a program displaying investment returns will give the wrong impression of actual average returns.

7. Correct

A program which is correct in every way. It is accurate, robust, and is capable of dealing with any input. In reality very few of these programs exist. Here’s one:

int main(void)
{
   return 0;
}

[1] Scowen, R.S., Ciechanowicz, Z., “Seven sorts of programs”, SIGPLAN Notices, 17(3), pp.74-79 (1982)

Coding Ada: Random numbers

Random numbers are easy in Ada. That is largely because things are intrinsic to the language. The package Ada.Numerics.Discrete_Random provides a random number generator suitable for generating random integer values from a specified range.

Firstly, a with clause is used to associate the package:

with ada.numerics.discrete_random;

Then we create a range of numbers:

type randRange is range 1..100;

Now a new instance of the package is created, using the specified range:

package Rand_Int is new ada.numerics.discrete_random(randRange);
use Rand_Int;

Then a random-value generator gen can be initialized:

gen : Generator;

It is then prepared for use by calling the procedure reset(gen). Random values can then be generated by calling the function random(gen), which will produce a new random value of type randRange from the generator gen each time it is called.

reset(gen);
num := random(gen);

Here’s a complete piece of code to do this:

with ada.Text_IO; use Ada.Text_IO;
with ada.numerics.discrete_random;

procedure randomN is
   type randRange is new Integer range 1..100;
   package Rand_Int is new ada.numerics.discrete_random(randRange);
   use Rand_Int;
   gen : Generator;
   num : randRange;
begin
   reset(gen);
   num := random(gen);
   put_line(randRange'Image(num));
end randomN;

 

 

The world of unstructured programming, i.e. spaghetti code

Much of the programming prior to 1970 was what is now considered unstructured. It is often synonymous with the pejorative term spaghetti code, meaning source code with a complex and tangled control structure. Languages such as Fortran were heavily reliant on programming structures such as goto, which provided much of their control framework. It provided jumps from one part of the program to another part – often in an extremely unstructured manner. By the mid-1960s, programming pioneers such as Edsger W. Dijkstra were starting to become concerned about programming clarity, correctness and structure. He contended that modular, goto-less programs are more easily created and understood by humans. In 1968 Dijkstra published a letter, “Go To statement considered harmful” in which he condemned the indiscriminate use of the goto statement [1]. One statement stands out: “The go to statement as it stands is just too primitive, it is too much an invitation to make a mess of one’s program”. Instead he advocated the use of procedures, and conditional and repetition for program control.

Fortran may have been one of the greatest offenders when it came to spaghetti code due to all the control structures which relied on jumping a label. Fortran contained the DO looping structure, the computed GO TO, the unconditional branch, the assigned GO TO, and the arithmetic if. One of the worst offenders may be the ability to use a goto to enter and leave a still executing do loop. Consider the small piece of code below, which was once legal (now the compiler will throw a warning, but still compile. It basically jumps from inside the DO loop to outside the DO loop, to back inside the DO loop again, ad infinitum. The value of the index variable I, actually never changes from its initial value 1.

      PROGRAM DO_LOOP
      J=0
  10  DO 15,I = 1,10
      GO TO 13
  12  J=J+1
  13  WRITE(*,*) "inside loop ", J
      GO TO 20
  15  CONTINUE
  20  J=J+1
      WRITE(*,*) "outside loop ", J
      GO TO 12
      END

Fortran didn’t introduce structured programming until 1978, gradually removing or restricting the offending jump statements. The goto statement is therefore intrinsically linked to the notion of spaghetti code. Adding one goto statement to a program does not make it spaghetti code. Spaghetti code is best understood by analyzing a legacy Fortran or Basic program. If it is impossible to follow the algorithmic flow of a program, then it is likely spaghetti code.

Consider an algorithm which takes a year and numeric value for a day in the year and calculates the day of the particular month the number represents.

doy_unstructured

Now compare the structured code to the spaghetti version. Which program is more readable, or can be easily traced?

doy_structured

Now if you have a haphazardly structured, sprawling, sloppy, duct-tape-and-baling-wire, spaghetti-code jungle it’s called a big ball of mud [2].

 

[1] Dijkstra, E., “Go To Statement Considered Harmful”, Communications of the ACM, Vol.11:3, pp.147–148, 1968.

[2] Foote, B., Yoder, J., “Big ball of mud”, Conf. on Patterns Languages of Programs (1997).

 

Usability confusion: Streetcar buttons

Everyday you can find usability confusion with everyday objects. On the new TTC Flexity articulated streetcars, each of the doors has a button like the one below. This button has two functions. The first is that it can be used to indicate that you want the streetcar to stop at the next stop. The issue is that most other buttons found on streetcars and buses use the word STOP on the request button. This means that many people don’t actually know that the button can be used in this manner. The second use is to open the door, which by all indications is the buttons primary function. However when a streetcar stops, all doors open automatically (or so it appears), so the button used in this context is not useful. The only time it is of relevance is when the streetcar remains in one position for a extended period, and the doors close. Then the button can be pushed, while the streetcar is stationary to exit the vehicle.

streetcarButton
TTC streetcar door button

The external side of the button achieves the same thing. These are fairly standardized door opening buttons, but they just don’t seem to work as they should, and it is easy to observe by watching people use it. It might be better if the outer red part of the button was press-able, maybe with the words “NEXT STOP”, and the inner button is for opening the door.