Saya terkejut menemukan bahwa sulit untuk menemukan informasi tentang masalah ini dalam bahasa Rusia, seolah-olah hanya sedikit orang yang peduli bahwa perpustakaan matematika yang digunakan dalam penyusun modern terkadang tidak memberikan hasil yang bulat dengan benar. Saya khawatir tentang situasi ini, karena saya hanya mengerjakan pengembangan perpustakaan matematika semacam itu. Dalam literatur asing, masalah ini tercakup dengan baik, jadi saya memutuskan untuk menyajikannya dalam bahasa Rusia dalam bentuk ilmiah yang populer, dengan mengandalkan sumber-sumber Barat dan masih sedikit pengalaman pribadi.

---

Teman-teman, demi kenyamanan Anda, artikel ini juga tersedia sebagai presentasi dalam bentuk video (hampir 34 menit), format ini lebih cocok untuk para pembaca yang merasa kesulitan untuk membangun gambar matematika yang diperlukan di kepala mereka, karena banyak materi ilustrasi dalam presentasi. Informasi dalam video tersebut sepenuhnya identik dengan konten artikel. Harap bertindak sesuai keinginan Anda.

---

Saya ulangi bahwa ini bukan artikel ilmiah, tetapi artikel sains populer, setelah membacanya, Anda akan secara singkat mengenalnya.

• Fungsi dasar transendental (exp, sin, log, cosh dan lain-lain) yang bekerja dengan aritmatika floating point dibulatkan secara tidak benar, terkadang membuat kesalahan pada bit terakhir.
• Alasan kesalahan tidak selalu terletak pada kemalasan atau kualifikasi pengembang yang rendah, tetapi dalam satu keadaan mendasar, yang belum dapat diatasi oleh sains modern.
• «», - .
• , , , , exp2(x) pow(2.0, x).

Untuk memahami artikel ini, Anda harus terbiasa dengan format titik mengambang IEEE-754. Sudah cukup jika Anda setidaknya hanya memahami bahwa, misalnya, ini adalah: 0x400921FB54442D18 - angka pi dalam format presisi ganda (binary64, atau ganda), yaitu, Anda hanya memahami apa yang saya maksud dengan rekaman ini; Saya tidak menuntut untuk dapat melakukan transformasi seperti itu dengan cepat. Dan saya akan mengingatkan Anda tentang mode pembulatan dalam artikel ini, ini adalah bagian penting dari cerita. Juga diinginkan untuk mengetahui bahasa Inggris "programmer", karena akan ada istilah dan kutipan dari literatur Barat, tetapi Anda dapat bertahan dengan penerjemah online.



Contoh dulu, supaya kamu langsung paham apa pokok pembicaraannya. Sekarang saya akan memberikan kodenya dalam C ++, tetapi jika ini bukan bahasa Anda, maka saya yakin Anda masih akan dengan mudah memahami apa yang tertulis. Silakan lihat kode ini:

#include <stdio.h>
#include <cmath>

int main() {
  float x = 0.00296957581304013729095458984375f;  // ,  .
  float z;
  z = exp2f(x);  // z = 2**x  .
  printf ("%.8f\n", z);  //      8   .
  z = powf(2.0f, x);  // z = 2**x  
  printf ("%.8f\n", z);  //   .
  return 0;
}

Bilangan x sengaja ditulis dengan sejumlah digit signifikan sehingga tepat dapat direpresentasikan dalam tipe float, yaitu, sehingga compiler akan mengubahnya menjadi kode biner tanpa pembulatan. Bagaimanapun, Anda tahu betul bahwa beberapa kompiler tidak dapat membulatkan tanpa kesalahan (jika Anda tidak tahu, sebutkan di komentar, saya akan menulis artikel terpisah dengan contoh). Selanjutnya dalam program ini, kita perlu menghitung 2 ^x , tetapi mari kita lakukan dengan dua cara: fungsi exp2f (x), dan eksponen eksplisit dari dua powf (2.0f, x). Hasilnya, tentu saja, akan berbeda, karena saya katakan di atas bahwa fungsi dasar tidak dapat berfungsi dengan benar di semua kasus, dan saya secara khusus memilih contoh untuk menunjukkan ini. Inilah hasilnya:

1.00206053
1.00206041

Empat kompiler memberi saya nilai-nilai ini: Microsoft C ++ (19.00.23026), Intel C ++ 15.0, GCC (6.3.0) dan Clang (3.7.0). Mereka berbeda dalam satu hal yang paling tidak signifikan. Berikut adalah kode heksadesimal untuk angka-angka ini:

0x3F804385  // 
0x3F804384  // 

Harap diingat contoh ini, di atasnya kita akan melihat esensi masalah sedikit nanti, tetapi untuk saat ini, agar Anda mendapatkan kesan yang lebih jelas, silakan lihat contoh untuk tipe data presisi ganda (ganda, binary64) dengan beberapa fungsi dasar lainnya. Saya menyajikan hasilnya di tabel. Jawaban yang benar (jika tersedia) memiliki * di bagian akhir.







