Penghitungan yang akurat dan cepat untuk bilangan floating point menggunakan contoh fungsi sinus. Bagian 2: libm

Saya melanjutkan seri artikel tentang bekerja dengan floating point. Dalam artikel pertama, saya memberikan pengantar matematika singkat dan menunjukkan cara paling sederhana dan paling jelas untuk menghitung sinus dengan contoh program dengan berbagai "perangkap". Artikel hari ini akan sedikit berbeda gayanya. Tidak akan ada latihan di sini, tapi kita akan menggali lebih dalam matematika dan masuk ke maha kudus - kode dari perpustakaan standar. Saya juga akan memberikan jawaban atas pertanyaan di akhir artikel pertama. Jadi ayo pergi.

Beberapa matematika lagi

Jelas, semakin cepat deret Taylor menurun dalam nilai absolut, semakin sedikit istilah yang dibutuhkan untuk mencapai akurasi yang dibutuhkan. Sehingga tampaknya hasilnya akan lebih akurat (di bawah ini akan dibahas lebih detail). Sebagai perbandingan, misalnya, ambil suku dari deret Taylor derajat ketujuh ($$$\frac{x^7}{7!}$) di $$$x=1.0$ dan $$$x=0.1$... Nilai ekspresi akan menjadi$$$1.198\cdot 10^{-4}$ dan $$$1.198\cdot 10^{-11}$masing-masing. Perbedaan yang sangat besar, bukan? Jadi, mari kita coba mencari cara untuk mengurangi batas atas interval kalkulasi fungsi sinus.

Perluasan seri di sekitar nilai yang diberikan

Untuk memahami metode ini, kita perlu kembali ke tahun pertama institut dan mengingat definisi deret Taylor ( wiki ). Singkatnya: mengetahui fungsi dan turunannya di beberapa titik, Anda dapat menemukan nilai fungsi di sekitar titik ini dengan mengembang menjadi deret Taylor. Untuk fungsi sinus, artinya sebagai berikut

$$

$$\sin(x_0+\Delta x) \approx \sin(x_0)+\sin'(x_0)\frac{\Delta x}{1!}+\sin''(x_0)\frac{\Delta x^2}{2!}+\sin'''(x_0)\frac{\Delta x^3}{3!}+\cdots$$

Apa yang diberikan pendekatan ini dari sudut pandang praktis? Bayangkan kita memiliki jeda dari$$$0$ sebelum $$$\pi/2$... Mari kita pilih 10 titik terdistribusi linier pada interval ini (pilihannya tidak optimal):$$$x_0=0$, $$$x_1=\pi/20$, $$$x_2=2\pi/20$, $$$x_i=i\cdot \pi/20$... Untuk setiap titik, hitung pelat dengan sinus dan turunannya pada titik ini. Sekarang Anda dapat memodifikasi fungsinya sehingga ketika Anda mendapatkan nilainya$$$x$ fungsi mengambil nilai terdekat $$$x_i$ dan menjabarkan secara berurutan di sekitar titik $$$x_i$, bukan sekitar nol ($$$\Delta x = x - x_i$).

Menggunakan transformasi trigonometri

Jika kita kembali lebih jauh, ke kelas senior sekolah, maka kita dapat mengingat satu formula yang sangat penting:

$$

$$\sin(x_0+\Delta x)=\sin(x_0)\cos(\Delta x)+\cos(x_0)\sin(\Delta x)$$

Dan kemudian semuanya sama seperti di paragraf sebelumnya. Kami memilih titik dalam interval, menghitung sinus dan kosinus untuk mereka, dan saat memanggil fungsi sinus, kami mencari yang terdekat dan, dengan menggunakan rumus di atas, menghitung sinus menggunakan nilai kecil$$$\Delta x$...

Pikirkan tentang mana dari kedua metode ini yang lebih baik untuk dipilih, tetapi untuk saat ini kita akan beralih dari matematika ke perhitungan praktis.

