Il y a quelque temps, j'avais écrit un code qui tentait de calculer sans utiliser les fonctions de bibliothèque. Hier, je passais en revue l'ancien code et j'ai essayé de le rendre aussi rapide que possible (et correct). Voici ma tentative jusqu'à présent:
const double ee = exp(1);
double series_ln_taylor(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 )
n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(1 - x) = -x - x**2/2 - x**3/3... */
n = 1 - n;
now = term = n;
for ( i = 1 ; ; ){
lgVal -= now;
term *= n;
now = term / ++i;
if ( now < 1e-17 ) break;
}
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Ici, j'essaie de trouver pour que soit un peu plus de n, puis j'ajoute la valeur de logarithme de , qui est inférieure à 1. À ce stade, l'extension Taylor de peut être utilisée sans souci.
J'ai récemment développé un intérêt pour l'analyse numérique, et c'est pourquoi je ne peux m'empêcher de poser la question, à quelle vitesse ce segment de code peut-il être exécuté dans la pratique, tout en étant suffisamment correct? Dois-je passer à d'autres méthodes, par exemple, en utilisant la fraction continue, comme celle-ci ?
La fonction fournie avec la bibliothèque standard C est presque 5,1 fois plus rapide que cette implémentation.
MISE À JOUR 1 : En utilisant la série arctan hyperbolique mentionnée dans Wikipedia , le calcul semble être presque 2,2 fois plus lent que la fonction de journal de bibliothèque standard C. Cependant, je n'ai pas vérifié de manière approfondie les performances, et pour un plus grand nombre, mon implémentation actuelle semble être VRAIMENT lente. Je veux vérifier à la fois ma mise en œuvre pour les erreurs liées et le temps moyen pour une large gamme de nombres si je peux gérer Voici mon deuxième effort.
double series_ln_arctanh(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
for ( i = 3 ; ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Toute suggestion ou critique est appréciée.
MISE À JOUR 2: Sur la base des suggestions ci-dessous, j'ai ajouté ici des modifications incrémentielles, qui sont environ 2,5 fois plus lentes que l'implémentation de bibliothèque standard. Cependant, je ne l'ai testé que pour les entiers cette fois, pour des nombres plus importants, le temps d'exécution augmenterait. Pour l'instant. Je ne connais pas encore de techniques pour générer des nombres doubles aléatoires , donc ce n'est pas encore complètement référencé. Pour rendre le code plus robuste, j'ai ajouté des corrections pour les cas d'angle. L'erreur moyenne pour les tests que j'ai faits est d'environ . ≤ 1 e 308 4 e - 15
double series_ln_better(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n == 0 ) return -1./0.; /* -inf */
if ( n < 0 ) return 0./0.; /* NaN*/
if ( n < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
/* the cutoff iteration is 650, as over e**650, term multiplication would
overflow. For larger numbers, the loop dominates the arctanh approximation
loop (with having 13-15 iterations on average for tested numbers so far */
for ( term = 1; term < n && lgVal < 650 ; term *= ee, lgVal++ );
if ( lgVal == 650 ){
n /= term;
for ( term = 1 ; term < n ; term *= ee, lgVal++ );
}
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
/* limiting the iteration for worst case scenario, maximum 24 iteration */
for ( i = 3 ; i < 50 ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
la source
frexpLa réponse de Kirill a déjà touché un grand nombre de questions pertinentes. Je voudrais développer certains d'entre eux sur la base d'une expérience pratique de conception de bibliothèque de mathématiques. Remarque: les concepteurs de bibliothèques mathématiques ont tendance à utiliser toutes les optimisations algorithmiques publiées, ainsi que de nombreuses optimisations spécifiques aux machines, qui ne seront pas toutes publiées. Le code sera fréquemment écrit en langage assembleur, plutôt que d'utiliser du code compilé. Il est donc peu probable qu'une implémentation simple et compilée atteigne plus de 75% des performances d'une implémentation de bibliothèque mathématique existante de haute qualité, en supposant des ensembles de fonctionnalités identiques (précision, gestion des cas spéciaux, rapport d'erreurs, prise en charge du mode d'arrondi).
La précision est généralement évaluée par comparaison avec une référence (tierce) de plus grande précision. Les fonctions à précision unique à argument unique peuvent facilement être testées de manière exhaustive, d'autres fonctions nécessitent des tests avec des vecteurs de test aléatoires (dirigés). De toute évidence, on ne peut pas calculer des résultats de référence infiniment précis, mais la recherche sur le dilemme du fabricant de tables suggère que pour de nombreuses fonctions simples, il suffit de calculer une référence avec une précision d'environ trois fois la précision cible. Voir par exemple:
Vincent Lefèvre, Jean-Michel Muller, "Pires cas pour un arrondi correct des fonctions élémentaires en double précision". In Proceedings 15th IEEE Symposium on Computer Arithmetic , 2001,111-118). (préimpression en ligne)
En termes de performances, il faut distinguer l'optimisation de la latence (importante quand on regarde le temps d'exécution des opérations dépendantes), et l'optimisation du débit (pertinent lorsqu'on considère le temps d'exécution des opérations indépendantes). Au cours des vingt dernières années, la prolifération des techniques de parallélisation matérielle telles que le parallélisme au niveau des instructions (par exemple superscalaire, processeurs hors service), le parallélisme au niveau des données (par exemple les instructions SIMD) et le parallélisme au niveau des threads (par exemple l'hyper-threading, processeurs multicœurs) a conduit à mettre l’accent sur le débit de calcul en tant que métrique la plus pertinente.
L'opération de fusion-ajout fusionnée ( FMA ), introduite pour la première fois par IBM il y a 25 ans, et maintenant disponible sur toutes les architectures de processeurs majeures, est un élément essentiel des implémentations de bibliothèques mathématiques modernes. Il offre une réduction des erreurs d'arrondi, offre une protection limitée contre l' annulation soustractive et simplifie considérablement l' arithmétique double-double .
C99log()C99fma()la source