Saya harap Anda tidak mendapat kesan bahwa saya sengaja mengambil beberapa tes yang benar-benar unik yang sulit Anda temukan? Jika demikian, mari masak dengan lutut kita penghitungan lengkap semua argumen pecahan yang mungkin untuk fungsi 2 ^x untuk tipe data float. Jelas bahwa kita hanya tertarik pada nilai x antara 0 dan 1, karena argumen lain akan menghasilkan hasil yang hanya berbeda pada nilai dalam bidang eksponen dan tidak menarik. Anda sendiri mengerti:

$$

$$2^x = 2^{[x]}\cdot2^{\{x\}}.$$

Setelah menulis program seperti itu (teks tersembunyi akan berada di bawah), saya memeriksa fungsi exp2f dan berapa banyak nilai yang salah yang dihasilkannya pada interval x dari 0 hingga 1.







Jelas dari program di bawah ini bahwa jumlah argumen x yang diuji adalah 197612997. Ternyata, misalnya, Microsoft C ++ salah menghitung fungsi 2 ^x untuk hampir satu persen darinya. Jangan bersuka cita, para penggemar GCC dan Clang, hanya saja fungsi ini diimplementasikan dengan benar di compiler ini, tetapi penuh dengan error pada compiler lainnya.





Saya tidak akan membuat pembaca bosan dengan contoh-contoh ini, hal utama di sini adalah untuk menunjukkan bahwa implementasi modern dari fungsi transendental dapat melengkapi bit terakhir dengan tidak benar, dan kompiler yang berbeda membuat kesalahan di tempat yang berbeda, tetapi tidak satupun dari mereka akan bekerja dengan benar. Omong-omong, Standar IEEE-754 memungkinkan kesalahan ini di bit terakhir (yang akan saya bicarakan di bawah), tetapi masih tampak aneh bagi saya: ok ganda, ini adalah tipe data yang besar, tetapi float dapat diperiksa dengan kekerasan! Apakah itu sulit dilakukan? Tidak sulit sama sekali, dan saya sudah menunjukkan contoh.



Kode enumerasi kami berisi fungsi "menulis sendiri" dengan perhitungan yang benar 2 ^xmenggunakan pendekatan polinomial derajat 10, dan itu ditulis dalam beberapa menit, karena polinomial tersebut diturunkan secara otomatis, misalnya, dalam sistem aljabar komputer Maple. Ini cukup untuk menetapkan kondisi polinomial untuk memberikan presisi 54 bit (untuk fungsi ini, 2 ^x ). Mengapa 54? Tetapi Anda akan segera mengetahuinya, tepat setelah saya memberi tahu Anda esensi masalahnya dan memberi tahu Anda mengapa, pada prinsipnya, sekarang tidak mungkin untuk membuat fungsi transendental yang cepat dan benar untuk tipe data presisi empat kali lipat (binary128), meskipun sudah ada upaya untuk menyerang masalah ini dalam teori.

Pembulatan default dan masalah dengannya

Jika Anda tidak tenggelam dalam pengembangan perpustakaan matematika, maka tidak ada salahnya Anda melupakan aturan pembulatan default untuk bilangan floating point menurut Standar IEEE-754. Oleh karena itu, saya akan mengingatkan Anda tentang itu. Jika Anda mengingat semuanya dengan baik, lihat setidaknya akhir dari bagian ini, Anda akan terkejut: Saya akan menunjukkan kepada Anda situasi di mana membulatkan angka bisa sangat sulit.



Anda dapat dengan mudah mengingat apa itu "round up" (to plus infinity), "round down" (to minus infinity) atau "round to zero" dengan namanya (jika ada, ada Wikipedia). Kesulitan utama bagi programmer muncul dengan pembulatan "ke yang terdekat, tetapi dalam kasus jarak yang sama dari yang terdekat - ke yang dengan digit terakhir genap". Ya, begitulah cara mode pembulatan ini diterjemahkan, yang oleh literatur Barat disebut singkatnya: "Bulatkan ikatan terdekat dengan genap".



Mode pembulatan ini digunakan secara default dan berfungsi sebagai berikut. Jika, sebagai hasil perhitungan, panjang mantissa ternyata lebih besar dari yang dapat ditampung oleh tipe data yang dihasilkan, pembulatan dilakukan ke nilai terdekat dari dua kemungkinan nilai. Namun, situasi mungkin muncul ketika bilangan asli ternyata tepat di tengah antara dua yang terdekat, maka hasil yang dipilih adalah bit terakhir (setelah pembulatan) ternyata genap, yaitu sama dengan nol. Pertimbangkan empat contoh di mana Anda perlu membulatkan menjadi dua bit setelah koma desimal biner:

