Recursion – Fortran finally gets it!

Fortran gave recursion a wide berth until Fortran 90. More than 25 years after ALGOL first allowed recursion, Fortran was finally allowing. Modern Fortran is one of the few languages where the function must be explicitly defined recursive. The prefix recursive must be used before the keyword function.

recursive function fibonacci (n) result (fib_result)
   integer, intent (in) :: n
   integer :: fib_result
   if (n <= 2) then
      fib_result = 1
   else
      fib_result = fibonacci (n-1) + fibonacci (n-2)
   end if
end function fibonacci

In Fortran 2018, the specs changed again. Here procedures without an explicit recursive attribute behave as if recursive is specified. F2018 provides a new attribute non_recursive, which can be used to mark a procedure that may not be called recursively.

Fibonacci and pineapples

The Fibonacci sequence is found in the strangest of places. Take pineapples for example. The number of spiral rows of fruitlets (eyes) in pineapples was study as early as 1933 in an article by Linford [1] published in Pineapple Quarterly, however no reference was made to Fibonacci numbers. In a follow-up study by Onderdonk [2] in 1970 it was found that the majority of pineapples had 8-13-21 rows of fruitlets, with a few smaller ones at 5-8-13. It was suggested that a pineapple with more fruitlets for a given size would have a finer texture, and it was hoped that a pineapple with 13-21-34 rows could be found, however Onderdonk never found any pineapples exhibiting such a pattern. In general, pineapples have three series of spirals, derived from the roughly hexagonal pattern of its fruitlets, or scales.

Here is an example of the hexagonal scale patterns found on a pineapple.

[1] Linford, M.B., “Fruit quality studies II. Eye number and eye weight”, Pineapple Quarterly, 3, pp.185-188 (1933)
[2] Onderdonk, P.B., “Pineapples and Fibonacci numbers”, The Fibonacci Quarterly, 8(5), pp.507-508 (1970)

Why calculating Fibonacci using binary recursion is just plain moronic

The Fibonacci sequence is a series of numbers of the form:

0, 1, 1, 2, 3, 5, 8, 13, 21, 34,…

It is derived such that each number is the sum of the two preceding ones, starting from F(0)=0 and F(1)=1. Or in equation form: F(n) = F(n-1) + F(n-2).

The original formula for the Fibonacci series lends itself to a natural, if somewhat naïve, example of binary recursion, which is probably the most notorious implementation found in the literature. The algorithm works by returning one if n = 1 or 2, and the sum of  f(n-1)  and f(n-2) if n>2 . Represented as a recursive function in C, this looks like:

int fib_binary(int n) {
   if (n == 0)
      return 0;
   else if (n == 1)
      return 1;
   else
      return fib_binary(n-1) + fib_binary(n-2);
}

This algorithm certainly generates the correct answer, but what is its computational cost, i.e. how long does it take?  If n is less than two, the function halts almost immediately. However for larger values of n, there are two recursive calls of the algorithm. This implies that the running time of the algorithm grows at exponential time, or f(n)=2^0.694n, which for n=200 is 2^140 operations. 

To understand how this is possible, consider the Fibonacci call tree for calculating the 5th Fibonacci number, which just happens to be 5. Notice how many times F₁ is called?

fibonacci_tree

The biggest problem with using binary recursion to calculate the Fibonacci numbers is the time spent re-calculating already calculated Fibonacci numbers. For example, when calculating f(40) using recursion, f(39) is calculated once,  f(35) is calculated 8 times, f(1) is calculated 102,334,155 times, and  f(0) is calculated 63,245,986 times. All up there is a total of 331,160,281 function calls. The calculation of f(40) actually takes an average of 900 milliseconds (on my MacBook Pro with a 2.9 GHz Intel Core i5).

This is an interesting analysis, rarely made in textbooks. Few textbooks discuss alternatives to binary recursion for Fibonacci. Indeed, the use of Fibonacci numbers to illustrate binary recursion is a good example of when not to use recursion.

For anyone interested in Fibonacci numbers, there is a whole journal available online, The Fibonacci Quarterly.

Generating a recursive Romanesco broccoli

