fortran

7 posts

RATFOR & FLECS - Emirp primes

For the final post in this series, let's write a real program in RATFOR and FLECS and see how they compare with the original FORTRAN. We'll be implementing the reverse-primes emirp program we did before.

FLECS version

C     FLECS PROGRAM TO DISPLAY EMIRPS  
C  
C     *** TEST IF A NUMBER IS PRIME ***  
      LOGICAL FUNCTION PRIME(N)
      INTEGER N
C     DEAL WITH NUMBERS <= 3  
      IF (N .LE. 1) GOTO 200
      IF (N .EQ. 2 .OR. N .EQ. 3) GOTO 100
C     CHECK IF DIVISIBLE BY 2 OR 3  
      IF (MOD(N,2) .EQ. 0) GOTO 200
      IF (MOD(N,3) .EQ. 0) GOTO 200
C     SEE IF DIVISIBLE BY 5, 7, ..., UP TO APPROX SQRT(N)  
      DO (I=5,999999,2)
      IF (I*I .GT. N) GOTO 100
      IF (MOD(N,I) .EQ. 0) GOTO 200
      FIN
 100  PRIME = .TRUE.
      RETURN
 200  PRIME = .FALSE.
      RETURN
      END
C  
C     *** REVERSE AN INTEGER'S DIGITS ***  
      INTEGER FUNCTION REVRSE(N)
      INTEGER N
      INTEGER M,R
C     M IS COPY OF N FROM WHICH WE TAKE DIGITS  
C     R IS REVERSED DIGITS  
      M = N
      R = 0
C     LOOP UNTIL NO MORE DIGITS  
      UNTIL (M .LT. 1)
C     TAKE LAST DIGIT FROM M AND APPEND TO R  
      R = R * 10
      R = R + MOD(M, 10)
      M = M / 10
      FIN
      REVRSE = R
      RETURN
      END
C  
C     *** TEST IF AN INTEGER IS AN EMIRP ***  
      LOGICAL FUNCTION EMIRP(N)
      INTEGER N
C     EXTERNAL FUNCTIONS  
      INTEGER REVRSE
      LOGICAL PRIME
C     R CONTAINS REVERSED DIGITS OF N  
      INTEGER R
      R = REVRSE(N)
C     N AND R MUST BOTH BE PRIME AND NOT THE SAME VALUE  
      IF (N .NE. R)
      IF (PRIME(N))
      IF (PRIME(R))
      EMIRP = .TRUE.
      RETURN
      FIN
      FIN
      FIN
      EMIRP = .FALSE.
      RETURN
      END
C  
C     *** DISPLAY AN INTEGER ***  
      SUBROUTINE SHOW(N)
      INTEGER N
      WRITE(6,50) N
 50   FORMAT(I10)
      RETURN
      END
C  
C  
C     *** MAIN ENTRY POINT ***  
C     I IS COUNT OF EMIRPS FOUND  
C     N IS NUMBER TO TEST  
C     EXTERNAL FUNCTION  
      LOGICAL EMIRP
      INTEGER I,N
      TEST-1
      TEST-2
      TEST-3
      STOP
C  
C     *** SHOW FIRST 20 EMIRPS ***  
      TO TEST-1
      N = 0
      I = 0
      WHILE (I .LT. 20)
      N = N + 1
      IF (EMIRP(N))
      CALL SHOW(N)
      I = I + 1
      FIN
      FIN
      FIN
C  
C     *** SHOW EMIRPS BETWEEN 7,700 AND 8,000 ***  
      TO TEST-2
      DO (N=7700,8000)
      IF (EMIRP(N)) CALL SHOW(N)
      FIN
      FIN
C  
C     *** SHOW 10,000TH EMIRP ***  
      TO TEST-3
      N = 0
      DO (I=1,10000)
      REPEAT UNTIL (EMIRP(N)) N = N + 1
      FIN
      CALL SHOW(N)
      FIN
C  
      END

Apart from the FORMAT specification and the PRIME function we've eliminated all line numbers. PRIME could be written without line numbers but with the multiple paths out of the function that would need their own RETURN I think it's better this way.

The internal procedures come in handy, eliminating the need for subroutines for TEST1-3, though this does make N and I global which makes me a little uneasy if this was a larger program.

We use the block structure often, with UNTIL, WHILE and REPEAT ... UNTIL; this simplifies code, though without indentation it's a little hard to follow; the output of the preprocesor is useful here to show what it thinks the indentation should be, for example:

  86           TO TEST-1
  87           .  N = 0
  88           .  I = 0
  89           .  WHILE (I .LT. 20)
  90           .  .  N = N + 1
  91           .  .  IF (EMIRP(N))
  92           .  .  .  CALL SHOW(N)
  93           .  .  .  I = I + 1
  94           .  .  ...FIN
  95           .  ...FIN
  96           ...FIN

The compiler diagnostics also helped a lot with catching errors with missing FINs.

RATFOR

Now let's try writing the RATFOR version.

######################################################################
# Ratfor program to display emirps
######################################################################

######### Test if a number is prime #########
logical function prime(n)  
    integer n  # Number to test

    # Deal with numbers <= 3
    if (n < 1) goto 200
    if (n == 2 | n == 3) goto 100

    # Check if divisible by 2 or 3
    if (mod(n,2) == 0) goto 200
    if (mod(n,3) == 0) goto 200

    # See if divisible by 5, 7, ..., up to approx sqrt(n)
    for (i = 5; i < 1000000; i = i + 2) {
        if (I*I > n) goto 100
        if (mod(n,i) == 0) goto 200
    }

 100  prime = .true.
      return
 200  prime = .false.
      return
