Category Archives: Fortran

Code I’ve written in Fortran 95.

Calculating Logarithms of an Arbitrary Base

Posted on May 12, 2012 by Jonathan Landrum
Reply

Fortran includes a couple of intrinsic functions for calculating logarithms: LOG(), which calculates the loge — the natural log or ln — of the number, and LOG10(), which calculates the log10 of the number. Why the former isn’t LN() and the latter LOG() is beyond me, but I digress.

These two logarithm functions are normally all one would need for their calculations, but say you need the logarithm of a fractional number to a fractional base? You can do that, since the definition of a logarithm has been extended to the set of all positive numbers, and can be calculated with some simpler math:

Equation for solving logarithms of an arbitrary base

Since calculating logarithms to a base of 10 is trivial, we can simply use this formula within a function called LogB() to return the logarithm of any positive number to any positive base:

49
50
51
52
53
54
55
56
57
58
59
60
61
62

! ----------------------------------------------------------
! FUNCTION LogB:
! Calculates the logarithm of a number to an arbitrary base.
! ----------------------------------------------------------
REAL FUNCTION LogB(n, b)
	IMPLICIT NONE
 
	REAL :: n, b
 
	LogB = LOG10(n) / LOG10(b)
 
	RETURN
 
END FUNCTION LogB

Using this function within a program is as simple as calling it: LogB(number, base). So simple, in fact, that I’m not even going to reproduce the main program here. This function is all that matters in the scope of this article, anyway. But why do the numbers in the actual parameter list have to be positive? Because negative logarithms are undefined:

Graph of the natural logarithm of x

Graph of the natural logarithm of x

As the number to be calculated approaches 0, the value decreases without bound. The y-axis is actually a vertical asymptote for it. In the function I wrote above, I didn’t put bounds checking in place. But Fortran is smart enough to deal with poor input, yielding NaN for negative input, and -Infinity for 0. And people wonder why I love this language.

~Jonathan

Translating Programs Between Languages

Posted on April 29, 2012 by Jonathan Landrum
Reply

I’ve written on this site so far in a few of the languages I know, and today I want to demonstrate that any language you like is usually a perfectly valid language to use. Each has their own niche that sets them apart, but generally they all do the same thing.

I was recently given an assignment to write a program that prints an hourglass shape using asterisks, and the user input should specify how many rows of asterisks there are on each end of the figure. For instance, a user input of 3 should output the following figure:

*****
 ***
  *
 ***
*****

Three rows on top, and three on bottom (counting the center asterisk as one each time.) I wrote the assignment in C++, since that’s what the assignment specifically called for. But I have also translated the program into Java and Fortran to demonstrate my point.

Java is only a slight change from C++, so it’s obvious. But Fortran is absolutely foreign in comparison; can an algorithm from C++ really be translated to Fortran? If you keep the basic premise of your algorithm’s logic intact, then yes, it can easily be done. You just have to be familiar with the differences in each language’s logical structures. Let’s start with the original C++ code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

#include <iostream>
using namespace std;
 
int main () {
    int input = 0;
 
    cout << "Enter the size of your hourglass: ";
    cin  >> input;
 
    for (int c = input - 1; c > 0; --c) {
        for (int i = 0; i < (input - c) - 1; ++i)
            cout << " ";
        for (int i = (2 * c + 1); i > 0; --i)
            cout << "*";
        cout << endl;
    } // End top for loop
    for (int c = 0; c < input; ++c) {
        for (int i = 0; i < (input - c) - 1; ++i)
            cout << " ";
        for (int i = (2 * c + 1); i > 0; --i)
            cout << "*";
        cout << endl;
    } // End bottom for loop
} // End main

Two for loops with diligent attention given to the counter variables and the conditions that make the loops exit, and you’re done. Simple, right? Converting this to Java is just as simple:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

import java.util.Scanner;
 
public class hourglass {
	public static void main (String[] args) {
		int input = 0;
		Scanner scan = new Scanner(System.in);
 
		System.out.print ("Enter the size of your hourglass: ");
		input = scan.nextInt();
 
		for (int c = input - 1; c > 0; --c) {
			for (int i = 0; i < (input - c) - 1; ++i)
				System.out.print (" ");
			for (int i = (2 * c + 1); i > 0; --i)
				System.out.print ("*");
			System.out.println ();
		} // End top for loop
		for (int c = 0; c < input; ++c) {
			for (int i = 0; i < (input - c) - 1; ++i)
				System.out.print (" ");
			for (int i = (2 * c + 1); i > 0; --i)
				System.out.print ("*");
			System.out.println ();
		} // End bottom for loop
	}
}

So we’ve got our program working in two different languages now (which wasn’t all that difficult, really.) Now let’s give Fortran a try. This one is going to be a bit tricky; we can’t just copy-and-paste the code and change the differences. Fortran is altogether different. But the principles of making this function work are identical. Watch:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

PROGRAM hourglass
	IMPLICIT NONE
 
	INTEGER :: input = 0
	INTEGER :: c
	INTEGER :: i
 
	WRITE (*,'(a)',advance='no') 'Enter the size of your hourglass: '
	READ  (*,*) input
 
	c = input - 1
	DO WHILE (c > 0)
		i = 0
		DO WHILE (i < (input - c) - 1)
			WRITE (*,'(a)',advance='no') ' '
			i = i + 1
		END DO
		i = (2 * c + 1)
		DO WHILE (i > 0)
			WRITE (*,'(a)',advance='no') '*'
			i = i - 1
		END DO
		c = c - 1
		WRITE (*,*)
	END DO
	c = 0
	DO WHILE (c < input)
		i = 0
		DO WHILE (i < (input - c) - 1)
			WRITE (*,'(a)',advance='no') ' '
			i = i + 1
		END DO
		i = (2 * c + 1)
		DO WHILE (i > 0)
			WRITE (*,'(a)',advance='no') '*'
			i = i - 1
		END DO
		c = c + 1
		WRITE (*,*)
	END DO
END PROGRAM hourglass

It’s longer than the other two; nearly twice as long. And there’s a lot of line return control going on. But it can be done, and without any hacks or Fortran-fu. It’s a simple translation of the exact for loop structure from before into a Fortran DO WHILE structure. And that, my friend, is how you translate code from one language to another. Know both languages first, know how to work logical constructs in both, and you’ll do fine. Happy coding.

~Jonathan

Creating a Decimal-to-Binary Converter in Fortran

Posted on March 7, 2012 by Jonathan Landrum
Reply

Screenshot of the program converting 65,535 to binary

I’ve been re-reading my Fortran book from last semester and have been finding some jewels to code. One such assignment was to make a decimal-to-binary converter that accepts numbers up to 1,023. That wasn’t good enough for me, so when I wrote my version of it, I made it go all the way to 65,535 (2 bytes).

The algorithm works the same way no matter how large you go; the only limitation is how large an integer your Fortran compiler will allow. I am really proud of this one; I’ve never devised an algorithm for converting a base-10 number to base-2, so I had to think about it for a bit. I can even imagine turning it into an external method that has an array of arbitrary size, and uses a for loop to assign and write the binary digits. You would just need to pass the array size as an actual parameter. It could easily be part of a larger program then.

