Calculating the Padovan Sequence with Fortran

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

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

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

Solving the Quadratic Equation with Fortran (With Complex Roots)

An interesting assignment given last semester was to write a program that returned the results of the Quadratic Equation to find the roots of a quadratic polynomial. What’s so interesting about that, you ask? The professor wanted the imaginary roots. In the end, it only took an “if” block to get it, but it did take some thought to get it right. At the time, this was one of my more proud moments as a programmer.

~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
! ---------------------------------------------------------
! Programmer:   Jonathan Landrum
! Date:         15 September 2011
! Class:        CSC 204
! ---------------------------------------------------------
! Program Name: quadraticEquation.exe
! Purpose:      Calculates the roots of a quadratic
!               equation, returning whether they are real
!               roots or imaginary roots.
! Assumptions:  1.) The input is in the form of a real
!                   number, not text. Text will cause the
!                   program to return an error.
! ---------------------------------------------------------
 
PROGRAM quadraticEquation
 
	IMPLICIT NONE
 
	! ------------------------------------------------------
	! Variables:
	! The response variable is a single character variable,
	! and input longer than one character will be truncated.
	! So technically, any word starting with the letter "y"
	! will keep the program running. Feature, not a bug(TM).
	! ------------------------------------------------------
	REAL		:: a, b, c, d = 0, fraction, radical
	CHARACTER	:: response
 
	! Introduce the program and ask for initial response
	WRITE (*,*) '*****************************************'
	WRITE (*,*) '*                                       *'
	WRITE (*,*) '* FORTRAN Quadratic Equation Calculator *'
	WRITE (*,*) '*                                       *'
	WRITE (*,*) '*****************************************'
	WRITE (*,*)
	WRITE (*,*) 'Copyright (c) 2011 Jonathan Landrum'
	WRITE (*,*)
	WRITE (*,*) 'This program calculates the roots of the'
	WRITE (*,*) 'quadratic equation, and returns whether'
	WRITE (*,*) 'the roots are real or imaginary.'
	WRITE (*,*)
	WRITE (*,*) 'The calculation is in the form:'
	WRITE (*,*) '(-B +/- (SQRT(B^2 - 4*A*C))'
	WRITE (*,*) '---------------------------'
	WRITE (*,*) '            2A'
	WRITE (*,*)
	WRITE (*,*) 'Do you want to continue? [Y/N]'
	WRITE (*,*)
 
	READ  (*,*) response ! Anything besides Y terminates
 
	! -----------------------------------------------------
	! Main Loop:
	! Keeps the user in the program until they choose to
	! close it.
	! -----------------------------------------------------
	WHILE ((response == 'Y') .OR. (response == 'y')) DO
 
		! -------------------------------------------------
		! Input Loop:
		! Takes in three real numbers for calculation.
		! -------------------------------------------------
		WRITE (*,*)
		WRITE (*,*) 'Enter a value for "A":'
		READ  (*,*) a
		WRITE (*,*) 'Enter a value for "B":'
		READ  (*,*) b
		WRITE (*,*) 'Enter a value for "C":'
		READ  (*,*) c
 
		d = b*b - 4*a*c
 
		! -------------------------------------------------
		! Computation Loop:
		! Determines if the discriminant is negative, and
		! computes the results.
		! -------------------------------------------------
		IF (d < 0) THEN
			fraction = -b/(2.0*a)
			radical  = (SQRT(ABS(d)))/(2.0*a)
			WRITE (*,*)
			WRITE (*,*) 'IMAGINARY ROOTS'
			WRITE (*,*)
			WRITE (*,*) fraction, ' + ', radical, 'i'
			WRITE (*,*) fraction, ' - ', radical, 'i'
		ELSE IF (d .EQ. 0) THEN
			WRITE (*,*)
			WRITE (*,*) 'REPEATED REAL ROOTS'
			WRITE (*,*)
			WRITE (*,*) -b/(2.0*a)
		ELSE
			WRITE (*,*)
			WRITE (*,*) 'REAL ROOTS'
			WRITE (*,*)
			WRITE (*,*) (-b + SQRT(d))/(2.0*a)
			WRITE (*,*) (-b - SQRT(d))/(2.0*a)
		END IF ! End computation Loop
 
		WRITE (*,*)
		WRITE (*,*) 'Do you want to convert another set of numbers? [Y/N]'
		WRITE (*,*)
		READ  (*,*) response ! Anything besides Y terminates
	END DO ! End Main Program Loop
 
	WRITE (*,*)
	WRITE (*,*) '\\//_ Live long and prosper.'
