Sau khi xem xét các cách tiếp cận khác nhau, là thuật toán được đề xuất trong giấy tờ sau đây:
MM Shepherd và JG Laframboise, "Chebyshev xấp xỉ của (1 + 2 x) exp (x) erfc x trong 0 ≤ x < ∞. " Toán của Tính, Tập 36, số 153, tháng 1 năm 1981, pp. 249-253 (online copy)
Ý tưởng cơ bản của bài báo là tạo ra một xấp xỉ (1 + 2 x) exp (x) erfc (x), từ đó chúng ta có thể tính toán erfcx (x) bằng cách đơn giản chia cho (1 + 2 x), và erfc (x) bởi sau đó nhân với exp (- x). Phạm vi giới hạn chặt chẽ của hàm, với các giá trị hàm gần bằng [1, 1.3], và độ phẳng "chung" của nó cho vay chính nó đến xấp xỉ đa thức.tính số của phương pháp này đang được cải thiện hơn nữa bằng cách thu hẹp khoảng cách xấp xỉ: lập luận ban đầu x được biến đổi bởi q = (x - K)/(x + K), trong đó K là một lựa chọn phù hợp liên tục , tiếp theo là tính toán p (q), trong đó p là đa thức.
Kể từ erfc - x = 2 - erfc x, chúng ta chỉ cần xem xét khoảng [0, ∞] được ánh xạ tới khoảng [-1, 1] bởi sự chuyển đổi này. Đối với độ chính xác đơn IEEE-754, erfcf() biến mất (trở thành 0) cho x> 10.0546875, vì vậy, chỉ cần xem xét x ∈ [0, 10.0546875). Giá trị "tối ưu" của K cho phạm vi này là gì? Tôi biết không có phân tích toán học nào cung cấp câu trả lời, bài báo gợi ý K = 3,75 dựa trên thử nghiệm. Một xấp xỉ đa thức của phương trình minimax bằng 9 là đủ cho các giá trị khác nhau của K trong vùng lân cận chung đó. {2, 2.625, 3.3125} Trong đó, K = 2 là sự lựa chọn thuận lợi nhất, vì nó cho phép tính toán rất chính xác (x - K)/(x + K), như được hiển thị trong this question.
Giá trị K = 2 và tên miền đầu vào cho x sẽ gợi ý rằng cần sử dụng biến thể 4 từ my answer, tuy nhiên một khi có thể chứng minh bằng thực nghiệm rằng biến thể ít tốn kém hơn 5 đạt được độ chính xác tương tự ở đây, có khả năng là do độ dốc rất nông của hàm xấp xỉ cho q> -0,5, gây ra bất kỳ lỗi nào trong đối số q được giảm khoảng 10 lần.
Do tính toán erfc() yêu cầu các bước xử lý ngoài phép tính xấp xỉ ban đầu, rõ ràng là tính chính xác của cả hai tính toán này phải cao để đạt được kết quả cuối cùng chính xác. Sửa lỗi kỹ thuật phải được sử dụng.
Một nhận xét rằng hệ số quan trọng nhất trong xấp xỉ đa thức của (1 + 2 x) exp (x) erfc (x) có dạng (1 + s), nơi s < 0.5. Điều này có nghĩa là chúng tôi có thể đại diện cho hệ số hàng đầu chính xác hơn bằng cách tách 1 và chỉ sử dụng s trong đa thức. Vì vậy, thay vì tính toán đa thức p (q), sau đó nhân với nghịch đảo r = 1/(1 + 2 x), nó tương đương về mặt toán học nhưng có lợi thế về tính toán xấp xỉ lõi như p (q) + 1 và sử dụng p để tính fma (p, r, r).
Độ chính xác của các bộ phận có thể được tăng cường bằng cách tính toán một thương ban đầu q từ đối ứng r, tính dư e = p 1-q * (1 + 2 x) với sự giúp đỡ của một FMA, sau đó sử dụng e để áp dụng các chỉnh q = q + (e * r), một lần nữa bằng cách sử dụng một FMA.
Lũy thừa có tính phóng đại lỗi, do đó tính toán của e-x phải được thực hiện một cách cẩn thận. Tính khả dụng của FMA trivially cho phép việc tính toán - x như một đôi float s cao: s thấp. e x là đạo hàm riêng của mình, vì vậy người ta có thể tính toán e s cao: s thấp như e s cao + e s cao * s thấp. tính toán này có thể được kết hợp với việc nhân rộng các kết quả trung gian trước r để mang r = r * e s cao + r * e s cao * s thấp . Bằng cách sử dụng FMA, một đảm bảo rằng cụm từ quan trọng nhất r * e s cao được tính chính xác nhất có thể.
Kết hợp các bước trên với một vài lựa chọn đơn giản để xử lý các trường hợp ngoại lệ và lập luận tiêu cực, một đến vào mã C sau:
float my_expf (float);
/*
* 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_erfcf (float x)
{
float a, d, e, m, p, q, r, s, t;
a = fabsf (x);
/* Compute q = (a-2)/(a+2) accurately. [0, 10.0546875] -> [-1, 0.66818] */
m = a - 2.0f;
p = a + 2.0f;
r = 1.0f/p;
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, 0.66818] */
p = -0x1.a48024p-12f; // -4.01020574e-4
p = fmaf (p, q, -0x1.42a172p-10f); // -1.23073824e-3
p = fmaf (p, q, 0x1.585784p-10f); // 1.31355994e-3
p = fmaf (p, q, 0x1.1ade24p-07f); // 8.63243826e-3
p = fmaf (p, q, -0x1.081b72p-07f); // -8.05991236e-3
p = fmaf (p, q, -0x1.bc0b94p-05f); // -5.42047396e-2
p = fmaf (p, q, 0x1.4ffc40p-03f); // 1.64055347e-1
p = fmaf (p, q, -0x1.540840p-03f); // -1.66031361e-1
p = fmaf (p, q, -0x1.7bf612p-04f); // -9.27639678e-2
p = fmaf (p, q, 0x1.1ba03ap-02f); // 2.76978403e-1
/* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
t = a + a;
d = t + 1.0f;
r = 1.0f/d;
q = fmaf (p, r, r); // q = (p+1)/(1+2*a)
e = (p - q) + fmaf (q, -t, 1.0f); // (p+1) - q*(1+2*a)
r = fmaf (e, r, q);
/* Multiply by exp(-a*a) ==> erfc(a) */
s = a * a;
e = my_expf (-s);
t = fmaf (a, -a, s);
r = fmaf (r, e, r * e * t);
/* Handle NaN arguments to erfc() */
if (!(a <= 0x1.fffffep127f)) r = x + x;
/* Clamp result for large arguments */
if (a > 10.0546875f) r = 0.0f;
/* Handle negative arguments to erfc() */
if (x < 0.0f) r = 2.0f - r;
return r;
}
/* Compute exponential base e. Maximum ulp error = 0.87161 */
float my_expf (float a)
{
float c, f, r;
int i;
// exp(a) = exp(i + f); i = rint (a/log(2))
c = 0x1.800000p+23f; // 1.25829120e+7
r = fmaf (0x1.715476p+0f, a, c) - c; // 1.44269502e+0
f = fmaf (r, -0x1.62e400p-01f, a); // -6.93145752e-1 // log_2_hi
f = fmaf (r, -0x1.7f7d1cp-20f, f); // -1.42860677e-6 // log_2_lo
i = (int)r;
// approximate r = exp(f) on interval [-log(2)/2,+log(2)/2]
r = 0x1.6a98dap-10f; // 1.38319808e-3
r = fmaf (r, f, 0x1.1272cap-07f); // 8.37550033e-3
r = fmaf (r, f, 0x1.555a20p-05f); // 4.16689515e-2
r = fmaf (r, f, 0x1.55542ep-03f); // 1.66664466e-1
r = fmaf (r, f, 0x1.fffff6p-02f); // 4.99999851e-1
r = fmaf (r, f, 0x1.000000p+00f); // 1.00000000e+0
r = fmaf (r, f, 0x1.000000p+00f); // 1.00000000e+0
// exp(a) = 2**i * exp(f);
r = ldexpf (r, i);
// handle special cases
if (!(fabsf (a) < 104.0f)) {
r = a + a; // handle NaNs
if (a < 0.0f) r = 0.0f;
if (a > 0.0f) r = 1e38f * 1e38f; // + INF
}
return r;
}
Tôi đã từng thực hiện riêng của tôi về expf() trong đoạn code trên để cô lập của tôi hoạt động từ những khác biệt trong triển khai expf() trên các nền tảng tính toán khác nhau. Nhưng bất kỳ việc thực hiện nào của expf() có lỗi tối đa là gần 0,5 ulp sẽ hoạt động tốt. Như được hiển thị ở trên, tức là khi sử dụng my_expf(), my_erfcf() có lỗi tối đa là 2,65712 ulps. Cung cấp tính sẵn có của một vectorizable expf(), mã ở trên nên vectorize mà không có vấn đề. Tôi đã kiểm tra nhanh với trình biên dịch Intel 13.1.3.198.Tôi đặt một cuộc gọi đến my_erfcf() trong một vòng lặp, thêm #include <mathimf.h>, thay thế các cuộc gọi đến my_expf() với một cuộc gọi đến expf(), sau đó biên soạn sử dụng các dòng lệnh chuyển mạch:
/Qstd=c99 /O3 /QxCORE-AVX2 /fp:precise /Qfma /Qimf-precision:high:expf /Qvec_report=2
Trình biên dịch Intel thông báo rằng vòng lặp đã được vector hóa, mà tôi đã kiểm tra lại bằng cách kiểm tra mã nhị phân đã tháo rời.
Vì my_erfcf() chỉ sử dụng các đối ứng chứ không phải là các bộ phận đầy đủ, nó phù hợp với việc sử dụng các triển khai nhanh đối ứng, miễn là chúng cung cấp kết quả gần như chính xác. Đối với các bộ vi xử lý cung cấp một xấp xỉ tương đương chính xác đơn nhanh trong phần cứng, điều này có thể dễ dàng đạt được bằng cách ghép nối nó với một phép lặp Halley với sự hội tụ khối. Ví dụ A (vô hướng) của phương pháp này cho bộ vi xử lý x86 là:
/* Compute 1.0f/a almost correctly rounded. Halley iteration with cubic convergence */
float fast_recipf (float a)
{
__m128 t;
float e, r;
t = _mm_set_ss (a);
t = _mm_rcp_ss (t);
_mm_store_ss (&r, t);
e = fmaf (r, -a, 1.0f);
e = fmaf (e, e, e);
r = fmaf (e, r, r);
return r;
}