~Jonathan

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

PROGRAM decimalToBinary
! ------------------------------------------------------------------------------
! Programmer:     Jonathan Landrum
!
! Purpose:        This program takes Integer input from the user of a decimal
!                 number, and returns the corresponding binary version.
!
! Assumptions:    1.) Input is positive and less than 65,536. Numbers equal to
!                     or greater than this will cause a rollover, producing
!                     incorrect results.
!                 2.) Input is in Integer form. Real input will be rounded to
!                     the nearest Integer.
! 
! Revisions:      Date		Programmer		Description of Change
!                 ===========	=====================	========================
!                 07 Mar 2012	Jonathan Landrum	Original code.
! ------------------------------------------------------------------------------
 
	IMPLICIT NONE
 
	! DATA DICTIONARY: Declared variables
	REAL                     :: input	! Input received from user
	INTEGER                  :: number	! Input converted to integer
	INTEGER, DIMENSION(0:15) :: array	! Array to store the answer
 
	! Introduce the program
	WRITE (*,*)
	WRITE (*,*) '* * * * * * * * * * * * * * * * * * * * * * * * * * * * *'
	WRITE (*,*) '*                                                       *'
	WRITE (*,*) '*          Fortran Decimal to Binary Converter          *'
	WRITE (*,*) '*                                                       *'
	WRITE (*,*) '* * * * * * * * * * * * * * * * * * * * * * * * * * * * *'
	WRITE (*,*)
	WRITE (*,*) 'This program takes an integer in base-10 and returns the'
	WRITE (*,*) 'corresponding number in base-2.'
	WRITE (*,*)
	WRITE (*,*) 'Output is arranged into nybbles and bytes like so:'
	WRITE (*,*) '1001 0100  1101 1110'
	WRITE (*,*)
	WRITE (*,*) 'Assumptions:  1.) Input is positive and less than 65,536.'
	WRITE (*,*) '                  Numbers equal to or greater than this'
	WRITE (*,*) '                  will cause a rollover, producing'
	WRITE (*,*) '                  incorrect results.'
	WRITE (*,*) '              2.) Input is in Integer form. Real input'
	WRITE (*,*) '                  will be rounded to the nearest Integer.'
	WRITE (*,*)
	WRITE (*,*) '---------------------------------------------------------'
	WRITE (*,*)
 
	! Get data values from the user
	WRITE (*,*) 'What number would you like to convert to binary?'
	READ  (*,*) input
	WRITE (*,*)
 
	! Round to an integer if our input has a fractional component
	number = NINT(input)
 
	! Perform the calculations
	array(0)  = MOD(number, 2)	! This is actually MOD(number / 1, 2),
	array(1)  = MOD(number / 2, 2)	! but that's silly and wasteful.
	array(2)  = MOD(number / 4, 2)
	array(3)  = MOD(number / 8, 2)
	array(4)  = MOD(number / 16, 2)
	array(5)  = MOD(number / 32, 2)
	array(6)  = MOD(number / 64, 2)
	array(7)  = MOD(number / 128, 2)
	array(8)  = MOD(number / 256, 2)
	array(9)  = MOD(number / 512, 2)
	array(10) = MOD(number / 1024, 2)
	array(11) = MOD(number / 2048, 2)
	array(12) = MOD(number / 4096, 2)
	array(13) = MOD(number / 8192, 2)
	array(14) = MOD(number / 16384, 2)
	array(15) = MOD(number / 32768, 2)
 
	! Formats for output
1	FORMAT (I1)	! Format for the digits of the binary number
2	FORMAT (1X)	! Format for the space between nybbles and bytes
 
	! Return the result
	WRITE (*,1,advance='no') array(15)
	WRITE (*,1,advance='no') array(14)
	WRITE (*,1,advance='no') array(13)
	WRITE (*,1,advance='no') array(12)
	WRITE (*,2,advance='no')
	WRITE (*,1,advance='no') array(11)
	WRITE (*,1,advance='no') array(10)
	WRITE (*,1,advance='no') array(9)
	WRITE (*,1,advance='no') array(8)
	WRITE (*,2,advance='no')
	WRITE (*,2,advance='no')
	WRITE (*,1,advance='no') array(7)
	WRITE (*,1,advance='no') array(6)
	WRITE (*,1,advance='no') array(5)
	WRITE (*,1,advance='no') array(4)
	WRITE (*,2,advance='no')
	WRITE (*,1,advance='no') array(3)
	WRITE (*,1,advance='no') array(2)
	WRITE (*,1,advance='no') array(1)
	WRITE (*,1,advance='no') array(0)
 
	! Exit the program
	WRITE (*,*)
	WRITE (*,*)
	WRITE (*,*) '\\//_ Live long and prosper.'
	WRITE (*,*)
END PROGRAM decimalToBinary

Formatting Output with Fortran

Posted on March 5, 2012 by Jonathan Landrum
Reply

An interesting thing about Fortran is its inherent lack of ability to handle strings. Not that the language cannot be made to work with strings; it’s just not designed to do so. Fortran really shines when you are dealing with numbers on the command line. Aside from that, it’s not the language you should use.

However, even in those cases, you still need text to make prompts, or you need your numbers to be formatted in a certain way. With Fortran, you can declare a FORMAT statement and assign a certain format to it. A label is applied to the left of the FORMAT, and it looks kind of like the line numbers in BASIC, but it’s not. It’s a label, and you reference that label within a READ or WRITE statement.

The premise of this program is there is a bracketed income tax system in Australia, along with a 1.5% Medicare tax on all income. Your job is to write a program that calculates how much tax a person will owe based on their income. The brackets break down like so:

Bracketed Income Tax System

I’ve used a FORMAT statement on line 35 to show how they’re used (you can read more about them at Dr. Ching-Kuang Shene’s site). That line essentially says “Expect two spaces of text, and ten spaces of a floating point number truncated to two decimal places, and put one space before and between them.” (Those two text spaces are used for the dollar sign.) What you get is this:

Screenshot of the finished product

The program itself is actually quite easy to write. It’s the FORMAT’s that will get you. But once you’ve gotten comfortable with them, they will be second nature, as well.

~Jonathan

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

PROGRAM tax
! ------------------------------------------------------------------------------
! Programmer:     Jonathan Landrum
!
! Purpose:        This program calculates the total taxes owed by an Australian
!                 citizen, based on their taxable income.
!
! Assumptions:    1.) Input is in Real form, truncated to 2 decimal places.
!                     Integer input will be assumed to be even dollars.
! 
! Revisions:      Date		Programmer		Description of Change
!                 ===========	=====================	========================
!                 05 Mar 2012	Jonathan Landrum	Original code.
! ------------------------------------------------------------------------------
 
	IMPLICIT NONE
 
	! DATA DICTIONARY: Declared variables
	REAL    :: income = 0	! Income input by the user
	REAL    :: taxes = 0	! Tax owed based on input
 
	! Introduce the program
	WRITE (*,*)
	WRITE (*,*) '-----------------------------------------'
	WRITE (*,*)
	WRITE (*,*) '          Fortran Tax Calculator'
	WRITE (*,*)
	WRITE (*,*) '-----------------------------------------'
	WRITE (*,*)
	WRITE (*,*) 'This program calculates the taxes owed by'
	WRITE (*,*) 'a citizen of Australia, based on income.'
	WRITE (*,*)
 
	! Formats for output
