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.
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.
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.
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.
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.
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.
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.
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.
Recommended Posts