1. Putaran 1,00 1 001. Bit ketiga setelah koma desimal adalah 1, tetapi kemudian ada bit keenam lagi, yaitu 1, yang berarti pembulatannya akan naik, karena bilangan aslinya lebih dekat dengan 1,01 daripada 1,00.
2. 1,001000. , 1,00 1,01, .
3. 1,011000. 1,01 1,10. , .
4. 1,010111. , 1,01, 1,10.

Dari contoh-contoh ini, mungkin terlihat bahwa semuanya sederhana, tetapi sebenarnya tidak. Faktanya adalah terkadang kita tidak bisa mengatakan dengan pasti apakah kita benar-benar berada di tengah-tengah antara dua nilai. Lihat contohnya. Misalkan kita ingin membulatkan lagi menjadi dua bit setelah titik desimal:



1,00 1 000000000000000000000000000000000000001



Sekarang jelas bagi Anda bahwa pembulatan harus naik, yaitu ke angka 1,01. Namun, Anda melihat angka dengan 40 bit di belakang koma. Bagaimana jika algoritme Anda tidak dapat memberikan presisi 40 bit dan hanya mencapai 30 bit? Maka itu akan memberikan nomor lain:



1.00 1 000000000000000000000000000



Tidak menyadari bahwa pada posisi ke-40 (yang tidak dapat dihitung oleh algoritme) akan ada yang disukai, Anda membulatkan angka ini ke bawah dan mendapatkan 1,00, yang mana itu salah. Anda salah membulatkan bagian terakhir - itulah subjek diskusi kita. Dari penjelasan di atas, ternyata untuk mendapatkan bit ke-2 yang benar saja, Anda harus menghitung fungsinya hingga 40 bit! Wow! Dan jika "lokomotif" nol ternyata lebih panjang? Itulah yang akan kita bicarakan di bagian selanjutnya.



Ngomong-ngomong, ini adalah kesalahan yang dilakukan banyak penyusun saat mengonversi notasi desimal dari bilangan titik mengambang ke dalam format biner yang dihasilkan. Jika angka desimal asli dalam kode program terlalu dekat dengan tengah di antara dua nilai biner yang dapat direpresentasikan secara akurat, maka angka tersebut tidak akan dibulatkan dengan benar. Tapi ini bukan topik artikel ini, tapi alasan untuk cerita tersendiri.

Inti dari masalah pembulatan bit signifikan terakhir

Masalahnya memanifestasikan dirinya karena dua alasan. Yang pertama adalah penolakan yang disengaja terhadap kalkulasi yang memakan waktu demi kecepatan. Dalam hal ini, selama akurasi yang ditentukan diamati, dan bit apa yang akan ada dalam respons adalah masalah sekunder. Alasan kedua adalah Dilema Pembuat Tabel, yang merupakan topik utama percakapan kami. Mari pertimbangkan kedua alasan tersebut secara lebih rinci.

Alasan pertama

Anda, tentu saja, memahami bahwa kalkulasi fungsi transendental diimplementasikan dengan beberapa metode perkiraan, misalnya, dengan metode pendekatan polinomial atau bahkan (jarang) dengan ekspansi deret. Untuk membuat penghitungan terjadi secepat mungkin, pengembang setuju untuk melakukan iterasi sesedikit mungkin dari metode numerik (atau mengambil polinomial dengan derajat serendah mungkin), selama algoritme memungkinkan kesalahan tidak melebihi setengah nilai bit terakhir mantissa. Dalam literatur, ini ditulis sebagai 0,5ulp (ulp = unit di tempat terakhir ).



Misalnya, jika kita berbicara tentang bilangan x tipe float dalam interval (0,5; 1), nilai ulp = 2 ^-23 . Pada interval (1; 2) ulp = 2 ^-22 . Dengan kata lain, jika x berada pada interval (0; 1) maka 2 ^xakan berada pada interval (1,2), dan untuk memastikan akurasi 0,5ulp, Anda perlu, secara kasar, untuk memilih EPS = 2 ^-23 (jadi kami akan menunjukkan konstanta "epsilon", pada orang umum disebut "kesalahan", atau "akurasi", kepada siapa sesukamu, tolong jangan cari kesalahan).



Untuk kalkulasi yang diterapkan, ini sudah cukup, tetapi fakta bahwa bit terakhir mungkin tidak sesuai dengan hasil absolut tidak penting bagi hampir 100% pemrogram, karena tidak penting bagi mereka apa bit itu, tetapi apa akurasinya.



