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:

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

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.

Comments

comments powered by Disqus