END PROGRAM quadraticEquation

Plotting Radioactive Decay with Fortran and ES-Plot

One of the assignments in Fortran last semester was to plot the exponential population growth of the United States between the years 1790 and 1860. Since I had recently taken Chemistry that summer, I also wrote a program to model the decay of the radioactive isotope Carbon-14, the isotope used in Radiocarbon Dating.

This program calculates the data, exports it to a CSV file, and finally calls the ES-Plot program and feeds it the CSV for data. An example output is below:

Carbon-14 Decay Model

In the graphic above, the vertical axis is number of moles present, and the horizontal axis is the number of years since the start of decay. It is difficult to read the numbers, but the first labeled tick mark is “1E+4”, 10,000 years. The average half-life of Carbon-14 is 5,730±40 years, not labeled here. The last tick mark is 40,000 years.

The user can input any quantity in moles for the program to calculate, but calculation will stop at 0.01% of the calculated mass because anything less than that is trivial in this case and is difficult to see in the resulting graphic. You will notice, therefore, that subsequent calculations result in the same graphic. That is because the decay rate of Carbon-14 is the same, regardless of quantity.

~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
! ----------------------------------------------------------
! Programmer:   Jonathan Landrum
! Date:         22 October 2011
! Class:        CSC 204
! ----------------------------------------------------------
! Program:      decay.f95
! Purpose:      Calculates and plots the radioactive decay
!               of Carbon-14. The half-life of Carbon-14 is
!               5,730 +/- 40 years. I have chosen 5,730 for
!               this model's lambda, as it is the mean.
! Assumptions:  1.) ES-Plot is installed.
!                   http://goo.gl/wmhZw
! ----------------------------------------------------------
 
PROGRAM decay
 
    IMPLICIT NONE
 
    CHARACTER                   :: response
    INTEGER                     :: t
    REAL                        :: p0, mass
    DOUBLE PRECISION, PARAMETER :: e = 2.71828182845904523536 !     ln 0.5
    DOUBLE PRECISION, PARAMETER :: k = -0.0001209680943385594 ! k = ------
                                                              !     lambda
 
    WRITE (*,*) '*******************************************'
    WRITE (*,*) '*                                         *'
    WRITE (*,*) '*FORTRAN CARBON RADIOACTIVE DECAY MODELLER*'
    WRITE (*,*) '*                                         *'
    WRITE (*,*) '*******************************************'
    WRITE (*,*)
    WRITE (*,*) 'Copyright (c) 2011 Jonathan Landrum'
    WRITE (*,*)
    WRITE (*,*) 'Calculates and plots the radioactive decay'
    WRITE (*,*) 'of Carbon-14 to the stable Nitrogen isotope'
    WRITE (*,*) 'Nitrogen-14, its most common isotope.'
    WRITE (*,*)
    WRITE (*,*) 'The half-life of Carbon-14 is 5,730 +/- 40'
    WRITE (*,*) 'years. I have chosen 5,730 for this model,'
    WRITE (*,*) 'as it is the mean.'
    WRITE (*,*)
    WRITE (*,*) 'Assumptions: 1.) ES-Plot is installed.'
    WRITE (*,*) '                 http://goo.gl/wmhZw'
    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  (*,*) 'How many moles of Carbon-14 do you want'
        WRITE  (*,*) 'to test?'
        WRITE  (*,*)
        READ   (*,*) p0
 
        t      = 0
        mass   = p0
 
        OPEN   (UNIT=8, FILE='decay.csv', status='unknown')
 
