Hata işlevinin hatasız şekilde gerçekleştirilmesi hata fonksiyonu erff()

Question

Hata işlevi standart normal dağılım ile yakından ilişkilidir ve doğal bilimlerde olduğu gibi diğer alanlarda da sıkça görülür. Örneğin, fiyatlandırma seçeneklerinde finansmanda kullanılır. Bunun için destek ilk olarak ISO C99'a ve daha sonra C++ 'ya erf(), erff() gibi fonksiyonlar olarak eklenirken, yakın zamana kadar en az bir popüler C/C + takım zincirinden yok olmuştur. Birçok proje hala sıkça yaşlı literatürdeki yaklaşımlar dayalı hata fonksiyonunun kendi uygulamalarını, sırayla geri

Cecil Hastings Jr gitmek gibi Abramowitz and Stegun, , "Dijital Bilgisayarlar için Yaklaşımlar" kullanırlar. Princeton University Press, 1955

Modern hesaplamada, transandantal işlevlerin sadık bir şekilde gerçekleştirilmiş uygulamaları, çoğu zaman matematik kütüphaneleri için asgari bir doğruluk standardı olarak görülür; Böyle bir standart, yüksek performanslı uygulamalara izin verir. Tüm giriş etki alanındaki matematiksel değere kıyasla, 1 ulp'dan daha az bir maksimum hata ile sonuç döndürdüğünde, bir işlev, tam anlamıyla yuvarlatılmış olarak adlandırılır. Eski yayınlanmış algoritmalar, IEEE-754 tekli hassas işlemlerin özel kullanımı ile uygulandığında sadık bir şekilde sonuç vermez.

Modern bilgisayar donanımı, kayan noktalı çarpma işlemi yapan bir kayan nokta işlemi (kayan noktalı çarpma) ve ardından kayan noktalı çarpma işleminin ardından tamamlanmış ürün eklenecek şekilde bir kayan noktalı nokta eklenmesi sağlar. ve operasyonun sonunda sadece tek bir yuvarlama meydana gelir. 1990 yılında IBM tarafından sunulan bu ertelenmiş işlem, birçok hesaplamada daha yüksek doğruluk ve daha yüksek performans sağlar. Bugün en yaygın iki CPU mimarisinde (ARM ve x86) ve GPU'larda mevcuttur. C ve C++ 'da fmaf() ve fmaf() fonksiyonları aracılığıyla ortaya çıkarılmıştır.

FMA'nın donanım tarafından yerel olarak desteklendiğini varsayarsak, her ikisi de hem yuvarlak hem de verimli olan tek duyarlıklı bir hata işlevi erff() nasıl yapılır? Tercihen, kod muhtemelen küçük kod modifikasyonundan sonra vektör edilebilir olmalıdır.

Answer 1

graph of the error function'a bakıldığında, işlevin menşe hakkında simetrik olduğunu gözlemler; Bu nedenle, yaklaşımlar önemsiz olarak pozitif yarı düzlem ile sınırlandırılabilir. Ayrıca, grafik, x = 1 yakınında bir yerde sınırlama olmak üzere iki bölüme ayrılabilir. Kökeni daha yakın olan segmentte, hata fonksiyonu oldukça doğrusal ve dikey olarak domine edilirken, segmentte orijinden daha uzakta yatay olarak baskındır ve üssel olarak çürüyen modadaki birliği asimptotik olarak ele alır.

Makul bir sonuç, x * p (x) 'in basit bir polinom yaklaşımı sıfıra yakın segment için uygunken, diğer segmentin 1 - exp (x * q (x)) ile iyi bir şekilde yakınlaştırılacağıdır. q (x) ikinci bir polinom yaklaşımıdır. fonksiyonunun Taylor serisinin genişlemesine dayanarak, orijine yakın segmentin tahmini, aslında x * p (x) formunda olmalıdır. İlk görev, iki bölüm arasındaki geçiş noktasını bulmaktır. Bunun için, 0.875'teki geçiş noktasıyla başlayan ve onu adım adım 1.0'a doğru ilerleyen deneysel bir yaklaşım kullandım. Geçiş noktasının her bir değeri için, sıfır ile Remez algorithm arasındaki geçiş noktası arasındaki hata fonksiyonuna bir başlangıç ​​minimax yaklaşımı oluşturdum. Bunlar daha sonra, FMA işlemlerinin kullanılmasıyla Horner şeması aracılığıyla polinomun değerlendirilmesine dayanan katsayı değerleri için arama arama yöntemini kullanarak doğruluk bakımından daha da geliştirildi.Geçiş noktasının arttırılması, elde edilen yaklaşık maksimum hata 1 ulp'yi aşana kadar tekrarlandı. Bu süreçle, iki yaklaşım segmenti arasındaki optimal sınırı 0.921875 olarak belirledim. Bu, yaklaşık 1 ulp'den az bir maksimum hataya ulaşan bir x * p (x) yaklaşımıyla sonuçlanır. Remez algoritması ayrıca polinom q (x) için bir başlangıç ​​minimax hesaplaması sağlamak amacıyla da kullanılmıştır. Q (x) ile içsel olarak içersinde kullanılan yaklaştırmanın, hatalarının birbirini güçlendirdiği veya dengelediği şekilde exp() ile arasında çok az bir etkileşim olduğu açıktır. Bu, q (x) katsayıları için en iyi seçimin, exp()'un uygulanmasına sıkı sıkıya bağlı olacağı anlamına gelir ve başlangıç ​​katsayı kümelerinin sezgisel arıtmasının bir parçası olarak dikkate alınmalıdır. Bu nedenle, kendimi herhangi bir kütüphane uygulamasından ayırmak için kendi uygulamamı expf() kullanmaya karar verdim. Asgari gereklilik olarak, expf()'un kendisinin tam anlamıyla yuvarlatılmış olması gerekir ve tam olarak ne kadar sıkı olduğunu belirlemeye çalışmadım, ancak bu yaklaşımın çalışması için biraz daha sıkı bir hataya uyması gerekir. Bu durumda, kendi uyguladığım expf(), yeterli olduğu ortaya çıkarılan 0.87161 ulps'lik bir hata sağlar.

