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. maceps_flowchart1.png

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

That's why I rewrote the flowchart. maceps_flowchart2.png

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.

Recommended Posts

Try porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 1)
Try porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 3)
Try porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 2)
Try to solve the programming challenge book with python3
[Introduction to Python3 Day 1] Programming and Python
[Introduction to Python] I compared the naming conventions of C # and Python.
Try to write a program that abuses the program and sends 100 emails
Try to implement and understand the segment tree step by step (python)
Put Cabocha 0.68 on Windows and try to analyze the dependency with Python
C language to see and remember Part 2 Call C language from Python (argument) string
C language to see and remember Part 1 Call C language from Python (hello world)
Porting and modifying doublet-solver from python2 to python3.
Try to write python code to generate go code --Try porting JSON-to-Go and so on
Python amateurs try to summarize the list ①
C language to see and remember Part 4 Call C language from Python (argument) double
C language to see and remember Part 5 Call C language from Python (argument) Array
C language to see and remember Part 3 Call C language from Python (argument) c = a + b
How to use the C library in Python
How to generate permutations in Python and C ++
Try to solve the Python class inheritance problem
Try to solve the man-machine chart with Python
Try using the Python web framework Tornado Part 1
Try using the Python web framework Tornado Part 2
"Cython" tutorial to make Python explosive: Pass the C ++ class object to the class object on the Python side. Part 2
The story of porting code from C to Go and getting hooked (and to the language spec)
Think about how to program Python on the iPad
Try to make a Python module in C language
[Python] Try to read the cool answer to the FizzBuzz problem
Try to visualize the room with Raspberry Pi, part 1
Try to solve the internship assignment problem with Python
Try to operate DB with Python and visualize with d3
Try embedding Python in a C ++ program with pybind11
Tips and precautions when porting MATLAB programs to Python
Specifies the function to execute when the python program ends
I felt that I ported the Python code to C ++ 98.
How to enjoy Python on Android !! Programming on the go !!
Go language to see and remember Part 7 C language in GO language
[C / C ++] Pass the value calculated in C / C ++ to a python function to execute the process, and use that value in C / C ++.
"Cython" tutorial to make Python explosive: Pass a C ++ class object to a class object on the Python side. Part ①
How to start the PC at a fixed time every morning and execute the python program
Try to import to the database by manipulating ShapeFile of national land numerical information with Python
A story about porting the code of "Try and understand how Linux works" to Rust
[Python] A program to find the number of apples and oranges that can be harvested
[Introduction to cx_Oracle] (Part 6) DB and Python data type mapping
Try to make it using GUI and PyQt in Python
Try hitting the Twitter API quickly and easily with Python
Try to get the function list of Python> os package
I wanted to solve the Panasonic Programming Contest 2020 with Python
Try connecting to Supervisord via XMLRPC to start / stop the program
Measure the execution result of the program in C ++, Java, Python.
Write a python program to find the editing distance [python] [Levenshtein distance]
Consider If Programming Was An Anime in Python and C
[Part.2] Crawling with Python! Click the web page to move!
I tried to illustrate the time and time in C language
Try to bring up a subwindow with PyQt5 and Python
Build a Python environment and transfer data to the server
I tried programming the chi-square test in Python and Java.
Try to automate the operation of network devices with Python
Solving with Ruby and Python AtCoder ABC011 C Dynamic programming
I want to know the features of Python and pip
I tried to enumerate the differences between java and python