Bagi yang belum paham, saya akan memberikan contoh pada sistem bilangan desimal. Berikut adalah dua angka: 1.999999 dan 2.0. Katakanlah yang pertama adalah apa yang diterima programmer, dan yang kedua adalah standar dari apa yang seharusnya diperoleh jika kita memiliki kemungkinan yang tidak terbatas. Perbedaan di antara mereka hanya sepersejuta, artinya jawaban dihitung dengan kesalahan EPS = 10 ^-6 . Namun, tidak ada satupun angka yang benar dalam jawaban ini. Apa itu buruk? Tidak, dari sudut pandang program aplikasi, ini ungu, pemrogram akan membulatkan jawaban, katakanlah, ke dua tempat desimal dan akan menerima 2,00 (misalnya, tentang mata uang, $ 2,00), dia tidak perlu lebih, tetapi fakta bahwa dia letakkan EPS = 10 ^-6 di program saya , kemudian bagus, mengambil margin untuk kesalahan perhitungan menengah dan menyelesaikan masalah dengan benar.



Dengan kata lain, jangan bingung: presisi dan jumlah bit (atau digit) yang benar adalah dua hal yang berbeda. Mereka yang membutuhkan akurasi (ini hampir 100% programmer), masalah yang dibahas sama sekali tidak menyangkut mereka. Siapa pun yang membutuhkan urutan bit untuk mencocokkan referensi yang dibulatkan dengan benar sangat mengkhawatirkan masalah ini, misalnya, pengembang pustaka fungsi dasar. Meskipun demikian, penting bagi setiap orang untuk mengetahui tentang ini untuk pengembangan umum.



Izinkan saya mengingatkan Anda bahwa ini adalah arah pertama dari masalah: bit terakhir dari jawaban mungkin salah karena ini adalah solusi yang disengaja. Hal utama adalah menjaga akurasi 0,5ulp (atau lebih tinggi). Oleh karena itu, algoritma numerik dipilih hanya dari kondisi ini, jika saja ia bekerja dengan sangat cepat. Dalam hal ini, Standar memungkinkan penerapan fungsi dasar tanpa pembulatan yang benar dari bit terakhir. Saya mengutip [1, bagian 12.1] (Inggris):

Versi 1985 dari Standar IEEE 754 untuk Aritmatika Titik Mengambang tidak menentukan apapun tentang fungsi dasar. Ini karena telah dipercaya selama bertahun-tahun bahwa fungsi yang dibulatkan dengan benar akan terlalu lambat setidaknya untuk beberapa argumen masukan. Situasi berubah sejak saat itu dan versi standar 2008 merekomendasikan (namun tidak memerlukan) bahwa beberapa fungsi dibulatkan dengan benar.

Berikut ini adalah fungsi yang direkomendasikan tetapi tidak perlu dibulatkan dengan benar:

Alasan kedua

Akhirnya kita sampai pada topik pembicaraan: Dilema Pembuat Tabel (disingkat TMD). Saya tidak dapat menerjemahkan namanya secara memadai ke dalam bahasa Rusia, hal ini diperkenalkan oleh William Kahan (pendiri IEEE-754) dalam artikel [2]. Mungkin jika Anda membaca artikelnya, Anda akan mengerti mengapa namanya persis seperti itu. Singkatnya, inti dari dilema ini adalah bahwa kita perlu mendapatkan pembulatan yang benar-benar akurat dari fungsi z = f (x), seolah-olah kita memiliki catatan bit tak terhingga dari hasil perhitungan sempurna z yang kita miliki. Tetapi jelas bagi semua orang bahwa kita tidak bisa mendapatkan urutan yang tak terbatas. Berapa banyak bit yang harus diambil? Di atas, saya menunjukkan contoh ketika kita perlu melihat 40 bit hasil untuk mendapatkan setidaknya 2 bit yang benar setelah pembulatan. Dan inti dari masalah TMD adalah kita tidak mengetahuinya terlebih dahulu, terserah berapa banyak bit untuk menghitung nilai z untuk mendapatkan sebanyak mungkin bit yang benar setelah pembulatan sesuai kebutuhan. Bagaimana jika ada seratus atau seribu? Kami tidak tahu sebelumnya!



Misalnya, seperti yang saya katakan, untuk fungsi 2 ^x , untuk tipe data float, di mana bagian pecahan mantissa hanya memiliki 23 bit, kita perlu melakukan perhitungan dengan akurasi ^2-54 sehingga pembulatan terjadi dengan benar untuk semua argumen x yang mungkin tanpa kecuali. Tidak sulit untuk mendapatkan perkiraan ini dengan pencarian menyeluruh, tetapi untuk sebagian besar fungsi lainnya, terutama untuk tipe ganda atau ganda panjang (letakkan "kelas" jika Anda tahu apa itu), perkiraan semacam itu tidak diketahui .



Mari kita sudah memahami mengapa ini terjadi. Sengaja saya berikan contoh pertama pada artikel ini dengan tipe data float dan meminta Anda untuk mengingatnya, karena pada tipe ini hanya ada 32 bit dan akan lebih mudah untuk melihatnya, pada tipe data lain situasinya serupa.