Romanesco broccoli is an edible flower of the species Brassica oleracea which includes cabbage, broccoli, Brussel sprouts, cauliflower, collard greens, and numerous other “cultivars”. It was first documented in Italy (as broccolo romanesco) in the sixteenth century. In French it is known as chou Romanesco, which translates to “Romanesco cabbage”, but it is by no means a cabbage. The Germans call it Pyramidenblumenkohl, or “pyramid cauliflower”. Romanesco broccoli is what could be termed a self-similar vegetable. There are vegetables which are composed of smaller copies of themselves ad infinitum, or at least until some limit where the similarity breaks down.

Equiangular spirals (also known as logarithmic spiral, Bernoulli spiral, and logistique) describe a family of spirals. It is defined as a curve that cuts all radii vectors at a constant angle.

How do we build a visualization of this vegetable? In Processing of course! Here is the Processing code to generate the Romanesco broccoli. The first piece of code handles setting up the Processing environment, and defining constants like the Golden ratio, the number of florets. It also sets up the Processing draw() function, which executes drawRomanesco(), and takes a snapshot of the final picture.

The remainder of the code deals with the function drawRomanesco() [1].

This is a nice recursive function, which as the base case draws a circle, using the function ellipse(). The recursive case, sets up some parameters, then uses a loop to derive the florets, which internal to the loop, recursively calls itself. The depth of the Romanesco is controlled by the number of levels, which in the case of the example above is 3. The first level is the overall Romanesco with 60 florets, the second level builds each of 60 florets within those florets, and the third level builds the florets within the florets, within the florets. The resulting image of the Romanesco broccoli, is below, showing the Romanesco from above. Of course it is not perfect, because the florets on the real Romanesco are angled out from the main axis.

[1] Code derived from algorithm on http://joyofprocessing.com.

 

Building a visualization of a sunflower spiral (ii)

So, Processing seems to be the most relevant language to code the sunflower spiral in. The calculation uses the Golden Angle (137.5°), which can be calculated in radians as follows:

angle = PI * (3.0 - sqrt(5.0))

Here is the code in Processing:

public void setup() {
   size(600,600);
}

public void draw() {
   background(255);
   translate(width*0.5, height*0.5);
   stroke(0);
   float angle = PI*(3.0-sqrt(5.0)); //137.5 in radians
   float r, theta, x, y;
   int n = 250;
 
   for (int i=1; i < n; i=i+1) {
      r = i;
      theta = i * angle;
      // Convert polar to cartesian
      x = r * cos(theta);
      y = r * sin(theta);
      ellipse(x, y, 5, 5);
   }
   noLoop();
}

This code draws 250 points in the spiral, performing the calculations in polar coordinates, and then converting them before drawing them using the ellipse() function, using a 5×5 circle. Instead of calculating the radius as r=c√i, we chose to use r=i, as this provides for a better visualization (technically an Archimedean spiral, r(k) = ck, with c=1). Now this just provides a representation of what a sunflower spiral would look like.

The thing I especially like about Processing is that by adding this line of code at the end (after noLoop()), I can save the visualization to file without any hassle:

saveFrame("sunflower250.png");

Here is the spiral:

If we modify the radius calculation to r = 5*sqrt(i), we get a much tighter configuration, and something which better represents a sunflower. The result is the canonical Vogel spiral.

Change the configuration to n=500 points gives:

It is of course easy to experiment with different spirals, just by changing the equation for r. For example  r = pow(i,0.8), is a form of parabolic spiral.

 

 

 

 

Building a visualization of a sunflower spiral (i)

Sometimes you want to build a picture to illustrate something, but not necessarily do it by hand. The answer is to write an algorithm to do it, and code it in some visualization happy language… and there aren’t that many really. But using Processing works quite well. Imagine if we wanted to draw a sunflower spiral. What is a sunflower spiral you may ask? Okay, some background first.

Sunflower heads are interesting because of Fibonacci sequences. Fibonacci sequences appear quite commonly in nature, as two consecutive Fibonacci numbers, especially in plant life. Examples include the arrangement of leaves on a stem, the fruitlets of a pineapple, scales on a pine cone, artichoke flowers, fern fronds, or plant seed heads. This is an aspect of plant form known as phyllotaxis (the botanical study of the arrangement of phylla (leaves, petals, seeds, etc.) on plants)The seeds on a seed head are often distributed in the head in two distinct sets of spirals which radiate from the centre of the head to the outermost edge in clockwise and counterclockwise directions. The spirals are logarithmic in character. The number set exhibited by the double set of spirals is intimately bound up with Fibonacci numbers. It appears that the reason for this formation is to allow the seedheads to pack the maximum number of seed in the given area. The most perfect example of phyllotaxis are afforded by the common sunflower (Helianthus annuus,L.) [1]. The head of the sunflower which is approximately 3-5 inches across the disk will have exactly 34 long and 55 short curves. A larger head 5-6 inches in diameter will show 55 long curves crossing 89 shorter ones.

