Last updates:
Mon Mar 15 17:44:42 2004
...
Thu Oct 9 14:35:24 2008
The book The Java Programming Language, third edition, Addison-Wesley, 2000, ISBN 0-201-70433-1, contains a short example on p. 7 of a Java program for computing and printing the Fibonacci sequence: 0 1 1 2 3 5 8 13 21 34 ... Each term after the first two is the sum of its two predecessors. The programs described below ignore the zero-th order term.
This series has a venerable history going back to the year 1202, when Leonardo Pisano (Leonardo of Pisa), also called Leonardo Fibonacci (Filius Bonaccii, son of Bonaccio) discovered it. The Fibonacci series turns up in a surprising number of places in discrete mathematics and computer science, including the problem of building a computer network that has maximal reliability. In Leonardo's case, its partial sum is the solution to the problem:
How many pairs of rabbits can be produced from a single pair if each pair produces a new pair every month, each new pair reproduces starting at the age of one month, and rabbits never die?
In 2002, on the 800th anniversary of its original
publication in Latin, Fibonacci's book was republished in
English; I've created a separate Web page for it,
available
here.
A comprehensive new book on Fibonacci and his number sequence was recently published: Alfred S. Posamentier and Ingmar Lehmann's The fabulous Fibonacci numbers, Prometheus Books, New York, 2007, ISBN-10 1-59102-475-7, ISBN-13 978-1-59102-475-0.
The Fibonacci sequence terms grow exponentially, asymptotically as powers of 1.618033...; that special number is known as the golden ratio, exactly equal to (1+sqrt(5))/2. There are two recent books about that number: Roger Herz-Fischler's A Mathematical History of the Golden Number, Dover, New York, 1998, ISBN 0-486-40007-7, and Mario Livio's The Golden Ratio: The Story of Phi, the World's Most Astonishing Number, Broadway Books, New York, 2002, ISBN 0-7679-0815-5. The Fibonacci terms grow surprisingly rapidly: after 23 terms, 16-bit integers are insufficient. After 46 terms, 32-bit integers are inadequate. After 52 terms, 36-bit integers are insufficient. After 78 terms, 53-bit integers overflow. After 92 terms (less than 8 years in Leonardo's rabbit problem), 64-bit integers fail. After 186 terms, 128-bit integers overflow.
Switching to floating-point arithmetic only helps a little: 32-bit IEEE 754 arithmetic handles 181 terms, 64-bit IEEE 754 handles 1471 terms, and 128-bit IEEE 754 handles 23,597 terms, the last of which has 4932 decimal digits, although only the first 36 are representable.
The Fibonacci terms provide a vivid demonstration of the inadequacy of fixed-size computer integer arithmetic systems. Worse, almost no programming languages detect integer overflow: TeX is a notable exception, but it only detects overflow from multiplication, not from addition! The Algol 60 implementation on TOPS-20 also caught the overflow.
Several of the programs below exploit this failure-to-detect-integer-overflow deficiency by looping until a negative Fibonacci term is encountered. How does the sum of two large positive numbers become negative? The addition produces an overflow bit that is stored in the high-order bit position, which is the sign bit on all modern computer arithmetic systems. For example, in a 3-bit twos-complement arithmetic system 011(base 2) + 011(base 2) == 110(base 2) == -2(base 10).
For more on the interesting history and applications of Fibonacci numbers and Fibonacci sequences, see Donald Knuth's monumental multivolume work, The Art of Computer Programming, Addison-Wesley (1968--date).
For comparison with other programming languages, I (with contributions from two others, credited below) implemented an extension of this code in 49 languages. I then extended the Java program to print all terms that are representable in 64-bit integer arithmetic, and to format the output into a neat aligned table. This takes a surprising amount of extra code, and the resulting Java program is three times the length of Fortran and C equivalents. In fairness, the Java version can produce internationalized output, and offers digit grouping for improved reability. Of the other languages sampled below, only Common Lisp and guile (GNU Scheme) provide digit grouping, but it is possible to add that feature to awk and perl with only a modest effort (a 10-line function). Three versions of the Java program are provided below, one with default formatting, and two with fancy formatting.
In the examples below, some languages impose severe restrictions on the size of integers (Bigloo Scheme: 29 bits, Emacs Lisp: 28 bits on 32-bit systems, 32 bits on 64-bit systems [the Emacs Lisp manual says that integers have a minimum of 29 bits, but that is easily proved wrong by experiment], Metafont: 12 bits, PostScript: 32 bits, TeX: 16 bits). Some (axiom, bc, Common Lisp, maple, mathematica, Objective Caml (ocaml), reduce, Scheme (but not the Bigloo implementation)) support arbitrary-precision integers. Others (awk, hoc, JavaScript, matlab, perl, PHP, python, ruby) masquerade integers as floating-point values, which imposes a 53-bit limit in IEEE 754 arithmetic. Ada and Java support 64-bit integers natively, while only some C, C++, and Fortran compilers have been extended with support for such values. TeX has no floating-point arithmetic at all. Instead, it has fixed-point dimensions, but not fixed-point numbers, so columns 3 and 4 of the output are suffixed by pt (printer's points).
Arbitrary-precision floating-point arithmetic is standardly provided by axiom, bc, Common Lisp, maple, mathematica, python, reduce, and Rexx, and there are separate BigFloat packages available for perl and ruby.
A Web search turns up several multiple-precision arithmetic packages for C, C++, and Fortran, but none of them are standardly integrated into those languages.
Java also has two standard support packages, BigInteger and BigDecimal, for arbitrary precision integer and fixed-decimal arithmetic. Unfortunately, the latter contains no elementary functions, and since it does not provide a floating exponent, it is useful primarily for currency calculations. One of the test programs below exhibits use of these packages.
Unlike Ada and C++, Java lacks operator overloading, so big-number arithmetic operations must be done through function calls, and because those functions are class methods, the name order is unfamiliar to people from a more traditional programming-language background. The simple expression
ratio = hi/lo
must be rewritten for big-number arithmetic as
BigDecimal ratio = new BigDecimal(hi);
BigDecimal den = new BigDecimal(lo);
ratio = ratio.divide(den, Fractional_Digits,
ratio.ROUND_HALF_DOWN);
This requires two object creations (and subsequent garbage collection), specification of precision, and rounding.
Obviously, if such operations were required frequently, one would define a function to hide these details, allowing the user to write instead something like this:
BigDecimal ratio = bd_div_bi(hi,lo);
Although this is more compact, object creation and garbage collection are still required.
Java does not (yet) have a BigFloat package: at the Sun Java Web site, http://java.sun.com/products/jdk/1.2/docs/guide/idl/mapping/idl.fm.html, I found this statement:
5.4.9 Future Long Double Types
There is no current support in Java for the IDL long double type. It is not clear at this point whether and when this type will be added either as a primitive type, or as a new package in java.math.*, possibly as java.math.BigFloat.
This is left for a future revision.
On 10-Jun-2002, I learned about an effort spearheaded by
Mike Cowlishaw at IBM UK to extend Java's BigDecimal page
to include decimal (as opposed to binary) floating-point
arithmetic; see
http://www.jcp.org/jsr/detail/13.jsp
for details. If this becomes standardly available, it
will certainly be a welcome addition to Java!
There are standard big-number packages, Integers and Reals, for Oberon-2, illustrated in the second Oberon-2 example below. Unfortunately, Oberon-2 appears to have no standard provision for padding of output strings, or conversion of arbitrary-precision numbers to fixed (instead of floating) decimal forms, so the output does not line up nicely. There is also unwanted intermediate expression crud output before the program's own output.
The Oberon-2 Optimizing Compiler Reference Manual describes a 64-bit integer data type, HUGEINT, said to be available only in the 64-bit implementations of the system. I downloaded that distribution and unpacked it, but a search of the library code shows that there is (as yet) no library support whatever for HUGEINT.
Several of these languages (axiom, bliss, JavaScript, mathematica, matlab, Metafont, PostScript, reduce, S-Plus, TeX) lack formatted I/O statements, making it nearly impossible to produce neat tabular output. Many of them assume that there is only one possible decimal representation of a floating-point number, and that `all' digits should be shown: this is simply nonsensical!
The Fibonacci problem is simple enough that implementation in most programming languages is straighforward. However, in a few languages, unexpected issues arose. They are documented in the following subsections.
The axiom program was rather frustrating to write, because the interpreter complains about a missing mate when the code was written with one statement per line. I ultimately reduced the loop body to a single line to get the code accepted. The line count in the table below counts statements, rather than lines, to give a fair count. Also, I found that axiom insists that its input file names end in .input, even though the axiom source tree typically uses the extension .as, and input files cannot be supplied on the command line. There is probably some obvious way around this nonsense, but I could not figure out how, despite having the axiom book at hand, and looking at lots of examples in the axiom source tree.
Digital Equipment Corporation (DEC), and its later owners Compaq and Hewlett-Packard, used the Bliss system programming language on PDP-10, PDP-11, VAX, Alpha, and IA-64 systems, for multiple operating systems on some of these architectures. Bliss was developed at Carnegie-Mellon University about 1970 (see its Wikipedia article ), and received wide use inside DEC, mostly for utilities and language compilers.
Bliss is a typeless language: an identifier in the language means the address of a storage location, rather than its value. To get the value, the dereferencing operator prefix, dot, must be used. Thus, foo is the address, and .foo is the value at that address; this difference is a frequent cause of programming errors for Bliss novices. Because all objects are words in the underlying architecture, the size of an object is too often visible, and causes portability issues: foo + 6 means 6 words past foo on the PDP-10, but 6 bytes past foo on byte-addressed architectures. However, the language provides functions to help conceal wordsize and addressing differences.
DEC documentation refers to Common Bliss as a machine-independent language, with variants Bliss-11 (PDP-11), Bliss-32 (VAX), and Bliss-36 (PDP-10); Bliss-64 (Alpha and IA-64) was added later.
The language lacks I/O facilities, with the expectation that programmers would use suitable system-call interfaces for such access. DEC's tutorial I/O package, tutio.r36, is a simple example that works on TOPS-10 and TOPS-20, but not on other operating systems. It supports character, integer, and string I/O, but not floating-point.
Regrettably, DEC never made Bliss compilers a standard component of their operating systems; they were only available at additional licensing cost, and as a result, most DEC customers never had them. This, plus the lack of a standard run-time library, meant that little software development was done in Bliss outside DEC and CMU.
For the purposes of the Fibonacci problem, it was necessary to write a small floating-point output routine that calls the TOPS-20 Monitor (kernel in Unixese) via the JSYS() function, which takes a system-call skip count, a system-call number, and up to four arguments in registers 1 through 4.
Because Bliss operates on words, multiword values, such as a double-precision number, are awkward to deal with. Numerical constants can be represented by literals of the form %F'3.14159' (single precision) and %D'3.14159' (double precision), but the manuals contained no examples of how to work with multiword literals, and it took considerable experimentation to find a solution.
Because the language is typeless, all explicit arithmetic is fullword integer arithmetic. Thus, floating-point arithmetic requires descent into machine language, for which there are fortunately standard builtin functions. For subtraction and division, the argument operand order is the reverse of what most programmers of other languages would expect.
Bliss contains many powerful features, and its closeness to the underlying architecture allows compilers to generate efficient code. Interestingly, Bliss has a VOLATILE attribute that predates the introduction of the volatile type modifier in C89 by about two decades.
The C++ cout statement normally outputs values in a field of data-dependent minimal width. However, C++ has I/O manipulators that can be used to set field widths and request fixed or scientific floating-point form. Unfortunately, none of the 35+ C++ compilers available to me supports the 1998 ISO C++ Standard. I found that only one compiler (Sun Solaris 2.8 CC) recognized the fixed manipulator by default, and because there is considerable variation in header file names, it required some experimentation to find a set of #include statements that worked with all of those C++ compilers that I tested.
The C# language resembles Java sufficiently that it was straightforward to prepare a C# version of the Fibonacci program starting from the Java version, and changing all of the library calls. Like Java, C# has only limited IEEE 754 support, with only 32-bit float and 64-bit double types, both with support for subnormals, only one rounding mode (the IEEE 754 default of round to nearest even), only one kind of NaN, no signed zero, and a math library comparable to that or Fortran 77 or C89, Unlike Java, C# has a reasonably powerful formatted I/O capability, albeit somewhat different from that of other languages, and one of the format types provides for digit grouping, making it easy to produce a very readable table of Fibonacci numbers.
C# also has a decimal pseudo-floating-point type intended for currency calculations, with constants suffixed M, possibly for money. C# uses a 128-bit representation of a scaled decimal value offering about 28 significant digits, and a limited exponent range of 1.0e-28 to 7.9e+28. Because the storage format is binary, rather than decimal, the overflow limit is not the expected 9.9999...e+28. Small values underflow silently to zero, but large values overflow, throwing an exception that must be caught in a try {} catch {} block if premature termination is to be prevented.
There is no math library for the decimal type, but it is simple to make a decimal square root routine that uses a call to the double Math.Sqrt() routine to get a starting value; one or two Newton-Raphson iteration steps then produce an accurate result. The range of the decimal type suffices to compute the first 139 Fibonacci numbers.
The JavaScript program was not difficult to write: it is a curious mix of HTML, awk-like arithmetic, and Java-like syntax.
I do not have any standalone JavaScript interpreters (e.g., Netscape's JSRef, or Microsoft's JScript), but the program can be run in Java-enabled browsers, such as Firefox, Konqueror, Microsoft Internet Explorer, Mozilla, Netscape Navigator, or Opera.
Unfortunately, capturing the output is problematic: you can print the window contents, but if you save the window with the displayed output, you get just the original JavaScript program.
If you use cut-and-paste with Netscape Navigator, each cell is on a separate line, and the table has to be tediously edited back into shape (I did that for the output below).
With the Microsoft browser, cut-and-paste on Sun Solaris produces complete garbage: long streams of question marks, interspersed with occasional binary characters.
Cut-and-paste of the Opera 5.0 window produces just two characters: UT! That bug was fixed in Opera 6.0, but table cells are pasted on separate lines, making table reconstruction tedious.
The displayed table looks neat in a browser, because it has been carefully formatted with HTML aligned-table commands. However, JavaScript lacks control over number-to-string conversions, so it is not possible to get aligned decimal points without a lot of extra programming.
When you select the menu View -> Page Source in Netscape Navigator 4.xx, you see an HTML page in which the script has been replaced by its output, with the entire HTML table on one line. With version 6.xx, you see the JavaScript code. Attempts to insert line breaks with "\n" failed: they are discarded. Inserting an HTML line break, <BR>, after each row suppresses output entirely: there is no browser error window that I can find to see just why this happens. Version 6.xx has a menu path Tasks -> Tools -> JavaScript Console, but errors are not reported there either. Wrapping the table in a fixed-width font environment, <TT>...</TT>, fails to produce a table in that font.
With Microsoft Internet Explorer, the menu page View -> Source produces a popup window in which the JavaScript, rather than its output, is visible. An HTML line break kills the output, with no console window to report why. Using a fixed-width font has no effect, unless one changes the encoding via the menu path View -> Encoding -> User Defined.
With Opera, the View -> Source menu path produces a popup window in which the JavaScript, rather than its output, is visible. An HTML line break kills the output with no diagnostics, and fixed-width fonts are ignored.
The JavaScript programming environment appears distinctly unfriendly, since none of the tested browsers offers a way to get error reports. Writing and debugging of more complex code would then be nearly impossible. JavaScript is not supported by the other browsers (amaya, arena, chimera, elinks, grail, hotjava, links, lynx, netrik, w3m, and xmosaic) that I tried.
By default, mathematica uses a multiline representations of powers-of-ten, instead of the simple, and widespread, practice of E followed by an optionally-signed power-of-ten, but the latter can fortunately be obtained by wrapping the number with a CForm[] or FortranForm[] function call.
Although awk does not have an integer format item that supports comma group separators, a 10-line function adds that capability.
The perl example was produced using the awk-to-perl translator, a2p, rather than writing the code from scratch.
I would have included Modula-3 examples, but I have no installed compilers for that language. I made two installation attempts for the Compaq/DEC Modula-3 compiler, which was last updated on 17-Sep-1996; neither attempt produced a working compiler, so I abandoned the effort.
PHP (originally an acronym for Personal Home Pages) is a server-side scripting language used on millions of Web servers. Unlike Javascript, which is implemented on the client side, and thus is dependent on the environment at each end-user site, PHP processing is invisible to clients. Instead, a PHP-enabled Web server effectively replaces PHP program fragments by their run-time output in Web pages served to remote clients, although what really happens is that the PHP interpreter simply echoes material outside the program fragments, and the entire output is normally formatted as HTML.
The PHP program was straighforward, and closely resembles the C version, but without type declarations. For Web use, it is wrapped in an HTML file with file extension .php instead of .html, and the program is further wrapped in <?php ... ?>. Those delimiters are required even in standalone mode.
The availability of a C-like printf() function makes formatted output easy, but as in C and several other languages, precision is limited to the IEEE 754 64-bit format (C's double type).
Tests on our local Web server running PHP version 5 uncovered two problems:
With a bit more code, the second column in the output table could be made to line up properly, but I did not bother to do so for this simple program.
Scheme, as defined by the Revised(5) Report on Scheme and the IEEE Std 1178-1990 IEEE Standard for the Scheme Programming Language, is a minimalist language. It turns out that the GNU implementation, guile, has some extensions and deviations from the Standard: a (while ...) loop in addition to the Scheme (do ...) loop, no (printf ...) or (fprintf ...), and different arguments and behavior, and a much extended format item repertoire, for (format ...). PLT Scheme, with its mzscheme interpreter, seems closer to the Standard. The Bigloo Scheme implementation limits integers to 29 bits, and lacks (format ...) and (printf ...). Consequently, I have provided three incompatible implementations for these Scheme dialects.
During the development of these programs, I found that guile (versions 1.4, 1.5.6, and 1.7.2) could not correctly output subnormal floating-point numbers and large floating-point numbers, and had unusual nonstandard representations of Infinity and NaN. These have been reported to the guile developers, and fixes for all of them should be available in future releases. The Bigloo and PLT Scheme implementations appear to handle underflow, overflow, and subnormals correctly.
Regrettably, the two Common Lisp implementations (clisp version 2.34 and gcl versions 2.2.2 and 2.6.6) used for these tests below are even worse than the Scheme implementations: both trap overflow and zero divide, and clisp traps underflow to subnormals. These traps are irreparably embedded in the C-code interpreters for Common Lisp programs.
By default, clisp evaluates floating-point expressions in 32-bit (single-precision) arithmetic, whereas gcl uses 64-bit (double-precision) arithmetic. However, with suitable type coercions, and Fortran-style D-exponent letters on constants, clisp will use double-precision arithmetic.
A command shell is unexpected to be useful in a programming exercise such as this one, because shells traditionally have not supplied floating-point arithmetic. Among UNIX shells in the Bourne and C shell families (sh, ksh, bash, zsh, csh, tcsh, ...), only ksh93 (and later) and zsh provide such support. Unfortunately, ksh has proved difficult to build on many platforms, and none of the major UNIX vendors currently ships that version: most have ksh88 or earlier. Chapter 14 of our book Classic Shell Scripting gives instructions for building ksh93, and we have now been able to build and install it on most, but alas, not all, of our more than 20 flavors of Unix.
Floating-point arithmetic in both ksh93 and zsh uses the 64-bit IEEE 754 format, but alas, in ksh93, output truncates results to about nine digits, instead of the expected 16.
Surprisingly, GNU bash, which reimplements sh and ksh, does not provide floating-point arithmetic, even though ksh has had it for nearly a decade.
Compile and run these programs like this (on most UNIX systems):
Here is a comparison of the lines of code, excluding comments and blank lines, sorted by increasing line count.
| 12 | Fibonacci3.hoc |
| 13 | Fibonacci3.ml |
| 13 | Fibonacci3.math |
| 13 | Fibonacci3.py |
| 13 | Fibonacci3.red |
| 13 | Fibonacci3.tcl |
| 14 | Fibonacci3.maxima |
| 14 | Fibonacci3.sl |
| 15 | Fibonacci4.lsp |
| 15 | Fibonacci3.rb |
| 15 | Fibonacci3.rex |
| 15 | Fibonacci4.rexx |
| 16 | Fibonacci3.bc |
| 16 | Fibonacci3.lua |
| 16 | Fibonacci3.m |
| 17 | Fibonacci4.hoc |
| 17 | Fibonacci3.mf |
| 17 | Fibonacci3b.ml |
| 17 | Fibonacci3.mzs |
| 17 | Fibonacci3.splus |
| 18 | Fibonacci3.el |
| 18 | Fibonacci3.input |
| 18 | Fibonacci3.scm |
| 20 | Fibonacci3.f |
| 21 | Fibonacci3.lsp |
| 21 | Fibonacci3.maple |
| 21 | Fibonacci3p.p |
| 21 | fibona.pas |
| 21 | Fibonacci3.php |
| 22 | fibona.sai |
| 23 | Fibonacci3.c |
| 27 | Fibonacci3.awk |
| 27 | Fibonacci3.tex |
| 28 | Fibonacci3.gps |
| 28 | Fibonacci3.javascript |
| 28 | Fibonacci3.perl |
| 28 | fibona.alg |
| 29 | Fibonacci3.cs |
| 30 | Fibonacci2.java |
| 33 | Fibonacci3.adb |
| 33 | Fibonacci3.cc |
| 34 | Fibonacci4.java |
| 38 | Fibonacci4.c |
| 40 | Fibonacci3o.Mod |
| 41 | Fibonacci4.cs |
| 42 | Fibonacci3.javascript |
| 52 | Fibonacci4o.Mod |
| 60 | Fibonacci5.java |
| 65 | fibflt.b36 |
| 68 | fibdbl.b36 |
| 72 | Fibonacci3.java |