Kami mulai dengan angka x = 0,00296957581304013729095458984375, ini adalah angka yang benar-benar dapat direpresentasikan dalam tipe data float, yaitu, ditulis sehingga dapat dikonversi ke sistem float biner tanpa pembulatan. Kami menghitung 2 ^x , dan jika kami memiliki kalkulator dengan presisi tak terbatas, maka kami harus mendapatkan (Jadi Anda dapat memeriksa saya, perhitungan dilakukan di sistem online WolframAlpha ):



1.0020604729652405753669743044108123031635398201893943954577320057 ...



Mari terjemahkan angka ini ke dalam biner, katakanlah 64 bit sudah cukup:



1.00000000100001110000100 1 000000000000000000000000000001101111101



Bit pembulatan (bit ke-24 setelah titik desimal) digarisbawahi. Pertanyaan: pembulatan ke mana? Atas atau bawah? Jelas, Anda tahu ini karena Anda melihat cukup banyak bagian dan Anda dapat membuat keputusan. Tapi perhatikan baik-baik ...



Setelah bit pembulatan, kami memiliki 29 angka nol. Ini berarti bahwa kita berada sangat, sangat dekat dengan tengah antara dua bilangan terdekat dan cukup bergerak ke bawah sedikit, karena arah pembulatan akan berubah. Tetapi pertanyaannya adalah: di manakah pergeseran ini? Algoritme numerik dapat secara berurutan, selangkah demi selangkah, mendekati nilai yang tepat dari sisi yang berbeda, dan sampai kita melewati semua 29 angka nol ini dan mencapai akurasi yang melebihi nilai nol terakhir di "lokomotif" ini, kita tidak akan mengetahui arah pembulatan ... Bagaimana jika sebenarnya jawaban yang benar adalah:



1.00000000100001110000100 0 11111111111111111111111111111?



Kemudian pembulatannya akan turun.



Kami tidak mengetahui hal ini hingga presisi mencapai bit ke-54 setelah koma desimal. Ketika bit ke-54 diketahui dengan tepat, kita akan tahu dengan tepat yang mana dari dua bilangan terdekat yang mendekati kita. Angka-angka seperti itu disebut titik yang paling sulit untuk dibulatkan [1, bagian 12.3] (titik kritis untuk pembulatan), dan angka 54 disebut kekerasan-ke-putaran, dan dilambangkan dengan huruf m pada buku yang dikutip.



Kompleksitas pembulatan (m) adalah jumlah bit minimum yang diperlukan untuk memastikan bahwa untuk semua argumen dari fungsi tertentu f (x) dan untuk rentang yang telah dipilih sebelumnya, fungsi f (x) dibulatkan dengan benar ke bit terakhir (untuk mode pembulatan yang berbeda mungkin ada yang berbeda nilai m). Dengan kata lain, untuk tipe data float dan untuk argumen x dari kisaran (0; 1) untuk mode pembulatan "genap terdekat", kompleksitas pembulatan adalah m = 54. Ini berarti bahwa untuk benar-benar semua x dari interval (0; 1) kita dapat memasukkan ke dalam algoritme dengan presisi yang sama ESP = ^2-54 , dan semua hasil akan dibulatkan dengan benar menjadi 23 bit setelah titik desimal biner.



Faktanya, beberapa algoritme dapat memberikan hasil yang tepat dan berdasarkan 53 dan bahkan 52 bit, brute force menunjukkan hal ini, tetapi secara teoritis, Anda membutuhkan tepat 54. Jika bukan karena kemungkinan untuk mengeluarkan brute force, kami tidak dapat "menipu" dan simpan beberapa bit, seperti yang saya lakukan di program brute force di atas. Saya mengambil polinomial dengan derajat yang lebih rendah dari yang seharusnya, tetapi tetap berhasil, hanya karena saya beruntung.



Jadi, terlepas dari mode pembulatan, kita memiliki dua kemungkinan situasi: "lokomotif uap" nol muncul di area pembulatan, atau "lokomotif uap" salah satunya. Tugas dari algoritma yang benar untuk menghitung fungsi transendental f (x) adalah untuk memperbaiki nilai dari fungsi ini sampai keakuratannya melebihi nilai dari bit terakhir dari "lokomotif uap" ini, dan sampai menjadi jelas jelas bahwa sebagai hasil dari fluktuasi berikutnya dari algoritma numerik untuk menghitung f (x) nol tidak akan berubah menjadi satu, atau sebaliknya. Segera setelah semuanya stabil, dan algoritme telah mencapai akurasi yang melampaui batas "lokomotif uap", maka kita dapat membulatkan seolah-olah kita memiliki jumlah bit yang tak terbatas. Dan pembulatan ini akan dilakukan dengan bit terakhir yang benar. Tetapi bagaimana ini bisa dicapai?

"Kruk"

Seperti disebutkan, masalah utamanya adalah mendapatkan algoritma untuk mengatasi lokomotif nol atau lokomotif yang muncul segera setelah bit pembulatan. Ketika lokomotif diatasi dan kita melihatnya secara keseluruhan, maka ini setara dengan fakta bahwa nol atau satu ini sudah dihitung dengan tepat , dan kita sudah tahu persis ke arah mana pembulatan akan terjadi. Namun jika kita tidak mengetahui panjang lokomotif tersebut, lalu bagaimana kita dapat merancang suatu algoritma?