Properti distribusi perkalian di dunia floating point

Saya harus meminta nasihat dari Internet, apa namanya $$$a\cdot(b+c)=a\cdot b + a\cdot c$... Ternyata itu adalah properti distributif. Mari kita kembali ke pertanyaan yang saya tanyakan di akhir bagian pertama. Yakni, mengapa ekspresi matematis setara$$$a\cdot(b+c)$ dan $$$a\cdot b + a\cdot c$dapat memberikan hasil yang berbeda dalam perhitungan floating point? Cara termudah untuk mengilustrasikannya adalah dengan sebuah contoh. Mari kita ambil sistem hipotetis yang bekerja dengan bilangan floating point dalam format desimal dengan presisi 4 digit. Mari kita anggap itu$$$b=1.000\text{E}0$, $$$c=1.234\text{E-}2$, dan $$$a=5.678\text{E-}2$... Pertama, mari kita ambil ekspresi dengan tanda kurung dan menghitungnya selangkah demi selangkah, ingat untuk membulatkannya di setiap langkah:

1)$$$b+c=1.000\text{E}0 + 1.234\text{E-}2=1.012\text{E}0$

2) $$$a\cdot(b+c)=5.678\text{E-}2\times 1.012\text{E}0=5.746\text{E-}2$

Tanggapan yang diterima $$$y=5.746\text{E-}2$

Sekarang mari kita hitung ekspresi kedua dengan cara yang sama selangkah demi selangkah:

$$$a\cdot b + a\cdot c=5.678\text{E-}2\times 1.000\text{E}0 + 5.678\text{E-}2\times 1.234\text{E-}2=5.678\text{E-}2 + 7.007\text{E-}4=5.748\text{E-}2$

Tanggapan yang diterima $$$y=5.748\text{E-}2$

Jawaban yang benar adalah 0,0574806652.



Seperti yang Anda lihat, jawaban yang diperoleh dalam kasus kedua jauh lebih dekat dengan yang benar daripada di kasus pertama. Jika kita menjelaskan ini dengan jari, maka bayangkan ketika dalam kasus pertama kita menambahkan angkanya menjadi 1,0$$$c=1.234\text{E-}2$kami hanya membuang dua digit terakhir. Mereka tidak lebih. Dalam kasus kedua, pembuangan terjadi di bagian paling akhir, setelah perkalian. Itu. dalam kasus kedua, operasi perkalian lebih akurat.



Tampaknya Anda bisa menyelesaikan ini, tetapi lihat lebih dekat pada metode pertama dan beri tahu saya apa hasil perhitungannya.$$$b+c-b$... Dan ... kami memiliki cara untuk membulatkan angka floating point! Jangan lewatkan contoh ini. Beri diri Anda waktu untuk memikirkannya. Pembulatan angka akan sangat intensif digunakan oleh kami lebih lanjut dalam artikel ini dan selanjutnya.



Mari kita perhatikan satu fitur lagi dari ekspresi ini. Bayangkan bahwa presisi 4 digit variabel tidak cukup bagi kita. Apa yang harus dilakukan? Dan di sini kita sudah memiliki jawabannya - untuk mewakili angka dalam formulir$$$b+c$dan menyimpannya dalam memori sebagai jumlah dari dua digit. Dan, karenanya, lakukan operasi (misalnya perkalian) secara terpisah untuk kedua suku. Teknik ini dijelaskan lebih detail di artikel Menambahkan dua bilangan floating point tanpa kehilangan presisi .



Pada artikel sebelumnya, saya juga menulis metode tersebut$$$a\cdot(b+c)$ada satu fitur yang tidak menyenangkan. Dan itu adalah sebagai berikut. Jumlah$$$c$ selalu terpotong pada digit signifikan terakhir sebuah angka $$$b$... Artinya berapapun jumlahnya$$$c$, jika sebuah $$$b+c \neq b$, maka kesalahan pada tanda terakhir selalu mungkin terjadi bahkan untuk yang kecil $$$c$... Ini tidak diperbolehkan dalam pendekatan di bab berikutnya.

