Penghitungan yang akurat dan cepat untuk bilangan floating point menggunakan contoh fungsi sinus. Bagian 3: titik tetap

Kami melanjutkan siklus perkuliahan ( bagian 1 dan bagian 2 ). Pada bagian 2, kami melihat apa yang ada di dalam perpustakaan libm dan dalam pekerjaan ini kami akan mencoba sedikit mengubah fungsi do_sin untuk meningkatkan akurasi dan kecepatannya. Saya akan mengutip fungsi ini lagi do_sin ):







Seperti yang ditunjukkan pada artikel sebelumnya, bagian 132-145. Dieksekusi untuk x dalam kisaran [0,126, 0,855469]. Baik. Mari kita coba menulis fungsi yang, dalam batasan yang diberikan, akan lebih akurat dan, mungkin, lebih cepat.



Cara kita menggunakannya cukup jelas. Akurasi penghitungan perlu diperluas untuk memasukkan lebih banyak tempat desimal. Solusi yang jelas adalah memilih tipe ganda panjang, menghitung di dalamnya, dan kemudian mengubahnya kembali. Dari segi akurasi, solusinya harus bagus, tapi dari segi performa, mungkin ada masalah. Bagaimanapun, long double adalah jenis data yang agak eksotis dan dukungannya dalam prosesor modern bukanlah prioritas. Pada instruksi SSE / AVX x86_64 dengan tipe data ini tidak berfungsi. Koprosesor matematis akan "terhempas".



Lalu apa yang harus Anda pilih? Mari kita lihat lebih dekat argumen dan batas fungsinya.



Mereka berada di wilayah 1.0. Itu. sebenarnya, kita tidak membutuhkan floating point. Mari gunakan bilangan bulat 64-bit saat menghitung fungsi. Ini akan memberi kita 10-11 bit tambahan ke presisi asli. Mari kita cari tahu cara menggunakan angka-angka ini. Angka dalam format ini direpresentasikan sebagai a / d , di mana a adalah bilangan bulat, dan d adalah pembagi yang kita pilih konstanta untuk semua variabel dan simpan "di memori kita", dan bukan di memori komputer. Di bawah ini adalah beberapa operasi untuk nomor tersebut:

$$

$$\frac{c}{d}=\frac{a}{d}\pm\frac{b}{d}=\frac{a\pm b}{d}\\ \frac{c}{d}=\frac{a}{d}\cdot\frac{b}{d}=\frac{a \cdot b}{d^2}\\ \frac{c}{d}=\frac{a}{d}\cdot x=\frac{a \cdot x}{d}$$

Seperti yang Anda lihat, tidak ada yang rumit tentang itu. Rumus terakhir menunjukkan perkalian dengan bilangan bulat apa pun. Perhatikan juga hal yang cukup jelas bahwa hasil perkalian dua variabel integer unsigned dengan ukuran N lebih sering merupakan jumlah ukuran hingga 2 * N inklusif. Penambahan dapat menyebabkan luapan hingga 1 bit ekstra.



Mari kita coba memilih pembagi d . Jelas, di dunia biner, yang terbaik adalah memilihnya sebagai pangkat dua, agar tidak membagi, tetapi hanya memindahkan register, misalnya. Kekuatan apa dari dua yang harus Anda pilih? Temukan petunjuknya di instruksi mesin perkalian. Misalnya, instruksi standar MUL dalam sistem x86 mengalikan 2 register dan menulis hasilnya ke 2 register juga, di mana 1 register adalah "bagian atas" dari hasil, dan yang kedua adalah bagian bawah.



Misalnya, jika kita memiliki 2 bilangan 64-bit, maka hasilnya adalah bilangan 128-bit yang ditulis menjadi dua register 64-bit. Sebut saja RH - huruf "atas", dan RL - huruf "kecil" ¹ . Kemudian secara matematis hasilnya bisa ditulis$$$R=R_H \cdot 2^{64}+R_L$... Sekarang kita menggunakan rumus di atas dan menulis perkaliannya$$$d=2^{-64}$

$$

$$\frac{c}{d}=\frac{a}{2^{64}}\cdot\frac{b}{2^{64}}=\frac{a \cdot b}{2^{128}}=\frac{R_H \cdot 2^{64} + R_L}{2^{128}}=\frac{R_H + R_L \cdot 2^{-64}}{2^{64}}$$

Dan ternyata hasil perkalian dua bilangan fixed-point ini adalah registernya $$$R=R_H$...



Lebih mudah lagi untuk sistem Aarch64. Instruksi "UMULH" mengalikan dua register dan menulis bagian "atas" dari perkalian ke register ke-3.



Baiklah kalau begitu. Kami telah menentukan nomor fixed-point, tetapi masih ada satu masalah. Angka negatif. Dalam deret Taylor, ekspansi dilakukan dengan tanda variabel. Untuk mengatasi masalah ini, kami mengubah rumus untuk menghitung polinomial dengan metode Goner, menjadi bentuk berikut:

$$

$$\sin(x)\approx x(1-x^2(1/3!-x^2(1/5!-x^2(1/7!-x^2\cdot1/9!))))$$