"Kruk" pertama

Bagi pembaca, jawabannya jelas: ambil aritmatika dengan ketepatan tak terbatas dan sengaja masukkan jumlah bit yang berlebihan, dan jika tidak cukup, masukkan lagi dan hitung ulang. Secara umum, itu benar. Ini dilakukan ketika kecepatan dan sumber daya komputer tidak memainkan peran khusus. Pendekatan ini memiliki nama: Strategi bertingkat Ziv [1, bagian 12.3]. Esensinya sangat sederhana. Algoritme harus mendukung penghitungan pada beberapa tingkatan: penghitungan awal yang cepat (dalam banyak kasus ternyata final), penghitungan yang lebih lambat tetapi lebih akurat (menghemat dalam banyak kasus kritis), bahkan lebih lambat, tetapi bahkan penghitungan yang lebih akurat (ketika benar-benar "buruk "Harus) dan seterusnya.



Dalam sebagian besar kasus, cukup mengambil akurasi sedikit lebih tinggi dari 0,5ulp, tetapi jika "lokomotif" muncul, maka kami meningkatkannya. Selama "lokomotif uap" tetap ada, kami meningkatkan keakuratannya hingga jelas sekali bahwa fluktuasi lebih lanjut dari metode numerik tidak akan memengaruhi "lokomotif uap" ini. Jadi, misalnya, dalam kasus kita, jika kita telah mencapai ESP = 2 ^-54 , maka pada posisi ke-54 sebuah unit muncul, yang seolah-olah "melindungi" lokomotif dari nol dan menjamin bahwa tidak akan ada lagi pengurangan nilai yang lebih besar dari atau sama dengan 2 ^-53 dan nol tidak akan berubah menjadi satu, menyeret bit pembulatan ke nol dengannya.



Itu adalah presentasi sains yang populer, semua sama dengan tes pembulatan Ziv, di mana ditunjukkan seberapa cepat, dalam satu langkah, untuk memeriksa apakah kita telah mencapai akurasi yang diinginkan, dapat dibaca di [1, Bab 12], atau di [3, Bagian 10.5].



Masalah dengan pendekatan ini jelas. Algoritme perlu dirancang untuk menghitung setiap fungsi transendental f (x) sehingga dapat meningkatkan akurasi penghitungan selama bagian tersebut. Untuk implementasi perangkat lunak, ini masih tidak terlalu menakutkan, misalnya, metode Newton memungkinkan, secara kasar, untuk menggandakan jumlah bit tepat setelah titik desimal pada setiap iterasi. Anda dapat menggandakan hingga menjadi "cukup", meskipun ini adalah proses yang agak memakan waktu, saya harus mengakui bahwa metode Newton tidak selalu dapat dibenarkan, karena memerlukan penghitungan fungsi invers f ^-1(x), yang dalam beberapa kasus mungkin tidak sesederhana menghitung f (x) itu sendiri. Untuk implementasi perangkat keras, "strategi Ziva" sama sekali tidak cocok. Algoritme, yang terpasang pada prosesor, harus melakukan serangkaian tindakan dengan jumlah bit yang telah ditentukan sebelumnya, dan ini cukup bermasalah untuk diterapkan jika kita tidak mengetahui angka ini sebelumnya. Mengambil stok, mengambil persediaan? Berapa banyak?



Pendekatan probabilistik untuk memecahkan masalah [1, Bagian 12.6] memungkinkan kita untuk memperkirakan nilai m (ingat, ini adalah jumlah bit, yang cukup untuk pembulatan yang benar). Ternyata panjang "lokomotif" dalam arti probabilistik sedikit lebih besar dari panjang mantissa angka tersebut. Jadi, dalam kebanyakan kasus, mantissa akan cukup untuk mengambil m sedikit lebih dari dua kali lipat nilai mantissa, dan hanya dalam kasus yang sangat jarang mantissa perlu mengambil lebih banyak lagi. Saya mengutip penulis karya ini: "kami menyimpulkan bahwa dalam praktiknya, m harus sedikit lebih besar dari 2p" (mereka memiliki p - panjang mantissa bersama dengan bagian bilangan bulat, yaitu, p = 24 untuk float). Lebih lanjut dalam teks, mereka menunjukkan bahwa kemungkinan kesalahan dengan strategi semacam itu mendekati nol, tetapi masih positif, dan ini dikonfirmasi oleh eksperimen.



Namun demikian, masih terdapat kasus dimana nilai m harus diambil lebih banyak, dan kasus terburuk tidak diketahui sebelumnya. Perkiraan teoritis dari kasus terburuk ada [1, bagian 12.7.2], tetapi mereka menghasilkan jutaan bit yang tidak terpikirkan, yang tidak bagus. Berikut adalah tabel dari karya yang dikutip (ini untuk fungsi exp (x) pada interval dari -ln (2) hingga ln (2)):

