À propos d'une approximation plus rapide du journal (x)

10

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:log(X)

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.uneeuneneunelog(1 - X)

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.log(X)

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 - 151e81e3084e-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;
}
sarker306
la source

Réponses:

17

Ce n'est pas vraiment une réponse faisant autorité , plus une liste de problèmes que vous devriez considérer, et je n'ai pas testé votre code.

0. Comment avez-vous testé l'exactitude et la vitesse de votre code? Les deux sont importants, quelque peu difficiles à bien faire, et vous ne donnez pas de détails. En d'autres termes, si je compare votre fonction avec celle log de ma machine, est-ce que j'obtiendrai également les mêmes chiffres, , ? D'après mon expérience de la lecture des repères de synchronisation des autres dans la littérature universitaire, il faut beaucoup de soin et de précision pour obtenir des synchronisations reproductibles , qui sont les seules synchronisations dont on puisse se soucier. Les microbenchmarks en particulier sont notoirement peu fiables.5.12.15.1

1. Le problème courant de l'évaluation d'une fonction directement avec sa série de Taylor (non modifiée) est le nombre de termes nécessaires à la convergence. Il y a 52 bits dans la mantisse de a , donc quand au début de la boucle de la série Taylor, vous pouvez vous attendre à ce que la boucle prenne environ 50 itérations. C'est assez cher et devrait être optimisé.n 1f(x)doublen12

1.5. Avez-vous vérifié votre code pour le grand ? J'ai essayé , ce qui conduit à , puis à la boucle de la série Taylor, conduisant à une convergence extrêmement lente: elle converge comme la série Harmonic, c'est-à-dire qu'elle ne le fait pas mais il y aura au plus termes. En règle générale, vous devriez avoir une sorte de "max-iteration-count" lié à la boucle. Dans ce cas, il se comporte comme ceci car est fini, mais déborde vers dans la boucle de réduction des arguments. La bonne réponse est .10 17 n 709.78266108405500745 n1.7976e+308term=infn=1dix17nterm *= e709.78266108405500745

2. Les fonctions implémentées dans les bibliothèques standard devraient être extrêmement robustes. Renvoyer ( ) car le logarithme d'un nombre négatif (ou zéro) n'est pas correct. Le logarithme de doit être , le logarithme d'un nombre négatif doit être NaN.0 0 - dix-30000-

Je soupçonne qu'avec un petit effort vous pouvez sacrifier une partie de cette robustesse pour la performance, par exemple, en restreignant la plage d'arguments ou en renvoyant des résultats légèrement moins précis.

3. Les performances de ce type de code peuvent dépendre beaucoup de l'architecture du processeur sur laquelle il s'exécute. C'est un sujet profond et complexe, mais les fabricants de processeurs comme Intel publient des guides d'optimisation qui expliquent les différentes interactions entre votre code et le processeur sur lequel il fonctionne. La mise en cache peut être relativement simple, mais des choses comme la prédiction de branche, le parallélisme au niveau de l'instruction et les blocages de pipeline en raison des dépendances des données sont difficiles à voir précisément dans le code de haut niveau, mais comptent beaucoup pour les performances.