end

######### Reverse an integer's digits #########
integer function revrse(n)  
    integer n  # Number to reverse
    integer m  # Copy of n from which we take digits
    integer r  # Reversed digits
    m = n
    r = 0
    while (m >= 1) {
        # Take last digit from m and append to r
        r = r * 10
        r = r + mod(m, 10)
        m = m / 10
    }
    revrse = r
    return
end

######### Test if an integer is an emirp #########
logical function emirp(n)  
    integer n       # Number to test
    integer revrse  # External function
    logical prime   # External function
    integer r       # Reversed digits of n
    r = revrse(n)
    emirp = .false.
    # n and r must both be prime and not the same value
    if (n .ne. r & prime(n) & prime(r)) {
        emirp = .true.
    }
    return
end

######### Display an integer #########
subroutine show(n)  
    integer n
    write(6,50) n
50  format(i10)  
    return
end

######### Show first 20 emirps #########
subroutine test1  
    logical emirp   # External function
    integer i       # Count of emirps found
    integer n       # Number to test
    n = 0
    for (i = 1; i <= 20; i = i + 1) {
        repeat {
            n = n + 1
        } until (emirp(n))
        call show(n)
    }
    return
end

######### Show emirps between 7,700 and 8,000 #########
subroutine test2  
    logical emirp   # External function
    integer n       # Number to test
    for (n = 7700; n <= 8000; n = n + 1) {
        if (emirp(n)) {
            call show(n)
        }
    }
    return
end

######### Show 10,000th emirp #########
subroutine test3  
    logical emirp   # External function
    integer i       # Count of emirps found
    integer n       # Number to test
    n = 0
    for (i = 1; i <= 10000; i = i + 1) {
        repeat {
            n = n + 1
        } until (emirp(n))
    }
    call show(n)
    return
end

######### Main entry point #########
call test1  
call test2  
call test3  
stop  
end  

I feel right at home with the braces and the C style for loops, though I miss the increment operator ++. prime would be much better if I could just return (.true.) but that does not work on the version of RATFOR on MTS so we keep the line numbers and gotos.

With the above, plus the free form input (which was supported on MTS FORTRAN anyway) and the operators like < it was easy to write. However, I got precisely zero diagnostics from the RATFOR preprocessor, with all my typos caught by the FORTRAN compiler, from which I'd have to find the problem in the original source. Easy enough in a small program but would be painful in larger ones.

Final thoughts

RATFOR and FLECS both make writing FORTRAN easier and more pleasant at the cost of an extra step in the development process, and I found both succeed at that. RATFOR is clearer and easier to get started with (especially coming from a C background today); the implementation is almost aggressively simple, as the authors admit in their paper, and I wonder how well it would scale for writing larger programs. FLECS has a more robust implementation but a more diffuse design, such as two versions of switch; features like printing a neatly indented output would certainly help on MTS or its contemporaries but the language lacks the cosmetic features that make RATFOR easier to read.

Neither are much used today; FORTRAN 77 and beyond took some of these ideas and built them into the core language. The idea of translating a richer language into a widely used but less expressive language is still alive though: think of Coffeescript or Typescript producing Javascript.

Further information

Full source code for these programs can be found on github.

RATFOR & FLECS - Language Features

FORTRAN meeting From the UM Computing Center Newsletter, Volume 5 Number 14, 24 September 1975, via Google Books. Proposal 2) seems to indicate a different preprocessor was being considered for UM as well as FLECS, I wonder if this was RATFOR or something else?

Hi and welcome back. Today let's continue our exploration of RATFOR and FLECS by comparing the language features they add to vanilla FORTRAN. The quotes below are from the RATFOR paper and FLECS manual, links to which are provided at the end of this post. Code samples for FORTRAN and FLECS are shown in upper case, RATFOR in lower case.

Design

RATFOR attempts to retain the merits of FORTRAN (universality, portability, efficiency) while hiding the worst FORTRAN inadequacies. The language is FORTRAN except for two aspects - [control flow and syntactic sugar] ... Throughout, the design principle which has determined what should be in RATFOR and what should not has been RATFOR doesn’t know any FORTRAN.

RATFOR focuses on control flow - if statements, blocks, looping - and cosmetics such as free form input, comments and other features that make FORTRAN more pleasant to write. By not knowing any FORTRAN, the design limits what features can be made available but also keeps it simple to implement and reduces the temptation to change FORTRAN into a different language altogether.

FLECS is a language extension of FORTRAN which has additional control mechanisms . These mechanisms make it easier to write FORTRAN by eliminating much of the clerical detail associated with constructing FORTRAN programs. FLECS is also easier to read and comprehend than FORTRAN.

FLECS also tries ti improve FORTRAN's control statements, taking ideas from several different languages including Pascal and Lisp. It has less cosmetic additions than RATFOR but adds the concept of internal procedures and includes features in the translator that help the programmer see the structure of their program.

Structure

RATFOR allows blocks of statements to be introduced within braces where FORTRAN would only allow a single statement. The fixed column format in classic FORTRAN is relaxed so any indentation is allowed. Multiple statements can appear on the same line if they are separated by semicolons.

if (x > 100) {  
   call error(x)
   err = 1; return
}

FLECS also has blocks which extend from the start of a control statement to the keyword FIN. It retains the fixed formatting of FORTRAN but prints a nicely indented view of the program when translating. So the example above would be entered as this in FLECS:

      IF (X .GT. 100)
      CALL ERROR(X)
      ERR = 1
      FIN