"Kruk" kedua

Dalam praktiknya, m tidak akan terlalu besar. Dan untuk menentukan kasus terburuk, "kruk" kedua diterapkan, yang disebut "pra-perhitungan lengkap". Untuk tipe data float (32 bits), jika fungsi f memiliki satu argumen (x), maka kita dapat dengan mudah "menjalankan" semua kemungkinan nilai x. Masalahnya hanya akan muncul dengan fungsi yang memiliki lebih dari satu argumen (di antaranya pow (x, y)), yang tidak dapat kita pikirkan seperti itu. Setelah memeriksa semua kemungkinan nilai x, kami menghitung m konstanta kami untuk setiap fungsi f (x) dan untuk setiap mode pembulatan. Kemudian algoritma perhitungan yang perlu diimplementasikan pada perangkat keras dirancang untuk memberikan akurasi sebesar 2 ^-m . Kemudian pembulatan f (x) dijamin benar di semua kasus.



Untuk tipe ganda (64 bit), pencacahan sederhana hampir tidak mungkin. Namun, mereka memilah! Tapi bagaimana caranya? Jawabannya diberikan dalam [4]. Saya akan memberi tahu Anda tentang itu dengan sangat singkat.



Domain dari fungsi f (x) dibagi menjadi beberapa segmen yang sangat kecil sehingga di dalam setiap segmen dimungkinkan untuk mengganti f (x) dengan fungsi linier bentuk b-ax (koefisien a dan b, tentu saja, berbeda untuk segmen yang berbeda). Ukuran segmen-segmen ini dihitung secara analitis sehingga fungsi linier semacam itu memang hampir tidak dapat dibedakan dari aslinya di setiap segmen.



Kemudian, setelah beberapa operasi penskalaan dan pergeseran, kita sampai pada masalah berikut: dapatkah garis b-ax lurus "cukup dekat" ke titik integer?



Ternyata relatif mudah untuk memberikan jawaban ya atau tidak. Yaitu, "ya" - jika titik yang berpotensi berbahaya dekat dengan garis lurus, dan "tidak" - jika tidak ada titik seperti itu, pada prinsipnya, dapat mendekati garis. Keindahan dari metode ini adalah bahwa jawaban "tidak" dalam praktik diperoleh dalam sebagian besar kasus, dan jawaban "ya", yang jarang diperoleh, memaksa Anda untuk menelusuri segmen dengan kekuatan kasar untuk menentukan poin spesifik mana yang kritis.



Namun demikian, iterasi atas argumen f (x) berkurang berkali-kali dan memungkinkan untuk mendeteksi titik putus untuk bilangan seperti ganda (binary64) dan panjang ganda (80 bit!). Ini dilakukan pada superkomputer dan, tentu saja, kartu video ... di waktu luang Anda dari penambangan. Namun, belum ada yang tahu apa yang harus dilakukan dengan tipe data binary128. Izinkan saya mengingatkan Anda bahwa bagian pecahan dari mantissa dari bilangan tersebut adalah 112 bit . Oleh karena itu, dalam literatur asing mengenai hal ini selama ini, hanya dapat ditemukan penalaran semi filosofis yang dimulai dengan “we hope…” (“we hope…”).



Detail metode, yang memungkinkan Anda untuk dengan cepat menentukan jalur di dekat titik bilangan bulat, tidak tepat di sini. Bagi mereka yang ingin mempelajari prosesnya lebih teliti, saya sarankan untuk melihat ke masalah menemukan jarak antara garis lurus dan Z ² , misalnya, di artikel [5]. Ini menggambarkan algoritme yang ditingkatkan, yang dalam proses pembuatannya menyerupai algoritme Euclid yang terkenal untuk menemukan pembagi persekutuan terbesar. Saya akan memberikan gambar yang sama dari [4] dan [5], yang menunjukkan transformasi masalah lebih lanjut:







Ada tabel ekstensif yang berisi kasus terburuk pembulatan pada interval berbeda untuk setiap fungsi transendental. Mereka ditemukan di [1 bagian 12.8.4] dan di [3, bagian 10.5.3.2], serta di artikel terpisah, misalnya, di [6].



Saya akan memberikan beberapa contoh dengan mengambil baris acak dari tabel tersebut. Saya tekankan bahwa ini bukan kasus terburuk untuk semua x, tetapi hanya untuk beberapa interval kecil, lihat sumbernya jika Anda tertarik.







Bagaimana cara membaca tabel? Nilai x ditentukan dalam notasi ganda floating point heksadesimal. Pertama, seperti yang diharapkan, ada yang terdepan, kemudian 52 bit dari bagian pecahan dari mantissa dan huruf P. Huruf ini berarti "kalikan dengan dua pangkat" diikuti dengan derajat. Misalnya, P-23 berarti mantissa yang ditentukan perlu dikalikan dengan 2 ^-23 .



