Saya telah menulis sebuah kode waktu lalu yang berusaha untuk menghitung tanpa menggunakan fungsi perpustakaan. Kemarin, saya meninjau kode lama, dan saya mencoba membuatnya secepat mungkin, (dan benar). Inilah upaya saya sejauh ini:
const double ee = exp(1);
double series_ln_taylor(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 )
n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(1 - x) = -x - x**2/2 - x**3/3... */
n = 1 - n;
now = term = n;
for ( i = 1 ; ; ){
lgVal -= now;
term *= n;
now = term / ++i;
if ( now < 1e-17 ) break;
}
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Di sini saya mencoba untuk menemukan sehingga e a hanya lebih n, dan kemudian saya menambahkan nilai logaritma dari n , yang kurang dari 1. Pada titik ini, ekspansi Taylor darilog(1-x)dapat digunakan tanpa mengkhawatirkan.
Saya baru-baru ini semakin tertarik pada analisis numerik, dan itulah sebabnya saya tidak bisa tidak bertanya, seberapa cepat segmen kode ini dapat dijalankan dalam praktik, sementara cukup benar? Apakah saya perlu beralih ke beberapa metode lain, misalnya, menggunakan fraksi lanjutan, seperti ini ?
The fungsi yang disediakan dengan C standar perpustakaan hampir 5,1 kali lebih cepat dari implementasi ini.
UPDATE 1 : Menggunakan seri arctan hiperbolik yang disebutkan di Wikipedia , perhitungannya tampaknya hampir 2,2 kali lebih lambat daripada fungsi log library C standar. Padahal, saya belum memeriksa kinerja secara ekstensif, dan untuk jumlah yang lebih besar, implementasi saya saat ini tampaknya SANGAT lambat. Saya ingin memeriksa kedua implementasi saya untuk batas kesalahan dan waktu rata-rata untuk berbagai nomor jika saya dapat mengaturnya. Ini adalah usaha kedua saya.
double series_ln_arctanh(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
for ( i = 3 ; ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Setiap saran atau kritik dihargai.
UPDATE 2: Berdasarkan saran yang dibuat di bawah ini, saya telah menambahkan beberapa perubahan tambahan di sini, yang sekitar 2,5 kali lebih lambat dari implementasi perpustakaan standar. Namun, saya telah mengujinya hanya untuk bilangan bulat kali ini, untuk jumlah yang lebih besar runtime akan meningkat. Untuk sekarang. Saya belum tahu teknik untuk menghasilkan angka ganda acak ≤ 1 e 308 , jadi belum sepenuhnya menjadi patokan. Untuk membuat kode lebih kuat, saya telah menambahkan koreksi untuk kasus sudut. Kesalahan rata-rata untuk tes yang saya buat adalah sekitar 4 e - 15 .
double series_ln_better(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n == 0 ) return -1./0.; /* -inf */
if ( n < 0 ) return 0./0.; /* NaN*/
if ( n < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
/* the cutoff iteration is 650, as over e**650, term multiplication would
overflow. For larger numbers, the loop dominates the arctanh approximation
loop (with having 13-15 iterations on average for tested numbers so far */
for ( term = 1; term < n && lgVal < 650 ; term *= ee, lgVal++ );
if ( lgVal == 650 ){
n /= term;
for ( term = 1 ; term < n ; term *= ee, lgVal++ );
}
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
/* limiting the iteration for worst case scenario, maximum 24 iteration */
for ( i = 3 ; i < 50 ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
sumber
frexpuntuk mendapatkan mantissa dan eksponen nomor tersebutJawaban Kirill sudah menyentuh sejumlah besar masalah yang relevan. Saya ingin memperluas beberapa dari mereka berdasarkan pengalaman desain perpustakaan matematika praktis. Catatan di muka: perancang perpustakaan matematika cenderung menggunakan setiap optimasi algoritme yang dipublikasikan, serta banyak optimisasi spesifik mesin, tidak semuanya akan dipublikasikan. Kode sering akan ditulis dalam bahasa assembly, daripada menggunakan kode yang dikompilasi. Oleh karena itu tidak mungkin bahwa implementasi langsung dan dikompilasi akan mencapai lebih dari 75% dari kinerja implementasi perpustakaan matematika berkualitas tinggi yang ada, dengan asumsi set fitur yang identik (akurasi, penanganan kasus khusus, pelaporan kesalahan, dukungan mode pembulatan).
Akurasi biasanya dinilai dengan perbandingan dengan referensi presisi tinggi (pihak ketiga). Argumen tunggal Fungsi presisi tunggal dapat dengan mudah diuji secara mendalam, fungsi lainnya memerlukan pengujian dengan vektor uji acak (diarahkan). Jelas seseorang tidak dapat menghitung hasil referensi yang sangat akurat, tetapi penelitian terhadap Dilema Table-Maker's menunjukkan bahwa untuk banyak fungsi sederhana, cukup untuk menghitung referensi dengan ketepatan sekitar tiga kali ketepatan target. Lihat misalnya:
Vincent Lefèvre, Jean-Michel Muller, "Kasus Terburuk untuk Pembulatan Benar dari Fungsi Dasar dalam Presisi Ganda". Dalam Prosiding Simposium IEEE 15 tentang Aritmatika Komputer , 2001,111-118). (pracetak online)
Dalam hal kinerja, kita harus membedakan antara mengoptimalkan latensi (penting ketika kita melihat waktu pelaksanaan operasi tergantung), versus mengoptimalkan untuk throughput (relevan ketika mempertimbangkan waktu pelaksanaan operasi independen). Selama dua puluh tahun terakhir, proliferasi teknik paralelisasi perangkat keras seperti paralelisme tingkat instruksi (misalnya superscalar, prosesor out-of-order), paralelisme tingkat data (misalnya instruksi SIMD), dan paralelisme tingkat-benang (misalnya hyper-threading, prosesor multi-core) telah mengarah pada penekanan pada throughput komputasi sebagai metrik yang lebih relevan.
Operasi penggandaan multiply-add ( FMA ), pertama kali diperkenalkan oleh IBM 25 tahun yang lalu, dan sekarang tersedia di semua arsitektur prosesor utama, adalah blok bangunan penting dari implementasi perpustakaan matematika modern. Ini memberikan pengurangan kesalahan pembulatan, memberikan perlindungan terbatas terhadap pembatalan subtraktif , dan sangat menyederhanakan aritmatika ganda-ganda .
C99log()C99fma()sumber