and the translator would print

IF (X .GT. 100)  
.  CALL ERROR(X)
.  ERR = 1
...FIN

This is useful when entering programs via cards where it is difficult to get indentation right.

It's possible to have a single statement after a control structure in which case the FIN is not needed:

IF (X .GT. 100) CALL ERROR(X)  

RATFOR comments are introduced with # and apply from that point to the end of the line, less restrictive than C in FORTRAN and FLECS which must be in the first column.

% will stop RATFOR processing the rest of the line, passing it through to FORTRAN directly. FLECS will look for a FLECS statement in column 7 and if found will translate the line; if not found it will pass through the whole line to FORTRAN.

Textual substitution

RATFOR allows constants to be set with define SYMBOL VALUE; any use of SYMBOL in the RATFOR program will be replaced with VALUE in the generated FORTRAN program.

include FILE will insert a copy of FILE at that point in the program, just like C's #include.

Operators

RATFOR allows the now-familiar symbols <, <=, !=, | etc to be used instead of .LT., .LE., .NE., .OR. etc. FLECS retains the FORTRAN operators.

Strings

Text in RATFOR programs in single or double quotes is converted to FORTRAN nH strings. Backslash escapes the next character. FLECS keeps FORTRAN strings.

Conditionals

FORTRAN has a simple iF statement where only one statement can be executed if the condition is true. RATFOR extends this by allowing else and nested ifs. An else clause is attached to the nearest if.