Selanjutnya, bayangkan bahwa fungsi f (x) dihitung dengan presisi tak hingga dan 53 bit pertama dipotong darinya (tanpa pembulatan!). Ini adalah 53 bit ini (salah satunya hingga koma) yang ditunjukkan di kolom f (x). Bit selanjutnya ditunjukkan di kolom terakhir. Tanda "derajat" dari urutan bit di kolom terakhir berarti jumlah pengulangan bit, yaitu, misalnya, 10 ⁵³1011 artinya pertama ada bit yang sama dengan 1, lalu 53 nol dan kemudian 1011. Kemudian tiga titik, yang berarti bahwa kita secara umum tidak membutuhkan bit yang lain sama sekali.



Lebih lanjut, ini adalah masalah teknologi - kita mengetahui kasus terburuk untuk setiap interval dari fungsi yang diambil secara terpisah dan kita dapat memilih untuk interval ini perkiraan sedemikian sehingga mencakup kasus terburuk dengan akurasinya. Dengan hanya beberapa tahun komputasi superkomputer, dimungkinkan untuk membuat implementasi perangkat keras fungsi dasar yang cepat dan akurat. Masalahnya kecil: tetap mengajarkan setidaknya pengembang kompiler untuk menggunakan tabel ini.

Mengapa ini dibutuhkan?

Pertanyaan bagus! Lagi pula, saya telah berulang kali berbicara di atas bahwa hampir 100% programmer tidak perlu mengetahui fungsi dasar dengan akurasi hingga bit terakhir yang dibulatkan dengan benar (seringkali mereka bahkan tidak memerlukan setengah dari bit), mengapa ilmuwan menggerakkan superkomputer dan menyusun tabel untuk memecahkan masalah "tidak berguna"?



Pertama, tantangannya sangat mendasar. Agak menarik untuk tidak mendapatkan pembulatan yang tepat demi pembulatan yang akurat, tetapi pada prinsipnya untuk memahami bagaimana masalah yang menarik ini dapat diselesaikan, rahasia apa dari matematika komputasi yang akan diungkapkan oleh solusinya kepada kita? Bagaimana rahasia ini dapat digunakan dalam tugas lain? Ilmu-ilmu fundamental - memang seperti itu, Anda dapat melakukan semacam "omong kosong" selama beberapa dekade, dan kemudian seratus tahun kemudian, berkat "omong kosong" ini, terobosan ilmiah terjadi di beberapa bidang lain.



Kedua, masalah portabilitas kode. Jika suatu fungsi mampu menangani bit terakhir dari hasil dengan cara apa pun yang diinginkannya, itu berarti bahwa pada platform yang berbeda dan pada kompiler yang berbeda, hasil yang sedikit berbeda dapat diperoleh (bahkan jika dalam kesalahan yang ditentukan). Dalam beberapa kasus, ini tidak penting, tetapi dalam beberapa kasus mungkin signifikan, terutama ketika program memiliki kesalahan yang muncul di satu platform, tetapi tidak muncul di platform lain justru karena bit yang berbeda dari hasil. Tetapi mengapa saya menjelaskan kepada Anda sakit kepala terkenal yang terkait dengan perilaku program yang berbeda? Anda tahu semua ini tanpa saya. Akan sangat bagus untuk memiliki sistem matematika yang bekerja persis sama di semua platform, tidak peduli seberapa terkompilasinya. Itulah yang perlu Anda lakukan dengan benar mengakhiri bagian terakhir.

Daftar sumber

[1] Jean-Michel Muller, “Elementary Functions: Algorithms and Implementation”, 2016



[2] William Kahan, “ A Logarithm Too Clever by Half ”, 2004



[3] Jean-Michel Muller, “Handbook of floating-point arithmetic” , 2018



[4] Vincent Lefèvre, Jean-Michel Muller, “Toward Correctly Rounded Transcendentals”, IEEE TRANSACTIONS ON COMPUTERS, VOL. 47, TIDAK. 11, NOVEMBER 1998. hlm. 1235-1243



[5] Vincent Lefèvre. “Hasil Baru pada Jarak Antara Segmen dan Z ² ”. Aplikasi ke Pembulatan Tepat. 17th IEEE Symposium on Computer Arithmetic - Arith'17, Jun 2005, Cape Cod, MA,

Amerika Serikat. hlm. 68-75



[6] Vincent Lefèvre, Jean-Michel Muller, “Kasus Terburuk untuk Pembulatan yang Benar dari Fungsi Dasar dalam Presisi Ganda”, Rapport de recherche (INSTITUT NA TIONAL DE RECHERCHE EN INFORMA TIQUE ET EN AUTOMA TIQUE) n˚4044 - Novembre 2000 - 19 halaman.