Méthode de double précision et de division IEEE-754

Question

Lorsque vous calculez des fonctions élémentaires, vous appliquez une modification constante. Spécialement dans l'implémentation de exp (x). Dans toutes ces implémentations, toute correction avec ln (2) se fait en deux étapes. Dans (2) est divisé en deux chiffres:

static const double ln2p1 = 0.693145751953125; 
static const double ln2p2 = 1.42860682030941723212E-6; 
// then ln(2) = ln2p1 + ln2p2 

Ensuite, tout calcul avec ln (2) se fait par:

blablabla -= ln2p1 
blablabla -= ln2p2 

Je sais qu'il est d'éviter l'effet arrondi. Mais pourquoi ce deux chiffres en particulier? Certains d'entre vous ont une idée de comment obtenir ces deux chiffres?

Merci!

Suite au premier commentaire, je complète ce post avec plus de matériel et une question très étrange. J'ai travaillé avec mon équipe et nous sommes d'accord que l'accord est de doubler potentiellement la précision en divisant le nombre ln (2) en deux nombres. Pour cela, deux transformations sont appliquées, le premier:

1) c_h = floor(2^k ln(2))/2^k 
2) c_l = ln(2) - c_h 

k indique les precisions, en regard aime dans la bibliothèque Cephes (~ 1980) pour flotteur k a été fixé sur 9, 16 pour le double et aussi 16 pour longtemps long double (pourquoi je ne sais pas). Donc pour c_h double a une précision de 16 bits mais 52 bits pour c_l. À partir de là, j'écris le programme suivant et détermine c_h avec une précision de 52 bits.

#include <iostream> 
#include <math.h> 
#include <iomanip> 

enum precision { nine = 9, sixteen = 16, fiftytwo = 52 }; 

int64_t k_helper(double x){ 
    return floor(x/log(2)); 
} 

template<class C> 
double z_helper(double x, int64_t k){ 
    x -= k*C::c_h; 
    x -= k*C::c_l; 
    return x; 
} 

template<precision p> 
struct coeff{}; 

template<> 
struct coeff<nine>{ 
    constexpr const static double c_h = 0.693359375; 
    constexpr const static double c_l = -2.12194440e-4; 
}; 

template<> 
struct coeff<sixteen>{ 
    constexpr const static double c_h = 6.93145751953125E-1; 
    constexpr const static double c_l = 1.42860682030941723212E-6; 
}; 

template<> 
struct coeff<fiftytwo>{ 
    constexpr const static double c_h = 0.6931471805599453972490664455108344554901123046875; 
    constexpr const static double c_l = -8.78318343240526578874146121703272447458793199905066E-17; 
}; 


int main(int argc, const char * argv[]) { 

    double x = atof(argv[1]); 
    int64_t k = k_helper(x); 

    double z_9 = z_helper<coeff<nine> >(x,k); 
    double z_16 = z_helper<coeff<sixteen> >(x,k); 
    double z_52 = z_helper<coeff<fiftytwo> >(x,k); 


    std::cout << std::setprecision(16) << " 9 bits precisions " << z_9 << "\n" 
             << " 16 bits precisions " << z_16 << "\n" 
             << " 52 bits precisions " << z_52 << "\n"; 



    return 0; 

}

Si je calcule maintenant un ensemble de valeurs différentes que je reçois:

bash-3.2$ g++ -std=c++11 main.cpp 
bash-3.2$ ./a.out 1 
9 bits precisions 0.30685281944 
16 bits precisions 0.3068528194400547 
52 bits precisions 0.3068528194400547 
bash-3.2$ ./a.out 2 
9 bits precisions 0.61370563888 
16 bits precisions 0.6137056388801094 
52 bits precisions 0.6137056388801094 
bash-3.2$ ./a.out 100 
9 bits precisions 0.18680599936 
16 bits precisions 0.1868059993678755 
52 bits precisions 0.1868059993678755 
bash-3.2$ ./a.out 200 
9 bits precisions 0.37361199872 
16 bits precisions 0.3736119987357509 
52 bits precisions 0.3736119987357509 
bash-3.2$ ./a.out 300 
9 bits precisions 0.56041799808 
16 bits precisions 0.5604179981036264 
52 bits precisions 0.5604179981036548 
bash-3.2$ ./a.out 400 
9 bits precisions 0.05407681688 
16 bits precisions 0.05407681691155647 
52 bits precisions 0.05407681691155469 
bash-3.2$ ./a.out 500 
9 bits precisions 0.24088281624 
16 bits precisions 0.2408828162794319 
52 bits precisions 0.2408828162794586 
bash-3.2$ ./a.out 600 
9 bits precisions 0.4276888156 
16 bits precisions 0.4276888156473074 
52 bits precisions 0.4276888156473056 
bash-3.2$ ./a.out 700 
9 bits precisions 0.61449481496 
16 bits precisions 0.6144948150151828 
52 bits precisions 0.6144948150151526 

Il aime quand x devient supérieure à 300 apparaît une différence. J'ai eu un regard sur la mise en œuvre de gnulibc

http://osxr.org:8080/glibc/source/sysdeps/ieee754/ldbl-128/s_expm1l.c

actuellement il utilise les 16 bits prevision pour C_h (ligne 84)

Eh bien, je suis probablement manque quelque chose, avec la norme IEEE, et Je ne peux pas imaginer une erreur de précision dans la glibc. Qu'est-ce que tu penses ?

Best,

Answer 1

ln2p1 est exactement 45426/65536. Ceci peut être obtenu par round(65536 * ln(2)). ln2p2 est simplement le reste. Donc, ce qui est si spécial à propos de ces deux nombres est le dénominateur 65536 (2).

D'après ce que je trouve la plupart des algorithmes utilisant cette constante remonte à la bibliothèque cephes, qui a été d'abord publié en 1984, où l'informatique 16 bits était encore dominant, ce qui explique probablement pourquoi 2 est choisi.