1       FORMAT (1X,A2,1X,F10.2)
 
	! Get data values from the user
	WRITE (*,*) 'What was your income?'
	READ  (*,*) income
	WRITE (*,*)
 
	! Perform the income tax calculations
	IF      (income > 60000) THEN
		taxes = (15580 + (0.47 * (income - 60000)))
	ELSE IF (income > 50000) THEN
		taxes = (11380 + (0.42 * (income - 50000)))
	ELSE IF (income > 20000) THEN
		taxes = (2380 + (0.30 * (income - 20000)))
	ELSE IF (income > 6000) THEN
		taxes = (0.17 * (income - 6000))
	END IF
 
	! Add Medicare to the total
	taxes = taxes + (income * 0.015)
 
	! Return the results 
	WRITE (*,*) 'Your taxes are:'
	WRITE (*,1) 'A$', taxes
	WRITE (*,*)
	WRITE (*,*) '\\//_ Live long and prosper.'
	WRITE (*,*)
END PROGRAM tax

Using Einstein’s Relativity Formula to Find Mass Consumed by a Nuclear Reactor

Posted on February 28, 2012 by Jonathan Landrum
1

I was doing some reading about Fortran this afternoon and ran across this problem. The story goes that we have a 400-MW nuclear reactor that is 100% efficient in its conversion to energy of the Uranium isotope U-235. Obviously, no reactor is 100% efficient, but the purpose of this exercise is not to find mass consumed, but to practice developing algorithms under that premise.

The program uses Einstein’s mass-energy equivalence formula to find mass consumed by the reactor.

Einstein's mass-energy equivalence formula

After putting in the constants for energy produced (400,000,000 joules/second) and the speed of light squared (8.98755179 × 1016 meters/second), it is simply a matter of using the most basic Algebraic manipulation to isolate mass. I will admit, it was an astonishingly small amount of U-235 turned into energy. Perhaps there’s hope for this technology after all.

~Jonathan

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

PROGRAM relativity
! ------------------------------------------------------------------------------
! Programmer:     Jonathan Landrum
!
! Purpose:        This program calculates the amount of mass consumed by a
!                 nuclear power plant over a user-supplied length of time. The
!                 calculation uses Einstein's relativity formula, E = mc^2.
!
! Assumptions:    1.) The power station has a 400-MW (400,000,000 joules per
!                     second) nuclear reactor. This is our E constant.
!                 2.) The plant is 100% efficient (not realistic, but the losses
!                     are not part of the point of this exercise.)
! 
! Revisions:      Date		Programmer		Description of Change
!                 ===========	=====================	========================
!                 28 Feb 2012	Jonathan Landrum	Original code.
! ------------------------------------------------------------------------------
 
	IMPLICIT NONE
 
	! DATA DICTIONARY: Declared constants
	INTEGER, PARAMETER :: seconds = 31536000	! Seconds in a year
	INTEGER, PARAMETER :: E       = 400000000	! Output in joules/s
	REAL,    PARAMETER :: C       = 8.98755179E16	! Speed of light^2 in m/s
 
	! DATA DICTIONARY: Declared variables
	CHARACTER :: response		! User-supplied response to exit or not
	REAL      :: time		! User-supplied length of time in years
	REAL      :: mass		! Mass consumed by the generator in kg
 
	! Introduce the program
	WRITE (*,*) '* * * * * * * * * * * * * * * * * * * * * * * * * * * * *'
	WRITE (*,*) '*                                                       *'
	WRITE (*,*) '*       Fortran Nuclear Power Generator Simulator       *'
	WRITE (*,*) '*                                                       *'
	WRITE (*,*) '* * * * * * * * * * * * * * * * * * * * * * * * * * * * *'
	WRITE (*,*)
	WRITE (*,*) 'This program uses Einstein''s theory of relativity to'
	WRITE (*,*) 'simulate how much mass is consumed by a 400-MW nuclear'
	WRITE (*,*) 'reactor, given it is 100% efficient.'
	WRITE (*,*)
	WRITE (*,*) '---------------------------------------------------------'
	WRITE (*,*)
 
	! Prompt the user to continue or exit
	WRITE (*,*) 'Would you like to continue? [Y/N]'
	READ  (*,*) response
	WRITE (*,*)
 
	DO WHILE (response == 'y' .OR. response == 'Y')
 
		! Prompt the user for length of time to run the simulation
		WRITE (*,*) 'How many years do you want to simulate?'
		READ  (*,*) time
		WRITE (*,*)
 
		! Perform calculations
		mass = E / C * seconds * time
 
		! Return results
		WRITE (*,*) 'The amount of Uranium-235 turned into energy by'
		WRITE (*,*) 'this station in ', time, ' years is ', mass,' kg.'
		WRITE (*,*)
 
		! Prompt the user to continue or exit
		WRITE (*,*) 'Would you like to continue? [Y/N]'
		READ  (*,*) response
		WRITE (*,*)
 
	END DO ! End "continue" loop
	WRITE (*,*) '\\//_ Live long and prosper.'
 
! Exit
END PROGRAM relativity

Fortran Complex Arithmetic Calculator

Posted on February 27, 2012 by Jonathan Landrum
1

An assignment I was given last semester in Fortran asked us to make a program to calculate complex arithmetic problems. The point of the assignment was not complex arithmetic, but using external modules to solve a problem. We first had to make the module to include, called functions.f95, that contained all the complex arithmetic functions we would need. We were to call this module with a driver program. The Fortran compiler will link the two files together at compile time if we use include one specific line:

USE functions

the USE command pulls in the module named functions.f95 and allows the current program to use those functions as if they were contained in this program. That part was simple. Now on to the math.

There were five complex equations given in the assignment: addition, subtraction, multiplication, division, and exponentiation. The formulas for finding these answers were also supplied.

Complex Addition formula

Complex Addition formula

Complex Subtraction formula

Complex Subtraction formula

Complex Multiplication formula

Complex Multiplication formula

Complex Division formula

Complex Division formula

Complex Exponentiation formula

Complex Exponentiation formula

Now that we know the formulas, we can start making our functions.f95 file. The easiest way I saw to do it was to split up the Real and Imaginary parts of the formulas and handle them separately. Imaginary arithmetic operates just like Real arithmetic, with an i appended to the answer. Breaking this problem into pieces like this made a daunting task easily accomplished.

And, since I worked my way down the list, starting with addition, compiling, making sure it works, and then going on to subtraction, by the time I got to exponentiation I could simply call the multiplication functions I had already written to solve my problem. Here’s the finished functions.f95 file:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154

