Essayez de porter le programme «Programmation informatique numérique FORTRAN77» vers C et Python (partie 1)

introduction

Je vis en tant qu'ingénieur en développement d'applications depuis longtemps, mais quand je me suis demandé si je devais essayer quelque chose dans la saison en cours ou étudier les mathématiques. J'ai soudain remarqué. Au départ, j'avais des mathématiques dans la matière que je suivais, alors j'aurais dû le faire dans cette université. Intelligence artificielle, reconnaissance de formes, etc. Je l'ai complètement oublié. Donc, j'ai d'abord tiré cela des différents textes que j'avais laissés. Programmation de calculs numériques

Première édition mai 1986. Si vous recherchez sur Google pour la programmation de calculs numériques, c'est toujours un succès après 33 ans. Il n'y a pas beaucoup de livres comme celui-ci.

Pour le moment, j'essaierai de porter le code de ce livre vers C et Python après avoir examiné. Il y en a 32 au total, mais je vais y aller doucement.

Après cela, le code qui sort est écrit dans Visual Studio Code sur Mac. Fortran et C sont compilés à l'aide de gcc.

Génération d'epsilon informatique

Erreur d'arrondi

Résumé du chapitre 1, P.4-5. Les nombres dans le programme sont fondamentalement tronqués ou arrondis. L'erreur causée par la troncature ou l'arrondi est appelée ** erreur d'arrondi **. La limite supérieure de l'erreur relative qui se produit lors de la troncature en virgule flottante β-aire à n chiffres est

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

Erreur d'arrondi maximale par rapport à ce 1

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

Est une valeur spécifique au système à virgule flottante utilisé, et s'appelle ** machine epsilon **. Cette valeur est

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

Peut également être défini comme le plus petit nombre à virgule flottante positive satisfaisant. $ \ Oplus $ signifie effectuer une opération de troncature ou d'arrondi après l'ajout.

Code Fortran

Programme Fortran basé sur la définition de (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

Le code dans le livre n'était que la partie sous-programme, donc la partie principale (4 lignes du bas) a été ajoutée. La minuscule est le code que j'ai ajouté et la majuscule est le livre (Fortran est insensible à la casse).

Le résultat de l'exécution est

   1.19209290E-07

Sera.

Comprendre le code

Soudain, c'est une instruction GO TO du premier code. Non, ça ne peut pas être aidé. L'auteur (décédé) est un grand professeur d'analyse numérique, pas un ingénieur. Pour le moment, réveillons-le dans l'organigramme.

Ouaip. Pour le moment, si un nouveau venu apporte un tel flux, il se déchaînera.

• Structure
• Ne laissez pas le module au milieu du module même si ce n'est pas une fatalité

C'est pourquoi j'ai réécrit l'organigramme.

Dans les bases du diagramme de flux JIS, la condition de boucle est la condition de fin, mais en C et Python à porter à partir de maintenant, la condition de boucle est la condition de continuation, donc sous cette forme.

Code C

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;
}

Puisque le code Fortran est à virgule flottante simple précision, nous utilisons également float en C. Je voulais le faire correspondre avec le code de base Fortran, j'ai donc passé des variables du main avec un pointeur et l'ai affiché dans le main. 　 Ainsi, en C, l'ordinateur epsilon est défini par une constante dans float.h, donc j'ai également affiché cela (FLT_EPSILON) dans le main.

Le résultat de l'exécution est

1.19209290E-07
1.19209290E-07

Sera.

Code Python

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)

Les nombres à virgule flottante standard Python sont à double précision, utilisez donc numpy pour faire correspondre la virgule flottante à simple précision. Et puisque l'ordinateur epsilon est défini dans numpy ainsi que C, cela est également affiché.

Le résultat de l'exécution est

1.19209290E-07
1.19209290E-07

est. Ce code, @konandoiruasa, a souligné dans un commentaire que "s'il y a une opération avec une constante, ce sera float64." Merci beaucoup. Veuillez vous référer au remède.

Résumé

EPSMAC est divisé en deux lors de la répétition (= multiplié par 0,5). La formule utilisée pour la condition de répétition, 1.0 + EPSMAC > 1.0 Si cela ne tient pas, cela signifie que la valeur avant la dernière division est la valeur minimale qui satisfait cela, alors arrêtez de répéter et doublez pour annuler la dernière division. Autrement dit, la valeur positive minimale qui satisfait la formule (1) = la différence entre le nombre minimal supérieur à 1 et 1 = epsilon de l'ordinateur.