Periksa apakah secara matematis sama persis dengan rumus aslinya. Namun di setiap tanda kurung ada sejumlah formulir$$$1/(2n + 1)! - x^2\cdot(\cdots)$selalu positif. Itu. konversi ini memungkinkan ekspresi dievaluasi sebagai bilangan bulat tak bertanda.

constexpr mynumber toint    = {{0x00000000, 0x43F00000}};  /*  18446744073709551616 = 2^64     */
constexpr mynumber todouble = {{0x00000000, 0x3BF00000}};  /*  ~5.42101086242752217003726400434E-20 = 2^-64     */

double sin_e7(double xd) {
  uint64_t x = xd * toint.x;
  uint64_t xx = mul2(x, x);
  uint64_t res = tsx[19]; 
  for(int i = 17; i >= 3; i -= 2) {
    res = tsx[i] - mul2(res, xx);
  }
  res = mul2(res, xx);
  res = x - mul2(x, res);
  return res * todouble.x;
}

Dimana $$$tsx[i]=1/i!$dalam format titik tetap. Kali ini, untuk kenyamanan, saya memposting semua kode di fast_sine github , menyingkirkan quadmath untuk kompatibilitas dengan dentang dan lengan. Dan saya mengubah sedikit metode penghitungan kesalahan.



Perbandingan fungsi sinus standar dan fungsi titik tetap diberikan dalam dua tabel di bawah ini. Tabel pertama menunjukkan keakuratan penghitungan (untuk x86_64 dan ARM) sama sekali. Tabel kedua adalah perbandingan kinerja.







Selama pengujian, nilai sinus "sebenarnya" dihitung menggunakan perpustakaan MPFR... Nilai ULP maksimum dianggap sebagai deviasi maksimum dari nilai "sebenarnya". Persentase kesalahan - jumlah kasus ketika nilai yang dihitung dari fungsi sinus oleh kami atau libm tidak cocok dengan nilai sinus yang dibulatkan menjadi dua kali lipat. Nilai deviasi rata-rata menunjukkan "arah" kesalahan kalkulasi: perkiraan terlalu tinggi atau terlalu rendah dari nilai. Seperti yang Anda lihat dari tabel, fungsi kita cenderung melebih-lebihkan nilai sinus. Ini bisa diperbaiki! Siapa bilang nilai tsx harus persis sama dengan koefisien deret Taylor. Ide yang agak jelas menyarankan dirinya untuk memvariasikan nilai koefisien untuk meningkatkan keakuratan pendekatan dan menghilangkan komponen kesalahan yang konstan. Sangat sulit untuk membuat variasi seperti itu dengan benar. Tapi kita bisa mencoba. Mari kita ambil contohNilai ke-4 dari larik koefisien tsx (tsx [3]) dan ubah angka terakhir a menjadi f. Mari restart program dan lihat akurasinya (sin_e7a). Lihat, ini sedikit, tapi bertambah! Kami menambahkan metode ini ke celengan kami.



Sekarang mari kita lihat performanya. Untuk pengujian, saya mengambil apa yang ada di tangan i5 mobile dan raspberry keempat yang sedikit di-overclock (Raspberry PI 4 8GB), GCC10 dari distribusi Ubuntu 20.04 x64 untuk kedua sistem.







Saya tidak berpura-pura lebih akurat dalam pengukuran ini. Variasi beberapa puluh persen dimungkinkan tergantung pada beban prosesor. Kesimpulan utama bisa dibuat seperti ini. Beralih ke aritmatika integer tidak memberikan peningkatan kinerja pada prosesor modern ² . Jumlah transistor yang tak terbayangkan dalam prosesor modern memungkinkan untuk melakukan perhitungan kompleks dengan cepat. Tapi, menurut saya, pada prosesor seperti Intel Atom, serta pada pengontrol yang lemah, pendekatan ini dapat memberikan peningkatan kinerja yang signifikan. Bagaimana menurut anda?



Meskipun pendekatan ini memberikan peningkatan akurasi, peningkatan akurasi ini tampaknya lebih menarik daripada berguna. Dalam hal kinerja, pendekatan ini dapat ditemukan sendiri, misalnya, di IoT. Tetapi untuk komputasi kinerja tinggi, ini bukan lagi arus utama. Dalam dunia sekarang ini SSE / AVX / CUDA lebih suka menggunakan penghitungan fungsi paralel. Dan dalam aritmatika floating point. Tidak ada analog paralel dari fungsi MUL. Fungsi itu sendiri lebih merupakan penghargaan terhadap tradisi.



Di bab berikutnya, saya akan menjelaskan bagaimana Anda dapat menggunakan AVX secara efektif untuk kalkulasi. Sekali lagi, mari kita masuk ke kode libm dan mencoba memperbaikinya.



¹ Tidak ada register dengan nama seperti itu di prosesor yang saya kenal. Nama-nama telah dipilih misalnya.

²Perlu dicatat di sini bahwa ARM saya dilengkapi dengan versi terbaru dari coprocessor matematika. Jika prosesor meniru kalkulasi floating point, hasilnya bisa sangat berbeda.