! =========================================
! Programmer:  Jonathan Landrum
! Date:        11 December 2011
! Class:       CSC 204
! =========================================
! Program:     functions.f95
! Purpose:     Provides functions for
!              complex arithmetic.
! Assumptions: None.
! =========================================
MODULE functions
CONTAINS
 
	! -------------------------------------
	!              Addition
	! -------------------------------------
 
	! -------------------------------------
	! addReal:
	! Performs real part of complex
	! addition
	! -------------------------------------
	INTEGER FUNCTION addReal(a, b, x, y)
		IMPLICIT NONE
		INTEGER :: a, b, x, y
		addReal = a + x
	END FUNCTION addReal
	! -------------------------------------
	! addImaginary:
	! Performs imaginary part of complex
	! addition
	! -------------------------------------
	INTEGER FUNCTION addImaginary(a, b, x, y)
		IMPLICIT NONE
		INTEGER :: a, b, x, y
		addImaginary = b + y
	END FUNCTION addImaginary
 
	! -------------------------------------
	!              Subtraction
	! -------------------------------------
 
	! -------------------------------------
	! subtractReal:
	! Performs real part of complex
	! subtraction
	! -------------------------------------
	INTEGER FUNCTION subtractReal (a, b, x, y)
		IMPLICIT NONE
		INTEGER :: a, b, x, y
		subtractReal = a - x
	END FUNCTION subtractReal
	! -------------------------------------
	! subtractImaginary:
	! Performs imaginary part of complex
	! subtraction
	! -------------------------------------
	INTEGER FUNCTION subtractImaginary (a, b, x, y)
		IMPLICIT NONE
		INTEGER :: a, b, x, y
		subtractImaginary = b - y
	END FUNCTION subtractImaginary
 
	! -------------------------------------
	!             Multiplication
	! -------------------------------------
 
	! -------------------------------------
	! multiplyReal:
	! Performs real part of complex
	! multiplication
	! -------------------------------------
	INTEGER FUNCTION multiplyReal (a, b, x, y)
		IMPLICIT NONE
		INTEGER :: a, b, x, y
		multiplyReal = a * x - b * y
	END FUNCTION multiplyReal
	! -------------------------------------
	! multiplyImaginary:
	! Performs imaginary part of complex
	! multiplication
	! -------------------------------------
	INTEGER FUNCTION multiplyImaginary (a, b, x, y)
		IMPLICIT NONE
		INTEGER :: a, b, x, y
		multiplyImaginary = a * y + b * x
	END FUNCTION multiplyImaginary
 
	! -------------------------------------
	!               Division
	! -------------------------------------
 
	! -------------------------------------
	! divideReal:
	! Performs real part of complex
	! division
	! -------------------------------------
	INTEGER FUNCTION divideReal (a, b, x, y)
		IMPLICIT NONE
		INTEGER :: a, b, x, y
		divideReal = a * x + b * y
	END FUNCTION divideReal
	! -------------------------------------
	! divideImaginary:
	! Performs imaginary part of complex
	! division
	! -------------------------------------
	INTEGER FUNCTION divideImaginary (a, b, x, y)
		IMPLICIT NONE
		INTEGER :: a, b, x, y
		divideImaginary = b * x - a * y
	END FUNCTION divideImaginary
 
	! -------------------------------------
	!            Exponentiation
	! -------------------------------------
 
	! -------------------------------------
	! exponentReal:
	! Performs complex exponentiation
	! -------------------------------------
	INTEGER FUNCTION exponentReal (a, b, n)
		IMPLICIT NONE
		INTEGER :: a, b, n, c = 1
		INTEGER :: out(2)
		out = (/a, b/)
		IF (n > 1) THEN
			WHILE (c < n) DO
				out(1) = multiplyReal(out(1), out(2), a, b)
                		out(2) = multiplyImaginary(out(1), out(2), a, b)
				c = c + 1
			END DO
		END IF
		exponentReal = out(1)
	END FUNCTION exponentReal
	! -------------------------------------
	! exponentImaginary:
	! Performs complex exponentiation
	! -------------------------------------
	INTEGER FUNCTION exponentImaginary (a, b, n)
		IMPLICIT NONE
		INTEGER :: a, b, n, c = 1
		INTEGER :: out(2)
		out = (/a, b/)
		IF (n > 1) THEN
			WHILE (c < n) DO
				out(1) = multiplyReal(out(1), out(2), a, b)
                		out(2) = multiplyImaginary(out(1), out(2), a, b)
				c = c + 1
			END DO
		END IF
		exponentImaginary = out(2)
	END FUNCTION exponentImaginary
END MODULE functions

Now that we’ve got our functions, we need a driver program that can use them. That’s easily accomplished, and I’ve even thrown in a bit of CLI design into this one. The main WHILE loop keeps the user in the loop unless they choose “X” to exit. They can select which type of complex arithmetic they’d like to do, and are then asked what the values are for a, b, x, and y. In the case of exponentiation, the user is asked for values for a and b, and a value for n, the degree to which the complex expression will be raised.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138

! ----------------------------------------------------------
! Programmer:   Jonathan Landrum
! Date:         11 December 2011
! Class:        CSC 204
! ----------------------------------------------------------
! Program:      driver.exe
! Purpose:      Performs complex arithmetic using modules.
! Assumptions:  Module "functions" is in the same directory.
! ----------------------------------------------------------
 