Fig 1: Sunflower head showing seeds [2]

Consider the image shown in Fig. 1, depicting a sunflower head. The shallower spiral of this sunflower repeats 34 times, whilst the steeper spiral repeats 55 times.

So to draw this we have to produce a spiral, which is sometimes called a Fibonacci spiral. A very simple model of the patterns in a sunflower seed head was developed by Helmut Vogel in 1979 [3]. He calculated the position of the seeds (k) in the head using polar coordinates (r, δ):

r(k) = c√k

θ(k) = kδ

where δ is the golden angle, 137.5°, and c is a constant scaling factor (some algorithms replace c√k with k). As the value of k increases, the position rotates through the golden angle, and the radius increases as the square root of k. This allows use to draw a sunflower head. This spiral is essentially identical to the orbit of a “nonrelativistic charged particle in a cyclotron” [3].

[1] Sutton, C, “Sunflower Spirals Obey Laws of Mathematics.” New Scientist , 18 April 1992, p. 16.

[2] Church, A.H., On the Relation of Phyllotaxis to Mechanical Laws, Williams and Norgate, London.Plate VII (1904)

[3] Vogel, H., “A better way to construct the sunflower head”, Mathematical Biosciences 44, pp.179–189 (1979).

 

Fibonacci in Swift

So, I didn’t really need another distraction, but Swift has an interactive environment, so now I’m playing. Here are two functions for Fibonacci, first an iterative one, and then a recursive one.

Here’s the iterative version:

func fibonacciI(n: Int) {
    var f1=1, f2=1, fib=0
    for i in 3...n {
        fib = f1 + f2
        print("Fibonacci: \(i) = \(fib)")
        f1 = f2
        f2 = fib
    }
}

And now the recursive version:

func fibonacciR(n: Int) -> Int {
    if (n == 0){
        return 0
    } else if (n == 1) {
        return 1
    }
    return fibonacciR(n-1) + fibonacciR(n-2)
}

And now the version using arrays:

func fibonacciA(n: Int) {
    var fib: [Int] = []
    fib.append(1)
    fib.append(1)
    for i in 2..<n {
        fib.append(fib[i-1]+fib[i-2])
    }
    for i in 0..<n {
        print("\(fib[i])")
    }
}

Nothing terribly different about these implementations from other languages.  I do like the way it is possible to create an array and then append items to it. The range constraints in loops, i.e. and ..< are kind of interesting as well.

 

 

Timing an Ada function

How does one time a function in Ada?

Below is a simple program which tests the run time of a recursive Fibonacci function. It uses the ada package Ada.Calendar, which is able to return time stamps. The variables startTime, and endTime are capable of holding a Time value. A call to Clock(), reads the machine internal clock, and stores the current time of day in the time variable. The difference is stored in a variable of type Duration, expressed in seconds.

Note also the put_line statements which make use of the Image criteria – X’Image(Y) is an Ada attribute where X is any discrete type and Y is an instance of that type. This attribute returns a string representation of the value passed as input.

And if you want to time another language check out the Wikipedia site on system timing.

with Ada.Text_IO; use Ada.Text_IO;
with ada.Integer_Text_IO; use Ada.Integer_Text_IO;
with Ada.Calendar; use Ada.Calendar;
with Ada.Text_Io; use Ada.Text_Io;

procedure testFibonacci is

    function fibonacci(n : in integer) return integer is
    begin 
        if n = 0 then
            return 0;
        elsif n = 1 then
            return 1;
        else
            return fibonacci(n-1) + fibonacci(n-2);
        end if;
    end fibonacci;

    n : integer;
    startTime, endTime : Time;
    milliS : Duration;
    result : integer;

begin
    put_line("Fibonacci number: ");
    get(n);
    startTime := Clock;
    result := fibonacci(n);
    endTime := Clock;
    milliS := (endTime - startTime) * 1000;

    put_line("Result: " & Integer'Image(result));
    put_line("Runtime = " & Duration'Image(milliS) & " milliseconds.");
 
end testFibonacci;