Versuchen Sie, das Programm "FORTRAN77 Numerical Computing Programming" auf C und Python zu portieren (Teil 1).

Einführung

Ich habe lange Zeit als App-Entwicklungsingenieur gelebt, aber als ich mich fragte, ob ich in der aktuellen Saison etwas ausprobieren oder Mathematik studieren sollte. Ich bemerkte es plötzlich. Ich hatte ursprünglich Mathematik in dem Fach, das ich belegte, also hätte ich es an dieser Universität machen sollen. Künstliche Intelligenz, Mustererkennung usw. Ich habe es völlig vergessen. Also zog ich dies zuerst aus den verschiedenen Texten heraus, die ich hinterlassen hatte. Programmierung numerischer Berechnungen

Erstausgabe Mai 1986. Als ich die numerische Berechnungsprogrammierung googelte, war es ein Buch, das nach 33 Jahren immer noch ein Hit war. Es gibt nicht viele solche Bücher.

Vorerst werde ich versuchen, den Code in diesem Buch nach der Überprüfung auf C und Python zu portieren. Insgesamt gibt es 32, aber ich werde es ruhig angehen lassen.

Danach wird der ausgegebene Code in Visual Studio Code auf dem Mac geschrieben. Fortran und C werden mit gcc kompiliert.

Erzeugung von Computer-Epsilon

Rundungsfehler

Zusammenfassung von Kapitel 1, S. 4-5. Die Zahlen im Programm sind grundsätzlich abgeschnitten oder gerundet. Der durch das Abschneiden oder Runden verursachte Fehler wird als ** Rundungsfehler ** bezeichnet. Die Obergrenze des relativen Fehlers, der beim Abschneiden auf einen β-ary n-stelligen Gleitkomma auftritt, ist

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

Maximaler Rundungsfehler relativ zu diesem 1

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

Ist ein Wert, der für das verwendete Gleitkommasystem spezifisch ist und als ** Maschinen-Epsilon ** bezeichnet wird. Dieser Wert ist

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

Kann auch als die kleinste positive Gleitkommazahl definiert werden, die erfüllt. $ \ Oplus $ bedeutet, dass nach dem Hinzufügen eine Kürzungs- oder Rundungsoperation ausgeführt wird.

Fortran-Code

Fortran-Programm basierend auf der Definition von (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

Der Code im Buch war nur der Unterprogrammteil, daher wurde der Hauptteil (untere 4 Zeilen) hinzugefügt. Der Kleinbuchstabe ist der Code, den ich hinzugefügt habe, und der Großbuchstabe ist das Buch (Fortran unterscheidet nicht zwischen Groß- und Kleinschreibung).

Das Ausführungsergebnis ist

   1.19209290E-07

Wird sein.

Den Code verstehen

Plötzlich ist es eine GO TO-Anweisung aus dem ersten Code. Nein, es kann nicht geholfen werden. Der Autor (verstorben) ist ein großartiger Lehrer für numerische Analyse, kein Ingenieur. Lassen Sie es uns vorerst im Flussdiagramm wecken.

Jep. Wenn ein Neuling einen solchen Fluss bringt, wird er vorerst toben.

• Struktur
• Lassen Sie das Modul nicht in der Mitte des Moduls, auch wenn dies nicht unvermeidlich ist

Deshalb habe ich das Flussdiagramm umgeschrieben.

In den Grundlagen des JIS-Flussdiagramms ist die Schleifenbedingung die Endbedingung, aber in C und Python, die von nun an portiert werden sollen, ist die Schleifenbedingung die Fortsetzungsbedingung, also in dieser 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;
}

Da der Fortran-Code ein Gleitkomma mit einfacher Genauigkeit ist, verwenden wir auch Gleitkomma in C. Ich wollte es mit dem grundlegenden Fortran-Code abgleichen, also habe ich Variablen vom Main mit einem Zeiger übergeben und im Main angezeigt. 　 In C wird das Computer-Epsilon durch eine Konstante in float.h definiert, daher habe ich auch das (FLT_EPSILON) im Main angezeigt.

Das Ausführungsergebnis ist

1.19209290E-07
1.19209290E-07

Wird sein.

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-Gleitkommazahlen haben eine doppelte Genauigkeit. Verwenden Sie daher numpy, um Gleitkommazahlen mit einfacher Genauigkeit zu entsprechen. Und da Computer-Epsilon sowohl in Numpy als auch in C definiert ist, wird dies auch angezeigt.

Das Ausführungsergebnis ist

1.19209290E-07
1.19209290E-07

ist. Dieser Code, @konandoiruasa, wies in einem Kommentar darauf hin, dass "wenn es eine Operation mit einer Konstanten gibt, es float64 sein wird". Vielen Dank. Bitte beziehen Sie sich auf das Mittel.

Zusammenfassung

EPSMAC wird während der Wiederholung in zwei Hälften geteilt (= multipliziert mit 0,5). Die für die Wiederholungsbedingung verwendete Formel: 1.0 + EPSMAC > 1.0 Wenn dies nicht zutrifft, ist der Wert vor der letzten Division der Mindestwert, der dies erfüllt. Hören Sie also auf zu wiederholen und verdoppeln Sie, um die letzte Division abzubrechen. Das heißt, der minimale positive Wert, der die Formel (1) erfüllt = die Differenz zwischen der minimalen Anzahl größer als 1 und 1 = Computer-Epsilon.