PROGRAM driver
	USE functions
	IMPLICIT NONE
 
	CHARACTER :: response = 'y', kind
	INTEGER   :: a, b, x, y, n
 
	WRITE (*,*) '*******************************************'
	WRITE (*,*) '*                                         *'
	WRITE (*,*) '*  FORTRAN COMPLEX ARITHMETIC CALCULATOR  *'
	WRITE (*,*) '*                                         *'
	WRITE (*,*) '*******************************************'
	WRITE (*,*)
	WRITE (*,*) 'PURPOSE:'
	WRITE (*,*) 'Performs complex arithmetic using modules.'
	WRITE (*,*)
	WRITE (*,*) 'Copyright (c) 2011 Jonathan Landrum'
	WRITE (*,*)
	WRITE (*,*) '-------------------------------------------'
	WRITE (*,*)
 
	! Any response beginning with the letter "y" will work.
	WHILE (response == 'y' .OR. response == 'Y') DO
 
		WRITE (*,*) 'What kind of arithmetic would you like to do?'
		WRITE (*,*)
		WRITE (*,*) 'Addition.............[A]'
		WRITE (*,*) 'Subtraction..........[S]'
		WRITE (*,*) 'Multiplication.......[M]'
		WRITE (*,*) 'Division.............[D]'
		WRITE (*,*) 'Exponentiation.......[E]'
		WRITE (*,*)
		WRITE (*,*) 'EXIT.................[X]'
		WRITE (*,*)
		READ  (*,*) kind
		WRITE (*,*)
 
		! Addition
		IF (kind == 'a' .OR. kind == 'A') THEN
			WRITE (*,*) 'Complex addition takes the form:'
			WRITE (*,*) '(a + bi) + (x + yi) = (a + x) + (b + y)i'
			WRITE (*,*) 'Enter a value for a:'
			READ  (*,*) a
			WRITE (*,*) 'Enter a value for b:'
			READ  (*,*) b
			WRITE (*,*) 'Enter a value for x:'
			READ  (*,*) x
			WRITE (*,*) 'Enter a value for y:'
			READ  (*,*) y
			WRITE (*,*) '(',addReal(a, b, x, y),' + ',addImaginary(a, b, x, y),'i)'
		END IF
 
		! Subtraction
		IF (kind == 's' .OR. kind == 'S') THEN
			WRITE (*,*) 'Complex subtraction takes the form:'
			WRITE (*,*) '(a + bi) - (x + yi) = (a - x) + (b - y)i'
			WRITE (*,*) 'Enter a value for a:'
			READ  (*,*) a
			WRITE (*,*) 'Enter a value for b:'
			READ  (*,*) b
			WRITE (*,*) 'Enter a value for x:'
			READ  (*,*) x
			WRITE (*,*) 'Enter a value for y:'
			READ  (*,*) y
			WRITE (*,*) '(',subtractReal(a, b, x, y),' - ',subtractImaginary(a, b, x, y),'i)'
		END IF
 
		! Multiplication
		IF (kind == 'm' .OR. kind == 'M') THEN
			WRITE (*,*) 'Complex multiplication takes the form:'
			WRITE (*,*) '(a + bi) * (x + yi) = (a * x - b * y) + (a * y + b * x)i'
			WRITE (*,*) 'Enter a value for a:'
			READ  (*,*) a
			WRITE (*,*) 'Enter a value for b:'
			READ  (*,*) b
			WRITE (*,*) 'Enter a value for x:'
			READ  (*,*) x
			WRITE (*,*) 'Enter a value for y:'
			READ  (*,*) y
			WRITE (*,*) '(',multiplyReal(a, b, x, y),' + ',multiplyImaginary(a, b, x, y),'i)'
		END IF
 
		! Division
		IF (kind == 'd' .OR. kind == 'D') THEN
			WRITE (*,*) 'Complex division takes the form:'
			WRITE (*,*) '(a + bi)   (a + bi) * (x - yi)'
			WRITE (*,*) '-------- = -------------------'
			WRITE (*,*) '(x + yi)        x^2 + y^2'
			WRITE (*,*) 'Enter a value for a:'
			READ  (*,*) a
			WRITE (*,*) 'Enter a value for b:'
			READ  (*,*) b
			WRITE (*,*) 'Enter a value for x:'
			READ  (*,*) x
			WRITE (*,*) 'Enter a value for y:'
			READ  (*,*) y
			WRITE (*,*) '(',divideReal(a, b, x, y),' * ',divideImaginary(a, b, x, y),'i)'
			WRITE (*,*) '------------------------------'
			WRITE (*,*) '     ',x * x + y * y
		END IF
 
		! Exponentiation
		IF (kind == 'e' .OR. kind == 'E') THEN
			WRITE (*,*) 'Complex exponentiation takes the form:'
			WRITE (*,*) '(a + bi)^n = (a + bi) * (a + bi) * (a + bi)...'
			WRITE (*,*) 'Enter a value for a:'
			READ  (*,*) a
			WRITE (*,*) 'Enter a value for b:'
			READ  (*,*) b
			WRITE (*,*) 'Enter a power to raise this complex number to:'
			READ  (*,*) n
			WRITE (*,*) '(',exponentReal(a, b, n),' + ',exponentImaginary(a, b, n),'i)'
		END IF
 
		! EXIT
		IF (kind == 'x' .OR. kind == 'X') THEN
			EXIT
		END IF
 
		WRITE (*,*)
		WRITE (*,*) 'Would you like to continue? [Y/N]'
		WRITE (*,*)
		READ  (*,*) response
		WRITE (*,*)
 
	END DO ! Continue loop
	WRITE (*,*) '\\//_ Live long and prosper.'
END PROGRAM driver

All in all, it took me a bit, but it was well worth it. I am very happy with how this one came out.

~Jonathan

Fortran Temperature Converter

Posted on February 24, 2012 by Jonathan Landrum
Reply

This is a fairly complex temperature converter I wrote last semester in Fortran after studying Chemistry over the summer. It’s not complex because it’s difficult, but because it protects the user in a lot of cases. The only thing I’ve found that will cause this program to fail is text input for the temperature to be converted.

However, there are plenty of loops to check for other user input. You’ll notice there are a couple “Unexpected Response” loops and a “Same Scale Selected” loop. The user is quite shielded from messing up. There are also logic loops for every conceivable case. You can convert to and from Celsius, Fahrenheit, and Kelvin scales, and the answer is presented in the correct format in every scenario. It won’t do your taxes for you, though.

~Jonathan

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210

! ---------------------------------------------------------
! Programmer:  Jonathan Landrum
! Date:        11 September 2011
! ---------------------------------------------------------
! Program:     tempConverter.f95
! Purpose:     Converts temperatures from one temperature
!              scale to another.
! Assumptions: 1.) The temperature input is in the form of
!                  a real number, not text. Text will cause
!                  the program to return an error.
!              2.) The temperature to be converted is in
!                  Celsius, Fahrenheit, or Kelvin scales.
! ---------------------------------------------------------
 