4. L' implémentation d'une fonction comme celle-ci signifie généralement que vous garantissez que pour le nombre à virgule flottante d'entrée , la sortie trouve à une certaine distance du point flottant le plus proche nombre à la vraie valeur . Vérifier que ce n'est pas totalement trivial, il n'y a aucune preuve dans votre code que vous l'avez fait, donc je ne sais pas si votre fonction est correcte (je suis sûr qu'elle est assez précise, mais à quel point ?). Ce n'est pas la même chose que de montrer que la série Taylor converge, en raison de la présence d'erreurs d'arrondi à virgule flottante.˜ y = ˜ f ( ˜ x )y=f( ˜ x )X~y~=F~(X~)y=F(X~)

4.5. Un bon moyen de tester la précision d'une fonction non testée serait de l'évaluer sur chacun des quatre milliards (moins si vous faites correctement la réduction d'argument, comme ici) de flotteurs simple précision, et de comparer les erreurs avec le journal standard de libm. Prend un peu de temps, mais au moins c'est complet.

5. Parce que vous connaissez dès le départ la précision des doubles, vous n'avez pas besoin d'avoir une boucle illimitée: le nombre d'itérations peut être calculé à l'avance (c'est probablement environ 50). Utilisez-le pour supprimer les branches de votre code ou au moins définir le nombre d'itérations à l'avance.

Toutes les idées habituelles sur le déroulement de la boucle s'appliquent également.

6. Il est possible d'utiliser des techniques d'approximation autres que la série Taylor. Il existe également des séries Chebyshev (avec la récurrence de Clenshaw), des approximants Pade et parfois des méthodes de recherche de racine comme la méthode de Newton chaque fois que votre fonction peut être refondue comme racine d'une fonction plus simple (par exemple, le fameux astuce sqrt ).

Les fractions continues ne seront probablement pas trop importantes, car elles impliquent une division, ce qui est beaucoup plus cher que les multiplications / ajouts. Si vous regardez _mm_div_ssà https://software.intel.com/sites/landingpage/IntrinsicsGuide/ , la division a la latence des cycles 13-14 et le débit de 5-14, en fonction de l'architecture, par rapport à 3-5 / 0,5-1 pour multiplier / ajouter / madd. Donc, en général (pas toujours), il est logique d'essayer d'éliminer autant que possible les divisions.

Malheureusement, les mathématiques ne sont pas comme un grand guide ici, parce que les expressions avec de courtes formules ne sont pas nécessairement les plus rapides les. Les mathématiques ne pénalisent pas les divisions, par exemple.

7. Les nombres à virgule flottante sont stockés en interne sous la forme (mantisse , , exposant ). Le journal naturel de est beaucoup moins naturel que le journal de base 2, pour lequel la première partie de votre code peut être remplacée par un appel à . m 1X=m×2em12<m1eXfrexp

8. Comparez votre logavec le login libmou openlibm(par exemple: https://github.com/JuliaLang/openlibm/blob/master/src/e_log.c ). C'est de loin le moyen le plus simple de découvrir ce que les autres ont déjà compris. Il existe également des versions spécialement optimisées de libm spécifiques aux fabricants de CPU, mais celles-ci n'ont généralement pas leur code source publié.

Boost :: sf a quelques fonctions spéciales, mais pas les fonctions de base. Il peut être instructif de regarder la source de log1p, cependant: http://www.boost.org/doc/libs/1_58_0/libs/math/doc/html/math_toolkit/powers/log1p.html

Il existe également des bibliothèques arithmétiques open-source de précision arbitraire comme mpfr, qui peuvent utiliser des algorithmes différents de libm en raison de la plus grande précision requise.

9. La précision et la stabilité des algorithmes numériques de Higham est une bonne introduction de niveau supérieur à l'analyse des erreurs des algorithmes numériques. Pour les algorithmes d'approximation eux-mêmes, la pratique d'approximation de la théorie de l'approximation par Trefethen est une bonne référence.

10. Je sais que cela est dit un peu trop souvent, mais les projets logiciels assez volumineux dépendent rarement du temps d'exécution d'une petite fonction appelée à plusieurs reprises. Ce n'est pas si courant d'avoir à se soucier des performances du journal, sauf si vous avez profilé votre programme et vous êtes assuré qu'il est important.

Kirill
la source
264-14e-15
1.13e-13term
 1e8
1
k=1dix7-1lnk
2
frexp X=m×2elnX=eln2+lnm
5

La 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).

eXplogerFcΓ

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.

log(1+X)=p(X)log(X)=2unetunenh((X-1)/(X+1))=p(((X-1)/(X+1))2)p

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()233

#include <math.h>

/* compute natural logarithm

   USE_ATANH == 1: maximum error found: 0.83482 ulp @ 0.7012829191167614
   USE_ATANH == 0: maximum error found: 0.83839 ulp @ 1.2788954397331760
*/
double my_log (double a)
{
    const double LOG2_HI = 0x1.62e42fefa39efp-01; // 6.9314718055994529e-01
    const double LOG2_LO = 0x1.abc9e3b39803fp-56; // 2.3190468138462996e-17
    double m, r, i, s, t, p, f, q;
    int e;

    m = frexp (a, &e);
    if (m < 0.70703125) { // 181/256
        m = m + m;
        e = e - 1;
    }
    i = (double)e;

    /* m in [181/256, 362/256] */

#if USE_ATANH
    /* Compute q = (m-1) / (m+1) */
    p = m + 1.0;
    m = m - 1.0;
    q = m / p;

    /* Compute (2*atanh(q)/q-2*q) as p(q**2), q in [-75/437, 53/309] */
    s = q * q;
    r =             0x1.2f1da230fb057p-3;  // 1.4800574027992994e-1
    r = fma (r, s,  0x1.399f73f934c01p-3); // 1.5313616375223663e-1
    r = fma (r, s,  0x1.7466542530accp-3); // 1.8183580149169243e-1
    r = fma (r, s,  0x1.c71c51a8bf129p-3); // 2.2222198291991305e-1
    r = fma (r, s,  0x1.249249425f140p-2); // 2.8571428744887228e-1
    r = fma (r, s,  0x1.999999997f6abp-2); // 3.9999999999404662e-1
    r = fma (r, s,  0x1.5555555555593p-1); // 6.6666666666667351e-1
    r = r * s;

    /* log(a) = 2*atanh(q) + i*log(2) = LOG2_LO*i + p(q**2)*q + 2q + LOG2_HI*i.
       Use K.C. Ng's trick to improve the accuracy of the computation, like so:
       p(q**2)*q + 2q = p(q**2)*q + q*t - t + m, where t = m**2/2.
    */
    t = m * m * 0.5;
    r = fma (q, t, fma (q, r, LOG2_LO * i)) - t + m;
    r = fma (LOG2_HI, i, r);

#else // USE_ATANH

    /* Compute f = m -1 */
    f = m - 1.0;
    s = f * f;

    /* Approximate log1p (f), f in [-75/256, 106/256] */
    r = fma (-0x1.961d64ddd82b6p-6, f, 0x1.d35fd598b1362p-5); // -2.4787281515616676e-2, 5.7052533321928292e-2
    t = fma (-0x1.fcf5138885121p-5, f, 0x1.b97114751d726p-5); // -6.2128580237329929e-2, 5.3886928516403906e-2
    r = fma (r, s, t);
    r = fma (r, f, -0x1.b5b505410388dp-5); // -5.3431043874398211e-2
    r = fma (r, f,  0x1.dd660c0bd22dap-5); //  5.8276198890387668e-2
    r = fma (r, f, -0x1.00bda5ecdad6fp-4); // -6.2680862565391612e-2
    r = fma (r, f,  0x1.1159b2e3bd0dap-4); //  6.6735934054864471e-2
    r = fma (r, f, -0x1.2489f14dd8883p-4); // -7.1420614809115476e-2
    r = fma (r, f,  0x1.3b0ee248a0ccfp-4); //  7.6918491287915489e-2
    r = fma (r, f, -0x1.55557d3b497c3p-4); // -8.3333481965921982e-2
    r = fma (r, f,  0x1.745d4666f7f48p-4); //  9.0909266480136641e-2
    r = fma (r, f, -0x1.999999d959743p-4); // -1.0000000092767629e-1
    r = fma (r, f,  0x1.c71c70bbce7c2p-4); //  1.1111110722131826e-1
    r = fma (r, f, -0x1.fffffffa61619p-4); // -1.2499999991822398e-1
    r = fma (r, f,  0x1.249249262c6cdp-3); //  1.4285714290377030e-1
    r = fma (r, f, -0x1.555555555f03cp-3); // -1.6666666666776730e-1
    r = fma (r, f,  0x1.999999999759ep-3); //  1.9999999999974433e-1
    r = fma (r, f, -0x1.fffffffffff53p-3); // -2.4999999999999520e-1
    r = fma (r, f,  0x1.555555555555dp-2); //  3.3333333333333376e-1
    r = fma (r, f, -0x1.0000000000000p-1); // -5.0000000000000000e-1

    /* log(a) = log1p (f) + i * log(2) */
    p = fma ( LOG2_HI, i, f);
    t = fma (-LOG2_HI, i, p);
    f = fma ( LOG2_LO, i, f - t);
    r = fma (r, s, f);
    r = r + p;
#endif // USE_ATANH

    /* Handle special cases */
    if (!((a > 0.0) && (a <= 0x1.fffffffffffffp1023))) {
        r = a + a;  // handle inputs of NaN, +Inf
        if (a  < 0.0) r =  0.0 / 0.0; //  NaN
        if (a == 0.0) r = -1.0 / 0.0; // -Inf
    }
    return r;
}
njuffa
la source
(+1) Savez-vous si les implémentations open-source courantes (comme openlibm) sont aussi bonnes qu'elles peuvent l'être, ou leurs fonctions spéciales peuvent-elles être améliorées?
Kirill
1
@Kirill Enfin, j'ai regardé les implémentations open source (il y a de nombreuses années), elles n'exploitaient pas les avantages de FMA. À l'époque, IBM Power et Intel Itanium étaient les seules architectures qui comprenaient l'opération, la prise en charge matérielle est désormais omniprésente. De plus, les approximations table plus polynôme étaient à la pointe de la technologie à l'époque, maintenant les tables sont en disgrâce: l'accès à la mémoire entraîne une consommation d'énergie plus élevée, ils peuvent (et font) interférer avec la vectorisation, et le débit de calcul a augmenté plus que le débit de la mémoire entraînant un impact potentiel négatif sur les performances des tableaux.
njuffa