Точный расчет PDF стандартного нормального распределения с использованием C стандартной математической библиотеки

Question

Вероятность функции плотности стандартного нормального распределения, определяется как е ^{-x 2/2}/√ (2π). Это можно сделать простым способом в C-коде. Образец реализации одинарной точности может быть:

float my_normpdff (float a) 
{ 
    return 0x1.988454p-2f * my_expf (-0.5f * a * a); /* 1/sqrt(2*pi) */ 
} 

Хотя этот код свободен от преждевременного опустошения, существует проблема точности, так как ошибка понесены в вычислению /2 увеличивается при последующем экспоненциации , Это можно легко продемонстрировать с помощью тестов против высокоточных ссылок. Точная ошибка будет отличаться в зависимости от точности используемых применений exp() или expf(); для достоверно округленных функций возведения в степень обычно наблюдалась максимальная погрешность около 2 ulps для IEEE-754 binary32 одинарная точность, около 2 ulps для IEEE-754 binary64 двойная точность.

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

Answer 1

Вопрос точности поднят в вопросе может эффективно и эффективно, быть решены за счет использования ограниченных количеств двойного float или дважды double вычисления, облегченные использованием операции слиты умножения-сложения (ФМА).

Эта операция доступна, так как C99 через стандартные математические функции fmaf(a,b,c) и fma(a,b,c), которые вычисляют * B + C, без округления промежуточного продукта. Хотя функции непосредственно сопоставляются с быстрыми аппаратными операциями почти на всех современных процессорах, они могут использовать код эмуляции на более старых платформах, и в этом случае они могут быть очень медленными.

Это позволяет вычисление продукта с удвоенной точностью, используя только две операции, в результате чего в голове: хвост пара нативной точных цифр:

prod_hi = a * b   // head 
prod_lo = FMA (a, b, -hi) // tail 

высокого порядка бит результата может быть переданы в экспоненциацию, в то время как младшие биты используются для повышения точности результата с помощью линейной интерполяции, используя тот факт, что е ^х является его собственным производным:

e = exp (prod_hi) + exp (prod_hi) * prod_lo // exp (a*b) 

Это позволяет устранить большую часть ошибки наивной реализации. Другой, второстепенный, источник ошибки вычисления - ограниченная точность, с которой представлена ​​константа 1/√ (2π).Это может быть улучшено с помощью головки: хвост представления для константы, которая обеспечивает в два раза нативную точности и вычислительные:

r = FMA (const_hi, x, const_lo * x) // const * x 

В следующей статье указует на то, что этот метод может даже привести к правильно округленному умножению на некоторые произвольном константы точности:

Николас Брисебарре и Жан-Мишель Мюллер, «Правильно закругленное умножение на произвольные прецизионные константы», IEEE Transactions on Computers, Vol. . 57, № 2, февраль 2008, стр 165-174

Сочетание двух методов, и заботиться о нескольких случаях, связанных с угловыми пренебрежимо малых, мы приходим к следующему float реализации на основе IEEE-754 binary32:

float my_normpdff (float a) 
{ 
    const float RCP_SQRT_2PI_HI = 0x1.988454p-02f; /* 1/sqrt(2*pi), msbs */ 
    const float RCP_SQRT_2PI_LO = -0x1.857936p-27f; /* 1/sqrt(2*pi), lsbs */ 
    float ah, sh, sl, ea; 

    ah = -0.5f * a; 
    sh = a * ah; 
    sl = fmaf (a, ah, 0.5f * a * a); /* don't flip "sign bit" of NaN argument */ 
    ea = expf (sh); 
    if (ea != 0.0f) ea = fmaf (sl, ea, ea); /* avoid creation of NaN */ 
    return fmaf (RCP_SQRT_2PI_HI, ea, RCP_SQRT_2PI_LO * ea); 
} 

Соответствующий double реализации, на основе IEEE-754 binary64, выглядит почти идентично, для различных постоянных значений, используемых исключением:

double my_normpdf (double a) 
{ 
    const double RCP_SQRT_2PI_HI = 0x1.9884533d436510p-02; /* 1/sqrt(2*pi), msbs */ 
    const double RCP_SQRT_2PI_LO = -0x1.cbc0d30ebfd150p-56; /* 1/sqrt(2*pi), lsbs */ 
    double ah, sh, sl, ea; 

    ah = -0.5 * a; 
    sh = a * ah; 
    sl = fma (a, ah, 0.5 * a * a); /* don't flip "sign bit" of NaN argument */ 
    ea = exp (sh); 
    if (ea != 0.0) ea = fma (sl, ea, ea); /* avoid creation of NaN */ 
    return fma (RCP_SQRT_2PI_HI, ea, RCP_SQRT_2PI_LO * ea); 
} 

Точность этих реализаций зависит от точности стандартных математических функций expf() и exp(), соответственно. В тех случаях, когда математическая библиотека C предоставляет верные версии этих версий, максимальная ошибка одной из двух реализаций, описанных выше, обычно составляет менее 2,5 ulps.