PROGRAM tempConverter
 
	IMPLICIT NONE
 
	! Variables:
	REAL              :: input = 0, output = 0
	CHARACTER         :: scaleIn, scaleOut, response
	CHARACTER(LEN=10) :: scaleInPretty, scaleOutPretty
 
	! Introduce the program and ask for initial response
	WRITE (*,*) '*****************************************'
	WRITE (*,*) '*                                       *'
	WRITE (*,*) '*     FORTRAN Temperature Converter     *'
	WRITE (*,*) '*                                       *'
	WRITE (*,*) '*****************************************'
	WRITE (*,*)
	WRITE (*,*) 'Copyright (c) 2011 Jonathan Landrum'
	WRITE (*,*)
	WRITE (*,*) 'This program takes a degree measurement'
	WRITE (*,*) 'from one scale, and converts it to the'
	WRITE (*,*) 'same value in another scale.'
	WRITE (*,*)
	WRITE (*,*) 'Scales are: Celsius, Fahrenheit, Kelvin'
	WRITE (*,*)
	WRITE (*,*) 'Would you like to continue? [Y/N]'
	WRITE (*,*)
 
	! Anything besides Y terminates
	READ  (*,*) response
 
	! -----------------------------------------------------
	! Main Loop:
	! Keeps the user in the program until they choose to
	! close it.
	! -----------------------------------------------------
	DO WHILE ((response == 'Y') .OR. (response == 'y'))
 
		! -------------------------------------------------
		! scaleIn Loop:
		! Ensures a legitimate scale has been selected.
		! -------------------------------------------------
		WRITE (*,*)
		WRITE (*,*) 'What scale is your temperature in?'
		WRITE (*,*) 'Celsius: [C], Fahrenheit: [F], Kelvin: [K]'
 
		READ  (*,*) scaleIn
 
		DO WHILE ((scaleIn /= 'c') .AND. &
		          (scaleIn /= 'C') .AND. &
		          (scaleIn /= 'f') .AND. &
		          (scaleIn /= 'F') .AND. &
		          (scaleIn /= 'k') .AND. &
		          (scaleIn /= 'K'))
 
			WRITE (*,*)
			WRITE (*,*) 'UNEXPECTED RESPONSE'
			WRITE (*,*)
			WRITE (*,*) 'What scale is your temperature in?'
			WRITE (*,*) 'Celsius: [C], Fahrenheit: [F], Kelvin: [K]'
 
			READ  (*,*) scaleIn
		END DO ! End scaleIn
 
		WRITE (*,*)
		WRITE (*,*) 'What is the value of your temperature?'
		WRITE (*,*) 'Use numbers. Text input will fail.'
 
		READ  (*,*) input
 
		! -------------------------------------------------
		! scaleOut Loop:
		! Ensures a legitimate scale has been selected.
		! -------------------------------------------------
		WRITE (*,*)
		WRITE (*,*) 'What scale do you want to convert to?'
		WRITE (*,*) 'Celsius: [C], Fahrenheit: [F], Kelvin: [K]'
 
		READ  (*,*) scaleOut
 
		DO WHILE ((scaleOut /= 'c') .AND. &
		          (scaleOut /= 'C') .AND. &
		          (scaleOut /= 'f') .AND. &
		          (scaleOut /= 'F') .AND. &
		          (scaleOut /= 'k') .AND. &
		          (scaleOut /= 'K'))
 
			WRITE (*,*)
			WRITE (*,*) 'UNEXPECTED RESPONSE'
			WRITE (*,*)
			WRITE (*,*) 'What scale do you want to convert to?'
			WRITE (*,*) 'Celsius: [C], Fahrenheit: [F], Kelvin: [K]'
 
			READ  (*,*) scaleOut
		END DO ! End scaleOut
 
		! -------------------------------------------------
		! Same Scale Loop:
		! Ensures the conversion between the same scale
		! won't happen. (No use wasting resources.)
		! -------------------------------------------------
		DO WHILE ((((scaleIn == 'c') .OR. (scaleIn == 'C')) .AND. ((scaleOut == 'c') .OR. (scaleOut == 'C')))  &
		     .OR. (((scaleIn == 'f') .OR. (scaleIn == 'F')) .AND. ((scaleOut == 'f') .OR. (scaleOut == 'F')))  &
		     .OR. (((scaleIn == 'k') .OR. (scaleIn == 'K')) .AND. ((scaleOut == 'k') .OR. (scaleOut == 'K'))))
 
			WRITE (*,*)
			WRITE (*,*) 'SAME SCALE SELECTED'
			WRITE (*,*)
			WRITE (*,*) 'What scale do you want to convert to?'
			WRITE (*,*) 'Celsius: [C], Fahrenheit: [F], Kelvin: [K]'
 
			READ  (*,*) scaleOut
		END DO ! End Same Scale
 
		! -------------------------------------------------
		! Main Conversion Loop:
		! Does the converting for every possible case, and
		! stores the output in a variable. This will be
		! printed later.
		! -------------------------------------------------
 
		! Celsius to Fahrenheit
		IF      (((scaleIn == 'c') .OR. (scaleIn == 'C')) .AND. ((scaleOut == 'f') .OR. (scaleOut == 'F'))) THEN
			output = (1.8 * input + 32)
 
		! Celsius to Kelvin
		ELSE IF (((scaleIn == 'c') .OR. (scaleIn == 'C')) .AND. ((scaleOut == 'k') .OR. (scaleOut == 'K'))) THEN
			output = (input + 273.15)
 
		! Fahrenheit to Celsius
		ELSE IF (((scaleIn == 'f') .OR. (scaleIn == 'F')) .AND. ((scaleOut == 'c') .OR. (scaleOut == 'C'))) THEN
			output = ((input - 32) / 1.8)
 
		! Fahrenheit to Kelvin
		ELSE IF (((scaleIn == 'f') .OR. (scaleIn == 'F')) .AND. ((scaleOut == 'k') .OR. (scaleOut == 'K'))) THEN
			output = ((input + 459.67) / 1.8)
 
		! Kelvin to Celsius
		ELSE IF (((scaleIn == 'k') .OR. (scaleIn == 'K')) .AND. ((scaleOut == 'c') .OR. (scaleOut == 'C'))) THEN
			output = (input - 273.15)
 
		! Kelvin to Fahrenheit
		ELSE IF (((scaleIn == 'k') .OR. (scaleIn == 'K')) .AND. ((scaleOut == 'f') .OR. (scaleOut == 'F'))) THEN
			output = (1.8 * input - 459.67)
 
		END IF ! End Main Conversion
 
		! -------------------------------------------------
		! Output Prettifier:
		! Converts all the c's and f's to their respective
		! scales. No k's because they can be hard-coded.
		! -------------------------------------------------
		IF      ((scaleIn == 'c') .OR. (scaleIn == 'C')) THEN
			scaleInPretty = 'Celsius'
		ELSE IF	((scaleIn == 'f') .OR. (scaleIn == 'F')) THEN
			scaleInPretty = 'Fahrenheit'
		END IF
 
		IF      ((scaleOut == 'c') .OR. (scaleOut == 'C')) THEN
			scaleOutPretty = 'Celsius'
		ELSE IF	((scaleOut == 'f') .OR. (scaleOut == 'F')) THEN
			scaleOutPretty = 'Fahrenheit'
		END IF
 
		! -------------------------------------------------
		! Output:
		! Return results of the conversion, and offer to
		! do another, or terminate the program.
		! -------------------------------------------------
 
		! Degrees in, degrees out
		IF      (((scaleIn == 'c') .OR. (scaleIn == 'C') .OR. (scaleIn == 'f') .OR. (scaleIn == 'F')) &
		  .AND. ((scaleOut == 'c') .OR. (scaleOut == 'C') .OR. (scaleOut == 'f') .OR. (scaleOut == 'F'))) THEN
			WRITE (*,*) input, ' degrees ', scaleInPretty, ' is equal to ', output, ' degrees ', scaleOutPretty
 
		! Degrees in, Kelvin out
		ELSE IF (((scaleIn == 'c') .OR. (scaleIn == 'C') .OR. (scaleIn == 'f') .OR. (scaleIn == 'F')) &
		  .AND. ((scaleOut == 'k') .OR. (scaleOut == 'K'))) THEN
			WRITE (*,*) input, ' degrees ', scaleInPretty, ' is equal to ', output, 'Kelvin'
 
		! Kelvin in, degrees out
		ELSE IF (((scaleIn == 'k') .OR. (scaleIn == 'K')) &
		  .AND. ((scaleOut == 'c') .OR. (scaleOut == 'C') .OR. (scaleOut == 'f') .OR. (scaleOut == 'F'))) THEN
			WRITE (*,*) input, ' Kelvin is equal to ', output, ' degrees ', scaleOutPretty
 
		! No case for Kelvin in, Kelvin out; same scale, won't compute
		END IF ! End Output
 
		WRITE (*,*)
		WRITE (*,*) 'Do you have another temperature to convert? [Y/N]'
		WRITE (*,*)
		READ  (*,*) response ! Anything besides Y terminates
	END DO ! End Main Program Loop
 
	WRITE (*,*)
	WRITE (*,*) '\\//_ Live long and prosper.'
