IEEE-754 метод двойной точности и расщепления

Question

Когда вы вычисляете элементарные функции, вы применяете постоянную модификацию. Специально в реализации exp(x). Во всех этих реализациях любое исправление с помощью ln(2) выполняется в два этапа. ln(2) делится на два числа:

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

Тогда любое вычисление с помощью ln(2) выполняется с помощью:

 blablabla -= ln2p1
 blablabla -= ln2p2

Я знаю, это чтобы избежать эффекта округления. Но почему эти два числа специально? У некоторых из вас есть идея, как получить эти два числа?

Спасибо!

После первого комментария я дополняю этот пост большим количеством материала и очень странным вопросом. Я работал со своей командой, и мы согласны с тем, что потенциально можно удвоить точность, разделив число ln(2) на два числа. Для этого применяются два преобразования: первое:

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

k указывает на точность, во взглядах, подобных в библиотеке Cephes (~1980), для float k было установлено 9, 16 для double и также 16 для long long double (почему я не знаю). Таким образом, для двойного c_h имеет точность 16 бит, а для c_l 52 бита.

Исходя из этого, я пишу следующую программу и определяю c_h с точностью 52 бита.

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

}

Если я вычислю сейчас набор значений, я получу:

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

Это как когда x становится больше 300, появляется разница. Я посмотрел на реализацию gnulibc

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

в настоящее время он использует 16-битное предвидение для c_h (строка 84)

Ну, я, вероятно, что-то упустил, со стандартом IEEE, и я не могу представить ошибку точности в glibc. Как вы думаете?

Лучший,

Answer 1

ln2p1 точно 45426/65536. Это может быть получено round(65536 * ln(2)), ln2p2 это просто остаток. Итак, что же особенного в этих двух числах, так это знаменатель 65536 (2¹⁶).

Из того, что я обнаружил, большинство алгоритмов, использующих эту константу, можно проследить до библиотеки cephes, которая впервые была выпущена в 1984 году, когда 16-битные вычисления все еще доминировали, что, вероятно, объясняет, почему выбран 2¹⁶.