exp kullanımı() gerektiren kademeli yavaş yol olduğu

, bir komut-aşamalı paralellik ve böylece performansını arttırmak için polinom q (X) alt sıra terimleri için Estrin's scheme kullanmayı tercih ettik. Bunu yapmanın doğruluğu üzerindeki etkisi ihmal edilebilir. Doğruluk nedenleriyle, Horner şeması polinomun yüksek mertebeli terimleri için kullanılmalıdır. Her iki polinomun en az anlamlı katsayılarına baktığımızda, her ikisinin de 1.128 olduğunu görüyoruz ... ve bu nedenle, katsayıyı (1 + 0.128 ...) içine bölerek, doğruluğu biraz iyileştirebiliriz. x ile son çarpma. Sonunda, erff() numaralı bir uygulamayı gerçekleştirebildim. Burada iki kod yolunun her biri, daha yüksek hassasiyetli bir referansa karşı kapsamlı bir testle belirlendiği üzere, 1 ulp'nin altında bir maksimum hataya ulaştı. Bu nedenle işlev sadakatle yuvarlanır. FMA kullanımı bu başarının önemli bir bileşenidir. Alet zincirine bağlı olarak, aşağıda gösterilen C99 kodu, olduğu gibi vektör edilebilir olabilir veya her ikisi de, her iki kod yolu, sonuçta istenen sonucun bir seçimiyle eşzamanlı olarak hesaplanacak şekilde manuel olarak modifiye edilebilir. Yüksek performanslı matematik kitaplıkları, my_expf() özel işlevim yerine kullanılması gereken expf()'un bir vektörleştirilebilir sürümünü içerir. expf()'un tüm vektör edilmiş uygulamaları, yeterli bir doğruluk sunmaz ve diğerleri için, polinom q (x) içindeki katsayıların ayarlanması gerekli olacaktır.

expf()'un özel bir sürümü kullanılıyorsa, burada yaptığım gibi, muhtemelen bir ldexpf() numaralı aramayı performans nedenleriyle daha hızlı makine özel kodla değiştirmek isteyebilir.

/* 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; 
} 

/* compute error function. max ulp error = 0.99993 */ 
float my_erff (float a) 
{ 
    float r, s, t, u; 

    t = fabsf (a); 
    s = a * a; 
    if (t > 0.921875f) { // 0.99527 ulp 
     r = fmaf (0x1.222900p-16f, t, -0x1.91d2ccp-12f); // 1.72948930e-5, -3.83208680e-4 
     u = fmaf (0x1.fd1336p-09f, t, -0x1.8d6300p-06f); // 3.88393435e-3, -2.42545605e-2 
     r = fmaf (r, s, u); 
     r = fmaf (r, t, 0x1.b55cb0p-4f); // 1.06777847e-1 
     r = fmaf (r, t, 0x1.450aa0p-1f); // 6.34846687e-1 
     r = fmaf (r, t, 0x1.079d0cp-3f); // 1.28717512e-1 
     r = fmaf (r, t, t); 
     r = my_expf (-r); 
     r = 1.0f - r; 
     r = copysignf (r, a); 
    } else { // 0.99993 ulp 
     r =    -0x1.3a1a82p-11f; // -5.99104969e-4 
     r = fmaf (r, s, 0x1.473f48p-08f); // 4.99339588e-3 
     r = fmaf (r, s, -0x1.b68bd2p-06f); // -2.67667342e-2 
     r = fmaf (r, s, 0x1.ce1a46p-04f); // 1.12818025e-1 
     r = fmaf (r, s, -0x1.8126e0p-02f); // -3.76124859e-1 
     r = fmaf (r, s, 0x1.06eba6p-03f); // 1.28379151e-1 
     r = fmaf (r, a, a); 
    } 
    return r; 
}