END PROGRAM tempConverter

Calculating the Padovan Sequence with Fortran

Posted on February 17, 2012 by Jonathan Landrum
Reply

The Padovan Sequence is a recursive series I discovered recently, and it is similar to the Fibonacci Sequence I’ve written about before. The difference is the patterns the two recursive sequences make.

This is a graphic representation of the Fibonacci Sequence:

Graphic representation of the Fibonacci Sequence

And this is a graphic representation of the Padovan Sequence:

Graphic representation of the Padovan Sequence

The base case of the Padovan Sequence is defined as
f(0) = f(1) = f(2) = 1

and the general case is defined as
f(n) = f(n − 2) + f(n − 3)

So you can see that we can easily create an algorithm to return the Padovan Sequence using any language which supports recursion, and I have chosen Fortran to do so. We can do this with an if/else block, defining first the general case, and then the base caes (or vice versa, if you’re so inclined.)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

! ----------------------------------------------------------
! Programmer:   Jonathan Landrum
! Date:         17 February 2012
! ----------------------------------------------------------
! Program:      padovanSequence.f95
! Purpose:      Calculates the Padovan Sequence.
! Assumptions:  None.
! ----------------------------------------------------------
 
PROGRAM padovanSequence
 
	IMPLICIT NONE
 
	INTEGER :: counter = 0, input = 0, padovan
 
	WRITE (*,*) '* * * * * * * * * * * * * * * * * * * * * * *'
	WRITE (*,*) '*                                           *'
	WRITE (*,*) '*     FORTRAN PADOVAN SEQUENCE GENERATOR    *'
	WRITE (*,*) '*                                           *'
	WRITE (*,*) '* * * * * * * * * * * * * * * * * * * * * * *'
	WRITE (*,*)
	WRITE (*,*) 'This program calculates and returns members'
	WRITE (*,*) 'of the Padovan Sequence.'
	WRITE (*,*)
	WRITE (*,*) 'How many Padovan Numbers would like to return?'
	WRITE (*,*)
	READ  (*,*) input
 
	! If the user inputs a negative number, ask again
	DO WHILE (input < 0)
		WRITE (*,*)
		WRITE (*,*) 'ERROR: Negative input'
		WRITE (*,*)
		WRITE (*,*) 'How many Padovan Numbers would like to return?'
		WRITE (*,*)
		READ  (*,*) input
	END DO
	! END while negative
 
	WRITE (*,*)
	WRITE (*,*) '--------------------------------------------------'
 
	! We have a legitimate number, do the calculation
	DO WHILE (counter < input)
		WRITE (*,*) padovan(counter)
		counter = counter + 1
	END DO
	! END Processing loop
 
	WRITE (*,*) '--------------------------------------------------'
	WRITE (*,*)
	WRITE (*,*) '\\//_ Live long and prosper.'
	WRITE (*,*)
END PROGRAM padovanSequence
 
! ----------------------------------------------------------
! FUNCTION padovan:
! Calculates and returns the integer of the Padovan Sequence
! at the specified input position. The sequence is:
! 1, 1, 1, 2, 2, 3, 4, 5, 7, 9, 12, 16, 21, 28, 37, 49, etc.
! ----------------------------------------------------------
RECURSIVE FUNCTION padovan (input) RESULT (output)
	IMPLICIT NONE
 
	INTEGER :: input
	INTEGER :: output
 
	IF (input > 2) THEN
		output = padovan(input - 2) + padovan(input - 3)
	ELSE
		output = 1
	END IF
 
END FUNCTION padovan

Because this program uses recursion, caution should be observed when asking for more than 50 or 60 numbers. By putting the general case first in the padovan() function, I’ve saved a little time insofar that the if/else block only has one test now versus the three it would have if I had put the base case first.

A better way to get larger members of the set would be to either call for each member individually, or save previous output to an array or external file and have the function start from the previous number. Or you could even do as I did with my prime number calculator and have the user define a range. As it stands, the program re-calculates every number back to 1 each time, and it starts with 0 each time. So it is not useful for finding the higher members of the set. However, it was certainly a fun programming experience, as I hadn’t tried recursion with Fortran before.

~Jonathan

Using Fortran to Calculate the Sum of the Digits of a Very Large Number

Posted on February 11, 2012 by Jonathan Landrum
Reply

This one is a rather interesting conundrum. A website I found posits that the sum of the digits of 210 (1,024) is 7 (1 + 0 + 2 + 4 = 7). What is the sum of the digits of 21000? That number is exceedingly large. 302 digits long, to be exact. While it is nothing to add 302 digits, getting these digits can pose a problem. This number is well outside the limit of integer space, so using a language such as Fortran and returning 2**1000 won’t work. We need a better method, and as it turns out, this method is easier to work with, as well. We are going to use an external file and read the integers from it, line-by-line.

What we first need to do is hop over to Wolfram|Alpha and type in “2^1000” to get the numbers. After copying it to the clipboard, paste it into your text editor of choice, and separate each digit by a line feed. I put a couple of extra, empty lines at the bottom for good measure — don’t want to get a “Read past end of file” error. Name the file veryLarge.dat and put it in the same directory as the program. Now we’re ready to code.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

! ==================================================
! Programmer:  Jonathan Landrum
! Date:        10 February 2012
! ==================================================
! Program:     veryLarge.f95
! Purpose:     Adds the individual integers in the
!              number 2^1000. For example, 2^10 is
!              1,024. 1 + 0 + 2 + 4 = 7. For the
!              number 2^1000, there are 302 numbers
!              to add.
! Assumptions: 1.) The file "veryLarge.dat" is
!                  readable and is in the same
!                  directory as the program.
!              2.) The file "veryLarge.dat" is
!                  populated with the digits of the
!                  number 2^1000, each on its own
!                  line, and the last line is blank.
! ==================================================
 
PROGRAM veryLarge
 
	IMPLICIT NONE
 
	! ------------------------------------------
	! Variables
	! ------------------------------------------
	INTEGER :: counter = 0
	INTEGER :: number = 0
	INTEGER :: sum = 0
 
	! ------------------------------------------
	! Introduce the program
	! ------------------------------------------
	WRITE (*,*)
	WRITE (*,*) '--------------------------------------------------'
	WRITE (*,*) '-      FORTRAN Very Large Number Calculator      -'
	WRITE (*,*) '--------------------------------------------------'
	WRITE (*,*)
	WRITE (*,*) 'This program calculates the sum of each of the'
	WRITE (*,*) 'digits of the number 2^1000.'
	WRITE (*,*)
 
	OPEN (UNIT=8, FILE='veryLarge.dat')
 
	! ------------------------------------------
	! Main processing loop
	! ------------------------------------------
	DO WHILE (counter < 302)
		READ  (8,*) number
 
		sum = sum + number
		counter = counter + 1
	END DO
 
	CLOSE (8)
 
	! ------------------------------------------
	! Return the results
	! ------------------------------------------
	WRITE (*,*) sum
	WRITE (*,*) '--------------------------------------------------'
	WRITE (*,*) '\\//_ Live long and prosper.'