if (x > 0) {  
  if (x > 10)
    write(6, 1) x
  else
    write(6, 2) x
else  
  weite(6, 3)

FLECS has IF and for negative tests UNLESS. It also has WHEN ... ELSE for a single positive and negative test.

The switch statement added in RATFOR looks like C but does not have break; the switch is exited after each case or default is executed. FLECS's equivalent is SELECT, so comparing the two:

switch (x) {  
  case 1: y=3
  case 2, 3: y=5
  default y=0
}
      SELECT (X)
      (1) Y=3
      (2) Y=5
      (3) Y=5
      (OTHERWISE) Y=0
      FIN

FLECS has CONDITIONAL which looks a lot like LISP's cond:

      CONDITIONAL
      (X.LT.-5.0)  U = U+W
      (X.LE.1.0)   U = U+W+Z
      (X.LE.10.5)  U = U-Z
      (OTHERWISE)  U = 0
      FIN

Looping

The FORTRAN DO loop has to have a line number marking the point where the loop will restart:

      DO 10 i = 1, n
      x(i) = 0.0
      y(i) = 0.0
      z(i) = 0.0
 10   CONTINUE

RATFOR replaces this with a block:

do i = 1, n {  
  x(i) = 0.0
  y(i) = 0.0
  z(i) = 0.0
}

It also allows break to exit a loop early and next to restart the loop like C's continue. It can be followed by an integer to say how many levels to apply, so break 2 would move out of a two level do statement immediately.

RATFOR also adds a while and for statement that look like C's - these allow immediate exit from the statement if the condition is true on entry, unlike in FORTRAN DO where the statement is always executed at least once (in the IBM implementation at least) and the conditional is tested at the end of the statement. A version of C's do ... while is provided as repeat ... until.

The FLECS equivalent for the above do loop would be:

      DO (I = 1, N)
      X(I) = 0.0
      Y(I) = 0.0
      Z(I) = 0.0
      FIN

FLEC's WHILE construct is similar to RATFOR's, with the conditional tested before the loop starts. By using REPEAT WHILE the body of the loop is executed at least once and the test made at the end of the loop. UNTIL can be used instead of WHILE in both cases to indicate that the loop ends when the conditional becomes true

      X = 0
      UNTIL (X.EQ.5)
      X = X + 1
      FIN

Return

To return a value from a function in FORTRAN and FLECS you must assign a value to the name of the function:

INTEGER FUNCTION DECREMENT(I)  
INTEGER I  
DECREMENT = I - 1  
RETURN  
END  

In the RATFOR paper it sayd you can give return a value:

integer function decrement(i)  
integer i  
return (i-1)  
end  

However, note this is not supported in the version supplied with MTS - it will just pass through such a return statement causing an error from the FORTRAN compiler.

Internal procedures

FLECS allows a group of statements to be defined as a procedure with TO which can then be called by giving its name. No parameters are passed - it uses global variables to communicate. The below example will print 5.

      INTEGER X
      X = 1
      INCREMENT-IT
      DOUBLE-AND-INCREMENT
      WRITE(6,50) X
      STOP
 50   FORMAT(I10) 
      TO INCREMENT-IT X = X + 1
      TO DOUBLE-AND-INCREMENT
      X = X * 2
      INCREMENT
      FIN
      END

Procedure names must include at least one hyphen and recursion is not allowed.

Operation

RATFOR runs as a simple translator, taking a RATFOR input file and producing a FORTRAN output file that must then be fed to the FORTRAN compiler. FLECS, as modified at UM, will both translate and call the FORTRAN compiler, producing machine code output that can be run directly.

Error handling

RATFOR will catch some errors, such as missing closing braces, but will otherwise delegate problems with the program to the FORTRAN compiler to catch, as it does not understand FORTRAN syntax. This could be difficult to trace back to the source of the error as the FORTRAN compiler would show the error in the generated FORTRAN, not the RATFOR original.

FLECS will find syntax errors and remove them from the program, allowing translation to continue at the cost of possibly causing further errors; it will not move on to compilation in this case.

Implementation

Not surprisingly given its authors' roots, RATFOR was originally written in around 1000 lines of C using yacc. The authors say it took less than a week to implement. As C was not widely available in the mid 70's, a version of RATFOR in RATFOR was produced that would generate around 2500 lines of basic FORTRAN so it could be used anywhere.

The FLECS implementation comes in at around 2200 lines of FLECS and took around six months to develop according to comments in the source code.

Further information

See Kernighan's RATFOR paper or the FLECSUser's Manual (in component 673/22; I've uploaded a copy here) for more information on the languages.

RATFOR & FLECS - Introduction

Most programmers will agree that FORTRAN is an unpleasant language to program in, yet there are many occasions when they are forced to use it.

From the introduction to 'RATFOR — A Preprocessor for a Rational FORTRAN' by Brian W. Kernighan

FORTRAN was the lingua franca for mainframe programmers in the 1960s and 1970s, but as Kernighan states it's not always easy to program in - the main reasons are lack of good control structures and the fixed line format. As a result, a number of preprocessors were developed that translated enhanced code down to plain FORTRAN that could then be compiled anywhere a compiler was available.

In this series of posts, we'll look at two preprocessors available on MTS: RATFOR and FLECS. MTS also had OVERDRIVE, but this is not available on D6.0 due to copyright reasons.

Prerequisites

No special installation instructions to get these preprocessors running - just do the standard D6.0 setup as described in this guide and then sign on as a regular user such as ST01.

RATFOR

RATFOR was developed by Brian Kernighan at Bell Telephone Labs in 1974; its syntax was (not surprisingly) inspired by the C programming language, with keywords like for, while and until. It was used as the language for examples in Software Tools and became one of the most popular preprocessors in use. Versions are still available today that run on Unix systems.

Preprocessing using *RATFOR

The version on MTS is called *RATFOR and takes a RATFOR program as input on scards and writes FORTRAN source to spunch. The generated file can then be compiled with *FTN.

Hello world

Here's a terminal log of how to compile and run a simple hello world program in RATFOR. This assumes the source code is in file hello.r.

# $list hello.r

      1     # *** Simple hello world program ***
      2     #
      3     integer i
      4     for (i = 0; i < 5; i = i + 1)
      5     {
      6        write(6, 200)
      7     }
      8     stop
      9     200 format("Hello, world!")
     10     end

# $run *ratfor scards=hello.r spunch=-hello.f
Execution begins   21:56:08  
Execution terminated   21:56:08  T=0.004

# #list -hello.f

      1           INTEGERI
      2           CONTINUE
      3           I=0
      4     23000 IF(.NOT.(I.LT.5))GOTO 23002
      5           WRITE(6,200)
      6     23001 I=I+1
      7           GOTO 23000
      8     23002 CONTINUE
      9           STOP
     10     200   FORMAT(13HHello, world!)
     11           END

# $run *ftn scards=-hello.f spunch=-load
Execution begins   21:56:36  
 No errors in MAIN
Execution terminated   21:56:36  T=0.008

# #run -load
Execution begins   21:56:39  
 Hello, world!
 Hello, world!
 Hello, world!
 Hello, world!
 Hello, world!
Execution terminated   21:56:39  T=0.001  

FLECS

FLECS was written in the early 1970s by Terry Beyer at the University of Oregon. It provides a smaller set of control structures that RATFOR but the syntax is closer to FORTRAN. Keywords include IF...THEN...ELSE and CONDITIONAL and multi-line statements are supported. It does not appear to have been used much past the introduction of FORTRAN77, but a a version is still available today for HPUX.

Compiling using UNSP:FLX

At the time D6.0 was released, FLECS was unsupported at UM so is available as the file FLX in UNSP:. The preprocessor does not used scards and spunch; instead, all parameters need to be passed in to par. Unlike RATFOR, FLECS can call the FORTRAN compiler directly to generate object code. In the listing below, PAR=SOURCE=hello.fl,P=*SINK*,FTNSOURCE,LOAD=-load would read source from hello.fl, print diagnostics to *SINK including the FORTRAN source generated, and write compiled output to -load.

Hello world

Here's a terminal log of how to compile and run a simple hello world program in FLECS. This assumes the source code is in file hello.fl.

# $list hello.fl

      1     C *** SIMPLE HELLO WORLD PROGRAM ***
      2     C
      3           DO (I = 1,5)
      4           WRITE (6,20)
      5           FIN
      6           STOP
      7        20 FORMAT(13H HELLO, WORLD)
      8           END

# $run UNSP:FLX PAR=SOURCE=hello.fl,P=*SINK*,FTNSOURCE,LOAD=-load
Execution begins   21:47:21  
 FFI(CT206)


 (FLECS VERSION 22.38)  MTS Version CT155 21:47:21    JAN 21, 1916    Page   1

 MTS Line#        Indented Source Listing...

      1     C *** SIMPLE HELLO WORLD PROGRAM ***
      2     C
      3           DO (I = 1,5)
      4           .  WRITE (6,20)
      5           ...FIN
      6           STOP
      7        20 FORMAT(13H HELLO, WORLD)
      8           END

   0.001 seconds CPU time used.  Translation rate is 480000 lines per CPU minute.

 There were   NO MAJOR ERRORS and   NO MINOR ERRORS in the above module.
 No preprocessor errors in module  1.


  MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 01-21-16 21:47:21 PAGE P001

    0001              DO 99998 I = 1,5             3.000
    0002              WRITE (6,20)                 4.000
    0003        99998 CONTINUE                     5.000
    0004              STOP                         6.000
    0005           20 FORMAT(13H HELLO, WORLD)     7.000
    0006              END                          8.000
     *OPTIONS IN EFFECT*  ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP
     *OPTIONS IN EFFECT*  NAME = MAIN    , LINECNT =       57
     *STATISTICS*    SOURCE STATEMENTS =        6,PROGRAM SIZE =      344
     *STATISTICS*  NO DIAGNOSTICS GENERATED
 No errors in MAIN

 NO STATEMENTS FLAGGED IN THE ABOVE COMPILATIONS.
Execution terminated   21:47:21  T=0.018

# $run -load
Execution begins   21:47:31  
 HELLO, WORLD
 HELLO, WORLD
 HELLO, WORLD
 HELLO, WORLD
 HELLO, WORLD
Execution terminated   21:47:31  T=0.001  

Further information

The Wikipedia article on RATFOR has a basic introduction to language features and the history of its development. Kernighan's paper on RATFOR goes into more detail on the language.

Not much appears to exist on the Internet describing FLECS, but the D6.0 MTS tapes does include the complete User's Manual (in component 673/22) and the interface to MTS.

MTS Volume 6 describes the FORTRAN compilers on MTS, which are needed to compile the RATFOR preprocessor's output.

The great FORTRAN-LISP shootout

One of the reasons LISP did not become more popular for scientific computing was the perception that it was slower than FORTRAN. In this post I will compare the performance of the LISP and FORTRAN implementations on MTS using a simple problem.

The problem and some caveats

We'll use the finding emirp primes problem introduced in the FORTRAN posts on this blog as the subject. To make it simpler we will just look at the third test, finding the 10,000th emirp.

So some caveats before we start. Benchmarks are notoriously difficult to get right and this is only a single sample. I'm not an expert at optimising LISP or FORTRAN so there may well be better ways of implementing this. Finally, the UTILISP implementation is from the 1980s, MTS LISP is from the 1970s and the FORTRAN G and H compilers date from at least the 1960s so this is not a historically accurate recreation.

Having said that, it'll be interesting to see how they stack up, so let's start.

FORTRAN

We'll use the source code for the FORTRAN emirp implementation discussed before, taking out the first and second test.

We compile without optimisations first (ie using FORTRAN G) and then try again with FORTRAN H. The results look like this:

# run *ftn scards=emirp3.f
# Execution begins   11:39:08 
  No errors in PRIME 
  No errors in REVRSE 
  No errors in EMIRP 
  No errors in SHOW 
  No errors in TEST3 
  No errors in MAIN 
# Execution terminated   11:39:08  T=0.045 
# run -load
# Execution begins   11:39:10 
     948349
# Execution terminated   11:39:19  T=9.05 
#
# $run *ftn scards=emirp3.f par=opt=h
# Execution begins   11:39:49 

  **** No errors for  PRIME ****      SAT MAY 09/15 11:39:49 

  **** No errors for REVRSE ****      SAT MAY 09/15 11:39:49 

  **** No errors for  EMIRP ****      SAT MAY 09/15 11:39:49 

  **** No errors for   SHOW ****      SAT MAY 09/15 11:39:49 

  **** No errors for  TEST3 ****      SAT MAY 09/15 11:39:49 

  **** No errors for   MAIN ****      SAT MAY 09/15 11:39:49 
# Execution terminated   11:39:49  T=0.074 
# run -load
# Execution begins   11:39:54 
     948349
# Execution terminated   11:39:59  T=4.281 

So 9.05s for FORTRAN G and 4.28s for FORTRAN H.

MTS LISP

Let's try the original MTS LISP first. It's fairly easy to translate using the (GO) form to do looping.

;; Checks if n is prime using a simple division test
(DEFUN PRIME-P (N)
  (COND
   ;; Deal with numbers <= 3
   ((LESS N 2) NIL)
   ((LESS N 4) T)
   ;; Check if divisible by 2 or 3
   ((EQUAL (REMAIN N 2) 0) NIL)
   ((EQUAL (REMAIN N 3) 0) NIL)
   ;; See if divisible by 5, 7, ..., up to approx sqrt(n)
   (T (PROG (I RESULT)
            (SETQ I 5)
A           (SETQ RESULT (COND ((GREATER (TIMES I I) N) 'PRIME)  
                               ((EQUAL (REMAIN N I) 0) 'NOT-PRIME)
                               (T (SETQ I (ADD I 2)) NIL)))
            (COND ((NOT RESULT) (GO A)))
            (EQ RESULT 'PRIME)))))

;; Return a number that has the reversed digits in n
(DEFUN REVERSE-DIGITS (N)
  (PROG (REVERSED)
        (SETQ REVERSED 0)
        ;; Take last digit from n and append to reversed
        (PROG ()
B             (SETQ REVERSED (TIMES REVERSED 10))  
              (SETQ REVERSED (ADD REVERSED (REMAIN N 10)))
              (SETQ N (IDIVIDE N 10))
              (COND ((GREATER N 0) (GO B))))
        ;; CAR needed to evaluate REVERSED otherwise PROG will return
        ;; the symbol REVERSED
        (CAR REVERSED)))

;; Check if a number is an emirp, ie both it and its reversed digits
;; are prime.
(DEFUN EMIRP-P (N)
  (PROG (REVERSED)
        (SETQ REVERSED (REVERSE-DIGITS N))
        (AND (NOT (EQUAL N REVERSED)) (PRIME-P N) (PRIME-P REVERSED))))

(DEFUN EMIRP-10K ()
  (PROG (N COUNT LAST)
        (SETQ N 0)
        (SETQ COUNT 0)
        (SETQ LAST 0)
C       (COND ((EMIRP-P N)  
               (SETQ COUNT (ADD COUNT 1)) (SETQ LAST N)))
        (SETQ N (ADD N 1))
        (COND ((LESS COUNT 10000) (GO C)))
        (CAR LAST)))

(EMIRP-10K)
(STOP)

We'll first execute it through the interpreter:

# $run *lisp scards=emirp.l
# Execution begins   12:41:27
      WELCOME TO LISP/MTS
 >    PRIME-P
 >    REVERSE-DIGITS
 >    EMIRP-P
 >    EMIRP-10K
* GC: COLLECTED 23370 CELLS
...
* GC: COLLECTED 23370 CELLS
>    948349
# Execution terminated   12:47:15  T=347.517 

Next we'll compile the functions. Add (COMPILE PRIME-P) and so on for all the functions and rerun

# $run *lisp scards=emirp.l
     WELCOME TO LISP/MTS
>    PRIME-P
>    REVERSE-DIGITS
>    EMIRP-P
>    EMIRP-10K
* COMPILE        04-17-79 RESTORED
>    (PRIME-P)
>    (REVERSE-DIGITS)
>    (EMIRP-P)
>    (EMIRP-10K)
* GC: COLLECTED 23370 CELLS
...
* GC: COLLECTED 23370 CELLS
>    948349
# Execution terminated   12:56:39  T=72.347 

OK, so we are down from 348s to 72s (the compilation itself takes less than 0.5s). Not bad, but still 17 times slower than FORTRAN H.

UTILISP

Let's see how UTILISP compares. We do need to change the program in quite a few places to match the different built in functions (eg + instead of ADD) but at least we can now use the do macro to create more expressive loops.

;; Checks if n is prime using a simple division test
(defun prime-p (n)
  (cond
   ;; Deal with numbers <= 3
   ((<= n 1) nil)
   ((<= n 3) t)
   ;; Check if divisible by 2 or 3
   ((= (remainder n 2) 0) nil)
   ((= (remainder n 3) 0) nil)
   ;; See if divisible by 5, 7, ..., up to approx sqrt(n)
   (t (do ((i 5 (+ i 2)) (result nil))
          (result (eq result 'prime))
          (setq result (cond ((> (* i i) n) 'prime)
                             ((= (remainder n i) 0) 'not-prime)
                             (t nil)))))))


;; Return a number that has the reversed digits in n
(defun reverse-digits (n)
  (do ((reversed 0) (m n (// m 10)))
      ((< m 1) reversed)
      ;; Take last digit from n and append to reversed
      (setq reversed (* reversed 10))
      (setq reversed (+ reversed (remainder m 10)))))


;; Check if a number is an emirp, ie both it and its reversed digits
;; are prime.
(defun emirp-p (n)
  (let ((reversed (reverse-digits n)))
    (and (neq n reversed) (prime-p n) (prime-p reversed))))


;; Return the 10,000th emirp
(defun emirp-10k ()
  (do ((n 0 (+ n 1)) (count 0) (last 0))
      ((eq count 10000) last)
      (cond ((emirp-p n) (setq count (+ count 1)) (setq last n)))))

;; Run the test
(print (emirp-10k))
(quit)

And let's run it through the interpreter.

# $run *utilisp 0=emirp.ul
# Execution begins   18:18:02 
  948349
# Execution terminated   18:21:30  T=208.35 

So that's 40% quicker than the MTS LISP interpreted version. Let's try compiling the forms, which we can just do by putting

(compile prime-p reverse-digits emirp-p emirp-10k)

after the definitions and then adding the fix parameter to the command line to reserve memory for the compiler (again compilation takes less than 0.5s of the total):

$ run *utilisp 0=emirp.ul par=fix=100
# Execution begins   18:32:27 
  ; Loading Utilisp Compiler System Version(82-11-09)
  948349
# Execution terminated   18:32:54  T=27.788 

Much better - 3 times faster than the MTS LISP compiled version but 7 times slower still than FORTRAN H.

Summary of results and final thoughts

Program Time (s)
FORTRAN G 9.05
FORTRAN H 4.28
MTS LISP interpreted 347.52
MTS LISP compiled 72.35
UTILISP interpreted 208.35
UTILISP compiled 27.79

Compiled UTILISP puts up a good showing and it's definitely usable for light numerical tasks but in the end FORTRAN wins on pure performance. The UTILISP implementation is much faster than MTS LISP for both interpreted and compiled code, which you'd expect given the amount of knowledge creating fast LISP implementations that was built up between the release of the two programs.

Another factor evident from this test was the lack of LISP standardisation (pre-Common Lisp) makes porting programs between implementations a chore, whereas FORTRAN was relatively standardised.

Any suggestions for optimising the code? Please let me know in the comments.

In the next post, in a couple of weeks, we'll look at PL/1.

Further information

Full source code for the programs shown here is on github.

FORTRAN - Calculating emirps

Let's look at an example of a complete program in FORTRAN and compare running it using the different compilers available on MTS.

The problem

Emirps are prime numbers that when their digits are reversed they make a different prime. An example is 13, which when reversed is the prime 31. The task on Rosetta Code is to:

  • show the first twenty emirps.
  • show all emirps between 7,700 and 8,000
  • show the 10,000th emirp

Writing the program

I think it's interesting to describe how I wrote the program and got it to work. On modern systems I would tend to dive in and start writing small pieces of code and associated tests, TDD style. On MTS I favour writing the program offline first, import it into MTS and then get it working - this is partially due to being less comfortable with the MTS editor but also to match the usage pattern for MTS when it was in use; then terminal time was expensive so you would try to get more done before reaching for the computer.

So I started by sketching out what routines I'd need
* Three subroutines TEST1-3to test each of the tasks given in the problem * A function EMIRP to test if a number was an emirp. This would reverse the digits and if different from the original, test if both were prime. * A function PRIME to test if a number was a prime * A function REVRSE to reverse a number

I used Emacs on my host OS to write the routines - was pleasantly surprised to find it still supports FORTRAN 66 syntax. This took around 45 minutes.

The prime number test I used was a brute force approach so not particularly efficient:

C     *** TEST IF A NUMBER IS PRIME ***  
      LOGICAL FUNCTION PRIME(N)
      INTEGER N
C     DEAL WITH NUMBERS <= 3  
      IF (N .LE. 1) GOTO 200
      IF (N .EQ. 2 .OR. N .EQ. 3) GOTO 100
C     CHECK IF DIVISIBLE BY 2 OR 3  
      IF (MOD(N,2) .EQ. 0) GOTO 200
      IF (MOD(N,3) .EQ. 0) GOTO 200
C     SEE IF DIVISIBLE BY 5, 7, ..., UP TO APPROX SQRT(N)  
 10   DO 20 I=5,999999,2
         IF (I*I .GT. N) GOTO 100
         IF (MOD(N,I) .EQ. 0) GOTO 200
 20   CONTINUE
 100  PRIME = .TRUE.
      RETURN
 200  PRIME = .FALSE.
      RETURN
      END

Reversing the digits was done by taking the modulus of the number with 10 and then appending each digit to the end:

C     *** REVERSE AN INTEGER'S DIGITS ***  
      INTEGER FUNCTION REVRSE(N)
      INTEGER N
      INTEGER M,R
C     M IS COPY OF N FROM WHICH WE TAKE DIGITS  
C     R IS REVERSED DIGITS  
      M = N
      R = 0
C     LOOP UNTIL NO MORE DIGITS  
 10   IF (M .LT. 1) GOTO 100
C     TAKE LAST DIGIT FROM M AND APPEND TO R  
      R = R * 10
      R = R + MOD(M, 10)
      M = M / 10
      GOTO 10
 100  REVRSE = R
      RETURN
      END

The emirp function was then a simple case of calling the two functions:

C     *** TEST IF AN INTEGER IS AN EMIRP ***  
      LOGICAL FUNCTION EMIRP(N)
      INTEGER N
C     EXTERNAL FUNCTIONS  
      INTEGER REVRSE
      LOGICAL PRIME
C     R CONTAINS REVERSED DIGITS OF N  
      INTEGER R
      R = REVRSE(N)
C     N AND R MUST BOTH BE PRIME AND NOT THE SAME VALUE  
      IF (N .EQ. R) GOTO 200
      IF (.NOT. PRIME(N)) GOTO 200
      IF (.NOT. PRIME(R)) GOTO 200
 100  EMIRP = .TRUE.
      RETURN
 200  EMIRP = .FALSE.
      RETURN
      END

The test functions are simple loops: here I show just the second test as an example.

C     *** SHOW EMIRPS BETWEEN 7,700 AND 8,000 ***  
      SUBROUTINE TEST2
C     EXTERNAL FUNCTION  
      LOGICAL EMIRP
C     N IS NUMBER TO TEST  
      INTEGER N
 10   DO 100 N=7700,8000
         IF (EMIRP(N)) CALL SHOW(N)
 100  CONTINUE
      RETURN
      END

I factored out the printing to another routine as each test needs to do this:

C     *** DISPLAY AN INTEGER ***  
      SUBROUTINE SHOW(N)
      INTEGER N
      WRITE(6,50) N
 50   FORMAT(I10)
      RETURN
      END

Getting it running with *IF66

Next I imported the program into MTS (using the emulated card reader).

I then ran *IF66 and did /compile emirp.f. Lots of syntax errors - I forgot the maximum length for identifiers was 7 characters and that I needed to declare functions used externally.

Next I tried to test the individual routines - here *IF66 was very useful, as I could call eg the reverse routine with some test input and see its output immediately. Some more problems to fix here, eg function parameters seem to be pass by reference in FORTRAN. An error in the DO-loop in the prime testing function lead to an infinite loop that I could not get out of with PA2, so a trip to the operator's console to STOP the MTS session was needed.

Finally I got it running after about another hour. The first two tasks were quick to execute but the third, to find the 10,000th emirp, ran for over 10 minutes without returning so I killed it. Time to run this using the compilers and see how they perform.

Using the compiler

First I tried it using FORTRANG:

# $run *ftn scards=emirp.f
# Execution begins   20:05:30 
  No errors in PRIME 
  No errors in REVRSE 
  No errors in EMIRP 
  No errors in SHOW 
  No errors in TEST1 
  No errors in TEST2 
  No errors in TEST3 
  No errors in MAIN 
# Execution terminated   20:05:30  T=0.034
# $run -load
# Execution begins   20:05:34 
         13
...
     948349
# Execution terminated   20:05:43  T=8.82 

So just over 0.03 seconds of CPU time to compile and nearly 9s to run all tests.

Let's try it with FORTRANH by using the opt parameter to *ftn. Note that the output of the compiler is quite different as FORTRANH is a completely different product:

# $run *ftn scards=emirp.f par=opt=h
# Execution begins   20:06:12 
  **** No errors for  PRIME ****      WED DEC 17/14 20:06:13
  **** No errors for REVRSE ****      WED DEC 17/14 20:06:13
  **** No errors for  EMIRP ****      WED DEC 17/14 20:06:13
  **** No errors for   SHOW ****      WED DEC 17/14 20:06:13
  **** No errors for  TEST1 ****      WED DEC 17/14 20:06:13
  **** No errors for  TEST2 ****      WED DEC 17/14 20:06:13
  **** No errors for  TEST3 ****      WED DEC 17/14 20:06:13
  **** No errors for   MAIN ****      WED DEC 17/14 20:06:13
# Execution terminated   20:06:12  T=0.092 
# $run -load
# Execution begins   20:06:30 
         13
...
     948349
# Execution terminated   20:06:35  T=4.607 

Compilation took almost three times as long, but the program ran 48% faster so there's definitely some good optimisation taking place. The compilation time difference may seem insignificant, but in an environment where you are paying for CPU time and working on a large program that requires many recompiles during development the costs would mount up.

Next step would probably be to improve the algorithm for testing for primeness, but I'll leave that for now.

Final thoughts on FORTRAN

FORTRAN 66 has not aged well - its fixed format lines, lack of character/strings, unstructured design and implicit variable types make it hard to use for large scale programming. But it's a simple language to pick up and use, and you can see that if you want to solve a numerical problem it can be very productive. Coupled with the optimising compiler and support tools available you can understand why it was the lingua franca of its time, rather like C today.

Further information

Full source code for this program is on github.

FORTRAN - Language features

In this post I'll look at the FORTRAN 66 language as available on MTS in more detail.

Source code format

FORTRAN originally had a fixed format for source code lines and used all upper case letters; as *FTN does not have this restriction we'll use lower case and free format in the examples below to make them easier to read, but note that this means these examples may not work on other, non-MTS compilers.

Update 4-Dec-2015: As noted by Mike Stramba, *FTN does require source code to be in upper case, although *IF66 is happy with upper or lower case.

One unique feature is that whitespace is not important, so the variable BXX could also be referred to as B XX. Keywords are context sensitive, so you can write IF IF.GT.3 J=2 where the second IF is actually a variable.

Types and identifiers

FORTRAN supports INTEGER, REAL. LOGICAL and COMPLEX number types. Character types are not support directly, but can be accessed through FORMAT statements (see below).

Variable names can have up to six characters and must start with a letter. Variables can be explicitly declared with their type; if not declared, variables starting with letters I to N are treated as integers, otherwise they are treated as reals.

Variables can be initialised with the DATA statement.

integer count, count1  
real temp  
complex point  
data i,count,temp,point/0, 10, 32.0, cmplx(1,2)/  

Arrays and matrices can be created with DIMENSION and initialised with DATA then individual values can be accessed using (), first index 1, in column major order, so in the below X(2, 1) will be intialised with the value 2.0.

dimension x(2,2)  
data x/1.0, 2.0, 3.0, 4.0/  

Operators and expressions

The usual arithmetic operators appear, including ** for exponentiation. Logical operators look like .LT. for less than or .AND. for boolean-and.

x = 2.0 ** 4.0  
i = 5 + 3  
j = x / i  
y = x / i  

In the result above, J will be set to 2 and y to 2.0.

Statements

Branching can be done with GOTO; there are three forms:

goto 100  
goto l  
goto (10,20,30) ll  
  • The first will jump unconditionally to the line labelled with 100.
  • The second will jump to the label based on the value of L
  • The third will jump to line 10 if ll = 1, 20 if ll = 2 etc

The DO statement allows looping

j = 0  
do 100 i=10,20,5  
j = j + i  
100 continue  
k = j * 2  

This will initialise I with 10 and then run the contents of the loop for values 10, 15, 20. After that it will jump to the line after the label 100. So for this example, K will end up with the value 90.

It's possible to nest loops and jump out of loops, but not jump into them.

The IF statement has two forms, logical and arithmetic. The example of the logical version below will set k to 3 if i < j.

if (i.lt.j) k=3  

The arithmetic version will branch to one of three labels based on whether the expression being tested is negative, zero or positive. For the example below it will branch to label 30 as the expression is positive.

J = -3  
if (j + 4) 10,20,30  

Input/output and formats

Inputs (READ) and outputs (WRITE) are associated with a FORMAT statement, a logical I/O unit and a list of variables. For example:

read(5,100) i1,x1,y1  
100   format(i2,f5.2,f2.0)  

with the input

1234.5678.90  

will read from I/O unit 5 (the keyboard) i1=12, x1=34.56, y1=78.0.

Character constants, with a leading length identifier, can be used in READ or WRITE statements as we saw in the hello world example:

write (6,200)  
200 format(13h HELLO, WORLD)  

prints "HELLO, WORLD" on unit 6, the console.

Functions and subroutines

FORTRAN comes with a large number of built in functions for mathematical processing eg MOD, ABS, SIN. On MTS there is also a library of operating system specific functions that can be called but I will not go into these here.

You can define your own functions to pass parameters by reference and optionally return a value: for example to calculate the average of two values:

real function avg2(x,y)  
real x, y  
avg2 = (x + y) / 2.0  
return  
end  

This could be called with

z = avg2(3.0, x)  

There is also subroutines which run by using call. For example, the below would set x to 2.0:

x = 1.0  
call add1(x)

subroutine add1(x)  
real x  
x = x + 1.0  
return  
end  

Further information

FORTRAN IV reference page is a great summary of FORTRAN 66 syntax.

A contemporary IBM references for FORTRAN can be found at Bitsavers.

MTS Volume 6 describes *FTN and *IF66 in full detail, along with the other FORTRAN tools available.