Cara kerjanya menggunakan pustaka GNU sebagai contoh

Bagaimana itu? Sudahkah Anda memilih metode mana dari dua metode yang dijelaskan di awal artikel yang Anda pilih untuk penghitungan sinus yang akurat? Apapun metode yang Anda pilih, keduanya benar. Apalagi mereka benar-benar identik. Percayalah, lihatlah. Di bawah ini saya akan menggunakan rumus sekolah. Mereka lebih mudah dijelaskan.



Berbekal pengetahuan yang diperoleh di artikel sebelumnya dan di artikel ini, Anda dapat dengan mudah memahami kode pustaka standar. Mari buka file s_sin.c dan temukan fungsi __sin di sana :



kodenya cukup sederhana. Sangat mudah untuk memahami bahwa ia memanggil sekumpulan fungsi yang berbeda tergantung pada batasan variabel input. Pada artikel ini, kita akan membahas bagian kode 218-224 untuk sudut 2 ^ -26 <| x | <0.855469. Anda dapat melihat bahwa di bagian kode ini fungsi do_sin (x, 0) dipanggil. Kami akan membahas fungsi ini lebih detail:

1. , dx=0 .
2. 129-130 , abs(x)<0.126, .. x , . , , , .
3. 136-137. , . x 2 . u x. , 0.345678. u=0.34, 0.005678.
4. 140-142. ( s ) ( c ) x . , cos(x)=1-c, 1.0, (. ), .
5. 143. u. , u=0.34 34. sin(u)=sn+ssn, cos(u)=cs+ccs. sn cs — «» u, ssn ccs — .
6. 144-145. sin(u+x)=(sn+ssn)*(1-c)+(cs+ccs)*s. , , 144-145. — .

Faktanya, saya hanya menjelaskan bagian paling sederhana dari menghitung sinus dengan cara ini. Ada banyak matematika yang tertinggal. Misalnya, bagaimana Anda menghitung ukuran tabel dan elemen di dalamnya? Dari manakah angka ajaib 0,126 dan 0,855469 itu berasal? Kapan memotong penghitungan dengan bilangan Taylor? Koreksi koefisien deret Taylor untuk memperhalus hasil.



Semua ini, tentu saja, menarik, tetapi, secara obyektif, metode yang disajikan memiliki banyak kelemahan: perlu untuk menghitung sinus (s) dan cosinus (c) secara bersamaan, yang membutuhkan dua kali lebih banyak perhitungan dari deret Taylor ¹ . Perkalian dengan nilai tabel, seperti yang bisa kita lihat, juga tidak gratis. Selain itu, menyimpan tabel 3520 byte dalam RAM, tentu saja, tidak menjadi masalah, tetapi mengaksesnya (bahkan di cache) bisa mahal.



Oleh karena itu, pada bagian selanjutnya kita akan mencoba menyingkirkan pelat dan menghitung sinus dalam interval [0,126, 0,855469] secara langsung, tetapi lebih akurat daripada pada bab pertama.



Sebelum mengakhiri - pertanyaan tentang kecerdasan. Angka besar dalam contoh ini adalah 52776558133248 = 3 * 2 ⁴⁴ . Dari mana asal bilangan itu, bukan, misalnya, 2 ⁴⁵ ? Saya akan merumuskan pertanyaan lebih tepat. Mengapa angka 3 * 2 ^N optimal saat membulatkan angka , dan bukan, misalnya, 2 ^{N + 1} ? Pertanyaan lain, N mana yang harus Anda pilih untuk membulatkan angka menjadi integer?



¹ Perlu dicatat bahwa keuntungan signifikan dari pendekatan ini dapat muncul ketika sinus dan cosinus dihitung secara bersamaan dari sudut yang sama. Fungsi kedua dapat dihitung hampir gratis.