END PROGRAM veryLarge

The number returned is actually quite small, and you’ll get your answer in no time. This problem turned out to be more preparation than coding. And that happens sometimes. The only way to have made this easier would be to use a language that makes use of n-bit integers via bigint classes, so that you could return 2**1000 and use it, but then you would need to output that to an array, and that’s basically what you’ve done by hand here using the .dat file. But in this case, there wasn’t as much code to write. So, more-prep-less-code, or no-prep-more-code. Both get the answer in about the same time, and this way was easier.

~Jonathan

Fortran Prime Number Calculator

Posted on January 30, 2012 by Jonathan Landrum
2

This program is one of my favorites I’ve written so far. It takes input from the user for the lower and upper bounds of a range of integers, calculates which numbers in that range are prime numbers, and outputs a list of the primes and the total number of primes found.

I tested the output of this program against this list of primes, and it checked out perfectly. It’s a fairly fast program, returning the primes between 1 and 1,000,000 in about 6 seconds on my machine. Only assumptions are the input must be in Real form. Text input will cause the program to throw an exception.

~Jonathan

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

! ----------------------------------------------------------
! Programmer:   Jonathan Landrum
! Date:         21 September 2011
! Class:        CSC 204
! ----------------------------------------------------------
! Program:      primes.exe
! Purpose:      Calculates the number of prime numbers
!               there are in a given range, and returns
!               them in a list.
! Assumptions:  1.) Input from the user is in integer
!                   form. Text input will cause the
!                   program to fail. Real number input
!                   will be truncated into integers.
! ----------------------------------------------------------
 
PROGRAM primes
 
	IMPLICIT NONE
 
	INTEGER   :: total = 0, counter, lower, upper
	CHARACTER :: response
	LOGICAL   :: prime
 
	WRITE (*,*) '******************************************'
	WRITE (*,*) '*                                        *'
	WRITE (*,*) '*     FORTRAN PRIME NUMBER CALCULATOR    *'
	WRITE (*,*) '*                                        *'
	WRITE (*,*) '******************************************'
	WRITE (*,*)
	WRITE (*,*) 'Copyright (c) 2011 Jonathan Landrum'
	WRITE (*,*)
	WRITE (*,*) 'Calculates how many prime numbers are'
	WRITE (*,*) 'between the lower and upper bound given by'
	WRITE (*,*) 'the user, and returns the list of those'
	WRITE (*,*) 'prime numbers.'
	WRITE (*,*)
	WRITE (*,*) 'I wanted to use the Sieve of Eratosthenes'
	WRITE (*,*) 'algorithm, but could not figure out how to'
	WRITE (*,*) 'make it work. Perhaps later.'
	WRITE (*,*)
	WRITE (*,*) 'The output of this program has been checked'
	WRITE (*,*) 'against this list:'
	WRITE (*,*) 'http://primes.utm.edu/lists/small/100000.txt'
	WRITE (*,*)
	WRITE (*,*) '--------------------------------------------------'
	WRITE (*,*)
	WRITE (*,*) 'Would you like to continue? [Y/N]'
	WRITE (*,*)
	READ  (*,*) response
	WRITE (*,*)
 
	! Any response beginning with the letter "y" will work.
	WHILE (response == 'y' .OR. response == 'Y') DO
		WRITE (*,*)
		WRITE (*,*) 'What is the lower bound of numbers you'
		WRITE (*,*) 'would like to check?'
		WRITE (*,*)
		READ  (*,*) lower
		WRITE (*,*)
		WRITE (*,*) 'What is the upper bound of numbers you'
		WRITE (*,*) 'would like to check?'
		WRITE (*,*)
		READ  (*,*) upper
		WRITE (*,*)
 
		! If the user inputs a smaller number for the upper
		! bound, or if they input the same number, ask again
		WHILE ((upper - lower) <= 0) DO
			WRITE (*,*)
			WRITE (*,*) 'ERROR: Upper bound not higher than'
			WRITE (*,*) '       lower bound'
			WRITE (*,*)
			WRITE (*,*) 'What is the lower bound of numbers you'
			WRITE (*,*) 'would like to check?'
			WRITE (*,*)
			READ  (*,*) lower
			WRITE (*,*)
			WRITE (*,*) 'What is the upper bound of numbers you'
			WRITE (*,*) 'would like to check?'
			WRITE (*,*)
			READ  (*,*) upper
			WRITE (*,*)
		END DO ! END while upper is smaller
 
		! We have legitimate numbers
		counter = lower
		total = 0
		WHILE ((counter < upper).AND.(upper > 0)) DO
			IF (Prime(counter)) THEN
				WRITE (*, '(i0,1x)', ADVANCE = 'no') counter
				total = total + 1
			END IF ! Outputs the number
			counter = counter + 1
		END DO ! Processing loop
 
		WRITE (*,*)
		WRITE (*,*)
		WRITE (*,*) 'The number of primes between ',lower,' and ',upper,' is ',total,'.'
		WRITE (*,*)
		WRITE (*,*) '--------------------------------------------------'
		WRITE (*,*)
		WRITE (*,*) 'Would you like to continue? [Y/N]'
		WRITE (*,*)
		READ  (*,*) response
	END DO ! End Continue loop
	WRITE (*,*)
	WRITE (*,*) '\\//_ Live long and prosper.'
END PROGRAM primes
 
! ----------------------------------------------------------
! FUNCTION Prime:
! Determines if the input is a prime number. Returns T or F.
! ----------------------------------------------------------
LOGICAL FUNCTION Prime(n)
	IMPLICIT NONE
 
	INTEGER :: n, i
 
	IF (n <= 1) THEN                             ! 1 is not prime (Really, it's not)
		Prime = .FALSE.
	ELSE IF ((n == 2).OR.(n == 3)) THEN          ! Hardcode 2 and 3
		Prime = .TRUE.
	ELSE IF (MOD(n,2) == 0) THEN                 ! Get rid of evens
		Prime = .FALSE.
 
	! All other cases are out, so now we check to see if
	! n is divisible by the odd numbers from 3 on.
	ELSE
		i = 3
		Prime = .TRUE.                           ! Assume it's prime, then prove
		DO                                       ! If i^2 is not a root of n, or if n%i = 0
			IF (i*i > n .OR. MOD(n,i) == 0) EXIT !(Won't be larger than the square)
			i = i + 2                            ! Iterate the counter by 2 to get odds
		END DO
 
		! Record the answer
		IF (i*i > n) THEN
			Prime = .TRUE.
		ELSE
			Prime = .FALSE.
		END IF
	END IF ! We have our answer
END FUNCTION Prime