Try porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 1)

Introduction

I've been living as an app development engineer for a long time, but when I wondered if I should try out something in the current season or study mathematics. I suddenly noticed. I originally had mathematics in my exam subject, and I should have done it at that university. Artificial intelligence, pattern recognition, etc. I've completely forgotten it. So, I first pulled this out from the various texts I had left. Numerical calculation programming

First edition May 1986. Numerical calculation programming is a book that is still a hit after 33 years. There aren't many books like this.

For the time being, I will try to port the code in this book to C and Python after reviewing. There are 32 in total, but I'm going to take it easy.

After this, the code that comes out is written in Visual Studio Code on Mac. Fortran and C are compiled using gcc.

Generation of machine epsilon

Rounding error

Summary of Chapter 1, P.4-5. The numbers in the program are basically rounded down or rounded. The error caused by the truncation or rounding is called ** rounding error **. The upper limit of the relative error that occurs when truncating to a β-ary n-digit floating point number is

\frac{\beta^{-n}}{\beta^{-1}} = \beta^{-(n-1)}

Maximum rounding error relative to this 1

\epsilon_M = \beta^{-(n-1)}

Is a value specific to the floating point system used and is called ** machine epsilon **. This value is

1\oplus\epsilon_M> 1　\tag{1}

Can also be defined as the smallest positive floating point number that satisfies. $ \ Oplus $ means to perform truncation or rounding operation after addition.

Fortran code

A Fortran program based on the definition in (1).

maceps.f

      SUBROUTINE MACEPS(EPSMAC)
        EPSMAC = 1.0
100    CONTINUE
       IF (1.0 + EPSMAC .LE. 1.0) THEN
        EPSMAC = 2 * EPSMAC
        RETURN
       END IF
       EPSMAC = EPSMAC * 0.5
       GO TO 100
      END 
    
      program maceps_main
        call maceps(epsmac)
        write(*,*)epsmac
      end program maceps_main

The code in the book was only the subroutine part, so the main part (bottom 4 lines) was added. The lowercase letters are the code I added, and the uppercase letters are the same as the book (Fortran is case insensitive).

The execution result is

   1.19209290E-07

Will be.

Understanding the code

Suddenly it's a GO TO statement from the first code. No, it can't be helped. The author (deceased) is a great teacher of numerical analysis, not an engineer. For the time being, let's wake it up in the flowchart.

Yup. For the time being, if a newcomer brings such a flow, it will rampage.

• Structure
• Don't leave the module in the middle of the module even though it is not inevitable

That's why I rewrote the flowchart.

In the basics of JIS flow chart, the loop condition is the end condition, but in C and Python to be ported from now on, the loop condition is the continuation condition, so in that form.

C code

maceps.c

#include <stdio.h>
#include <float.h>

void maceps(float *epsmac){
    *epsmac = 1.0f;
    while((1.0f + *epsmac) > 1.0f)
    {
        *epsmac *= 0.5;
    }
    *epsmac *=2;
}

int main(void){
    float epsmac;

    maceps(&epsmac);

    printf("%.8E\n", epsmac);
    printf("%.8E\n", FLT_EPSILON);

    return 0;
}

Since the Fortran code is single precision floating point, we also use float in C. I wanted to match it with the basic Fortran code, so I passed variables from the main with a pointer and displayed it in the main. 　 So, in C, the machine epsilon is defined by a constant in float.h, so I also displayed that (FLT_EPSILON) in the main.

The execution result is

1.19209290E-07
1.19209290E-07

Will be.

Python code

maceps.py

import numpy as np
def maceps(epsmac):
    epsmac[0] = 1.0
    epsmac[1] = 1.0
    while epsmac[1]+epsmac[0] > epsmac[1]:
        epsmac[0] = epsmac[0] * 0.5
    epsmac[0] = epsmac[0] * 2    
    return

epsmac = np.array([0,0],dtype=np.float32)

maceps(epsmac)

print("%.8E" % epsmac[0])
print("%.8E" % np.finfo(np.float32).eps)

Python standard floating point numbers are double precision, so use numpy to match single precision floating point. And since the machine epsilon is defined in numpy as well as C, that is also displayed.

The execution result is

1.19209290E-07
1.19209290E-07

is. This code, @konandoiruasa, pointed out in a comment that "if there is an operation with a constant, it will be float64." Thank you very much. Please refer to the remedy.

Summary

EPSMAC is divided in half during repetition (= multiplied by 0.5). The formula used for the repetition condition, 1.0 + EPSMAC > 1.0 If this doesn't hold, it means that the value before the last division is the minimum value that satisfies this, so stop repeating and double to cancel the last division. That is, the minimum positive value that satisfies the formula (1) = the difference between the minimum number greater than 1 and 1 = machine epsilon.