용지 간단 다항식 의무가 단단히 경계 헬퍼 기능 스케일링 상보 오차 함수를 매핑 영리한 변환 제안 근사. 성능 향상을 위해 다양한 변형을 실험했지만,이 모두가 정확도에 부정적인 영향을 미쳤습니다. 변환 (x - K)/(x + K)에서의 상수 K의 선택은 핵심 근사의 정확도와 명확하지 않은 관계를 갖는다. 나는 경험적으로 종이와 다른 "최적"값을 결정했다.
핵심 근사값에 대한 인수와 중간 결과를 다시 erfcx으로 변환하면 추가 반올림 오류가 발생합니다. 정확성에 미치는 영향을 줄이려면 보상 단계를 적용해야합니다. 이전에 question & answer regarding erfcf에 자세하게 설명되어 있습니다.FMA의 가용성은이 작업을 크게 단순화합니다.
/*
* Based on: M. M. Shepherd and J. G. Laframboise, "Chebyshev Approximation of
* (1+2x)exp(x^2)erfc x in 0 <= x < INF", Mathematics of Computation, Vol. 36,
* No. 153, January 1981, pp. 249-253.
*
*/
float my_erfcxf (float x)
{
float a, d, e, m, p, q, r, s, t;
a = fmaxf (x, 0.0f - x); // NaN-preserving absolute value computation
/* Compute q = (a-2)/(a+2) accurately. [0,INF) -> [-1,1] */
m = a - 2.0f;
p = a + 2.0f;
#if FAST_RCP_SSE
r = fast_recipf_sse (p);
#else
r = 1.0f/p;
#endif
q = m * r;
t = fmaf (q + 1.0f, -2.0f, a);
e = fmaf (q, -a, t);
q = fmaf (r, e, q);
/* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1,1] */
p = 0x1.f10000p-15f; // 5.92470169e-5
p = fmaf (p, q, 0x1.521cc6p-13f); // 1.61224554e-4
p = fmaf (p, q, -0x1.6b4ffep-12f); // -3.46481771e-4
p = fmaf (p, q, -0x1.6e2a7cp-10f); // -1.39681227e-3
p = fmaf (p, q, 0x1.3c1d7ep-10f); // 1.20588380e-3
p = fmaf (p, q, 0x1.1cc236p-07f); // 8.69014394e-3
p = fmaf (p, q, -0x1.069940p-07f); // -8.01387429e-3
p = fmaf (p, q, -0x1.bc1b6cp-05f); // -5.42122945e-2
p = fmaf (p, q, 0x1.4ff8acp-03f); // 1.64048523e-1
p = fmaf (p, q, -0x1.54081ap-03f); // -1.66031078e-1
p = fmaf (p, q, -0x1.7bf5cep-04f); // -9.27637145e-2
p = fmaf (p, q, 0x1.1ba03ap-02f); // 2.76978403e-1
/* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
d = a + 0.5f;
#if FAST_RCP_SSE
r = fast_recipf_sse (d);
#else
r = 1.0f/d;
#endif
r = r * 0.5f;
q = fmaf (p, r, r); // q = (p+1)/(1+2*a)
t = q + q;
e = (p - q) + fmaf (t, -a, 1.0f); // residual: (p+1)-q*(1+2*a)
r = fmaf (e, r, q);
if (a > 0x1.fffffep127f) r = 0.0f; // 3.40282347e+38 // handle INF argument
/* Handle negative arguments: erfcx(x) = 2*exp(x*x) - erfcx(|x|) */
if (x < 0.0f) {
s = x * x;
d = fmaf (x, x, -s);
e = expf (s);
r = e - r;
r = fmaf (e, d + d, r);
r = r + e;
if (e > 0x1.fffffep127f) r = e; // 3.40282347e+38 // avoid creating NaN
}
return r;
}
부정적인 반면에서이 구현의 최대 오차는 expf()의 표준 수학 라이브러리의 구현의 정확성에 따라 달라집니다 다음과 같이
결과 단 정밀도 코드가 보인다. 인텔 컴파일러 버전 13.1.3.198을 사용하고 /fp:strict으로 컴파일하면 철저한 테스트에서 양수 반면에서 최대 2.00450 ulps 및 음수 반면에서 2.38412 ulps의 최대 오류를 관찰했습니다. 이 시점에서 내가 말할 수있는 가장 좋은 점은 expf()의 충실하게 둥근 구현은 최대 ulper가 2.5 ulp 미만이 될 것입니다.
코드가 두 작업 영역을 포함하고 있지만 작업 속도가 느릴 수 있지만 특수 형식의 역수대로 발생하므로 많은 플랫폼에서 빠른 상호 근사를 사용할 수 있습니다. 상호 근사가 충실하게 반올림되는 한, 실험에 기초하여 erfcxf() 정확도에 미치는 영향은 무시해도 될 것 같습니다. 고속 SSE 버전 (최대 오류 < 2.0 ulps)과 같이 약간 큰 오류조차도 사소한 영향을 미치는 것으로 보입니다.
/* Fast reciprocal approximation. HW approximation plus Newton iteration */
float fast_recipf_sse (float a)
{
__m128 t;
float e, r;
t = _mm_set_ss (a);
t = _mm_rcp_ss (t);
_mm_store_ss (&r, t);
e = fmaf (0.0f - a, r, 1.0f);
r = fmaf (e, r, r);
return r;
}
두 번 정밀 버전 erfcx()는 단 정밀도 버전 erfcxf() 구조적으로 동일하지만 더 많은 조건을 가진 최소 최대 다항식 근사를 필요로한다. 검색 공간이 매우 커질 때 많은 추론이 무너질 것이므로 핵심 근사값을 최적화 할 때 문제가됩니다. 아래의 계수는 현재까지 내 최선의 해결책을 나타내며, 확실히 개선의 여지가 있습니다. 인텔 컴파일러 및 /fp:strict으로 구축하고 2 랜덤 테스트 벡터를 사용하면 관찰 된 최대 오류는 양의 반쪽 평면에서는 2.83788 ulps이고 음의 반 평면에서는 2.77856 ulps입니다.
double my_erfcx (double x)
{
double a, d, e, m, p, q, r, s, t;
a = fmax (x, 0.0 - x); // NaN preserving absolute value computation
/* Compute q = (a-4)/(a+4) accurately. [0,INF) -> [-1,1] */
m = a - 4.0;
p = a + 4.0;
r = 1.0/p;
q = m * r;
t = fma (q + 1.0, -4.0, a);
e = fma (q, -a, t);
q = fma (r, e, q);
/* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1,1] */
p = 0x1.edcad78fc8044p-31; // 8.9820305531190140e-10
p = fma (p, q, 0x1.b1548f14735d1p-30); // 1.5764464777959401e-09
p = fma (p, q, -0x1.a1ad2e6c4a7a8p-27); // -1.2155985739342269e-08
p = fma (p, q, -0x1.1985b48f08574p-26); // -1.6386753783877791e-08
p = fma (p, q, 0x1.c6a8093ac4f83p-24); // 1.0585794011876720e-07
p = fma (p, q, 0x1.31c2b2b44b731p-24); // 7.1190423171700940e-08
p = fma (p, q, -0x1.b87373facb29fp-21); // -8.2040389712752056e-07
p = fma (p, q, 0x1.3fef1358803b7p-22); // 2.9796165315625938e-07
p = fma (p, q, 0x1.7eec072bb0be3p-18); // 5.7059822144459833e-06
p = fma (p, q, -0x1.78a680a741c4ap-17); // -1.1225056665965572e-05
p = fma (p, q, -0x1.9951f39295cf4p-16); // -2.4397380523258482e-05
p = fma (p, q, 0x1.3be1255ce180bp-13); // 1.5062307184282616e-04
p = fma (p, q, -0x1.a1df71176b791p-13); // -1.9925728768782324e-04
p = fma (p, q, -0x1.8d4aaa0099bc8p-11); // -7.5777369791018515e-04
p = fma (p, q, 0x1.49c673066c831p-8); // 5.0319701025945277e-03
p = fma (p, q, -0x1.0962386ea02b7p-6); // -1.6197733983519948e-02
p = fma (p, q, 0x1.3079edf465cc3p-5); // 3.7167515521269866e-02
p = fma (p, q, -0x1.0fb06dfedc4ccp-4); // -6.6330365820039094e-02
p = fma (p, q, 0x1.7fee004e266dfp-4); // 9.3732834999538536e-02
p = fma (p, q, -0x1.9ddb23c3e14d2p-4); // -1.0103906603588378e-01
p = fma (p, q, 0x1.16ecefcfa4865p-4); // 6.8097054254651804e-02
p = fma (p, q, 0x1.f7f5df66fc349p-7); // 1.5379652102610957e-02
p = fma (p, q, -0x1.1df1ad154a27fp-3); // -1.3962111684056208e-01
p = fma (p, q, 0x1.dd2c8b74febf6p-3); // 2.3299511862555250e-01
/* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
d = a + 0.5;
r = 1.0/d;
r = r * 0.5;
q = fma (p, r, r); // q = (p+1)/(1+2*a)
t = q + q;
e = (p - q) + fma (t, -a, 1.0); // residual: (p+1)-q*(1+2*a)
r = fma (e, r, q);
/* Handle argument of infinity */
if (a > 0x1.fffffffffffffp1023) r = 0.0;
/* Handle negative arguments: erfcx(x) = 2*exp(x*x) - erfcx(|x|) */
if (x < 0.0) {
s = x * x;
d = fma (x, x, -s);
e = exp (s);
r = e - r;
r = fma (e, d + d, r);
r = r + e;
if (e > 0x1.fffffffffffffp1023) r = e; // avoid creating NaN
}
return r;
}