1       FORMAT (1X,I6,A2,F6.3)
 
        WRITE  (*,*) 'Year,   Mass'
        WRITE  (*,*) '------------'
        WRITE  (8,*) '$ Carbon-14 Decay Model'
 
        WHILE  (mass >= (p0 * 0.01)) DO
            mass  = (p0 * (e ** (k * t)))
            WRITE (8,*) t, ', ', mass
            WRITE (*,1) t, ', ', mass
            t     = t + 1
        END DO
 
        CLOSE  (8)
 
        WRITE  (*,*)
        WRITE  (*,*) '------------'
        WRITE  (*,*) 'End of list.'
        WRITE  (*,*)
        WRITE  (*,*) 'File "decay.csv" successfully created.'
        WRITE  (*,*) 'Calling ESPlot...'
 
        CALL SYSTEM ('"C:\Program Files (x86)\ESPlot\esplot.exe" decay.csv')
 
        WRITE  (*,*) 'Process Complete'
        WRITE  (*,*)
        WRITE  (*,*) '------------------------------------------'
        WRITE  (*,*) 'Would you calculate another mass? [Y/N]'
        WRITE  (*,*)
        READ   (*,*) response
        WRITE  (*,*)
    END DO ! Continue loop
    WRITE (*,*) '\\//_ Live long and prosper.'
END PROGRAM decay

Birthday Problem

In Probability and Stats last semester, my professor assigned us the project of determining how many students there would need to be in a class in order to have a 50/50 chance of having a repeat birth date. It’s called the Birthday Problem or Birthday Paradox. He mentioned that it’s quite a popular problem in probability classes, and that we cannot look online for the answer.

He wanted us to use our intuition to form a method for finding the answer. Obviously, since I’m a CS major, I turned to the computer for my answer. I was taking Fortran that same semester, so it was the freshest language on my mind. Any programming language would probably do, but I chose this one. A lot of what I’ve seen on the web with this problem is written in C or C++, so I feel I have a slight advantage here. ;~] Upon compiling and running the program, I came to my answer in less than one second: 23.

~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
! ---------------------------------------------------------
! Programmer:      Jonathan Landrum
! Date:            20 September 2011
! Class:           MAT 353
! ---------------------------------------------------------
! Program:         birthdays.exe
! Purpose:         Calculates the number of people you
!                  would need to question in order to have
!                  a 50% chance in getting two or more
!                  duplicate birth dates among the set.
! Assumptions:     None.
! ---------------------------------------------------------
 
PROGRAM BIRTHDAYS
	IMPLICIT NONE
 
	REAL      :: probability = 0.0, compliment = 1.0, days = 365
	INTEGER   :: counter = 1
	CHARACTER :: response
 
	WRITE (*,*) '******************************************'
	WRITE (*,*) '*                                        *'
	WRITE (*,*) '* FORTRAN DUPLICATE BIRTHDATE CALCULATOR *'
	WRITE (*,*) '*                                        *'
	WRITE (*,*) '******************************************'
	WRITE (*,*)
	WRITE (*,*) 'Calculates the number of people you would'
	WRITE (*,*) 'need to question in order to have a 50/50'
	WRITE (*,*) 'chance in having a duplicate birthdate.'
	WRITE (*,*)
	WRITE (*,*) 'This program takes no arguments from the'
	WRITE (*,*) 'user, except to continue after this line.'
	WRITE (*,*)
	WRITE (*,*) 'Would you like to continue? [Y/N]'
	WRITE (*,*)
	READ  (*,*) response
 
	! Any response beginning with the letter "y" will do.
	WHILE (response == 'y' .OR. response == 'Y') DO
 
		! -------------------------------------------------
		! Main:
		! This loop begins at 1 with the calculation and
		! continues until probability is greater than 50%.
		! -------------------------------------------------
		WHILE (probability < 0.51) DO
			compliment  = compliment * (days/365)
			probability = 1.0 - compliment
 
			WRITE (*,*) '| ', counter, ' | ', probability, ' |'
 
			counter = counter + 1
			days    = days - 1
		END DO ! Main loop
 
		WRITE (*,*)
		WRITE (*,*) 'Would you like to continue? [Y/N]'
		WRITE (*,*)
		READ  (*,*) response
	END DO ! Continue loop
	WRITE (*,*)
	WRITE (*,*) '\\//_ Live long and prosper.'
END PROGRAM BIRTHDAYS