👨‍👨‍👦 🌑 👨🏾‍🚒 Penghitungan yang akurat dan cepat untuk bilangan floating point menggunakan contoh fungsi sinus. Pendahuluan dan bagian 1 👃🏻 🧑🏿‍🤝‍🧑🏻 🏏

Baca dengan cermat artikel yang sangat bagus dari ArtemKaravaevpada penambahan angka floating point. Topiknya sangat menarik dan saya ingin melanjutkannya dan menunjukkan dengan contoh bagaimana bekerja dengan bilangan floating point dalam praktek. Mari kita ambil pustaka glibc GNU (libm) sebagai referensi. Dan agar artikel tidak terlalu membosankan, mari kita tambahkan komponen kompetitif: kita akan mencoba tidak hanya mengulang, tetapi juga meningkatkan kode perpustakaan, membuatnya lebih cepat / akurat.

Sebagai contoh, saya telah memilih fungsi sinus trigonometri. Ini adalah fungsi yang tersebar luas yang matematikanya terkenal dari sekolah dan universitas. Pada saat yang sama, ketika diimplementasikan, akan ada banyak contoh nyata dari pekerjaan yang "benar" dengan angka. Saya akan menggunakan double sebagai angka floating point.

Dalam seri artikel ini, banyak hal yang direncanakan, mulai dari matematika hingga kode mesin dan opsi kompilator. Bahasa penulisan artikel adalah C ++, tetapi tanpa "embel-embel". Berbeda dengan bahasa C, contoh yang berfungsi akan lebih mudah dibaca bahkan untuk orang yang tidak terbiasa dengan bahasa tersebut dan menggunakan lebih sedikit baris.

Artikel akan ditulis dengan metode pencelupan. Subtugas akan dibahas, yang kemudian akan digabungkan dalam satu solusi untuk masalah tersebut.

Dekomposisi sinus dalam deret Taylor.

Fungsi sinus berkembang menjadi deret Taylor tak terhingga.

\sin (x) = x - \frac{x^{3}}{3!} + \frac{x^{5}}{5!} - \frac{x^{7}}{7!} + \frac{x^{9}}{9!} - \dots

$\sin(x)=x-{\frac {x^{3}}{3!}}+{\frac {x^{5}}{5!}}-{\frac {x^{7}}{7!}}+{\frac {x^{9}}{9!}}-\cdots$

Jelas bahwa kita tidak dapat menghitung deret tak hingga, kecuali untuk kasus ketika ada rumus analitik untuk jumlah tak hingga. Tapi ini bukan kasus kami))) Misalkan kita ingin menghitung sinus dalam interval

[0, \frac{π}{2}]

$[0, \frac{\pi}{2}]$ ... Kami akan membahas pekerjaan dengan interval secara lebih rinci di Bagian 3. Mengetahui itu

\sin (\frac{π}{2}) = 1

$\sin(\frac{\pi}{2})=1$ Perkirakan, temukan suku pertama yang bisa dibuang berdasarkan kondisi itu

\frac{(π / 2)^{n}}{n!} < e

$\frac{(\pi/2)^n}{n!}<e$ dimana

e

$e$ ini adalah perbedaan antara angka 1 dan angka terkecil yang lebih besar dari 1. Secara kasar, ini adalah bit terakhir dari mantissa ( wiki ). Persamaan ini lebih mudah diselesaikan dengan gaya brute force. Untuk

e \approx 2.22 \times 10^{- 16}

$e \approx 2.22\times10^{-16}$ ... saya mengatur

n = 23

$n=23$ sudah bisa dijatuhkan. Pilihan yang benar dari jumlah istilah akan dibahas di salah satu bagian selanjutnya, jadi untuk hari ini kita akan "aman" dan mengambil persyaratan hingga

n = 25

$n=25$ inklusif.

Istilah terakhir kira-kira 10.000 kali lebih kecil dari

e

$e$ ...

Solusi paling sederhana

Tangan sudah gatal, kami menulis:

Teks lengkap dari program untuk pengujian

#include <iostream>
#include <iomanip>
#include <cmath>
#include <array>
#include <bitset>
#include <quadmath.h>
//      clang
//#include "/usr/lib/gcc/x86_64-linux-gnu/10/include/quadmath.h"
#include <numeric>
#include <limits>
#include <vector>

#include <boost/timer/timer.hpp>
#include <boost/math/special_functions/factorials.hpp>

namespace bm = boost::math;

using namespace std;

typedef union { uint32_t i[2]; double x; } mynumber;

array<double, 25> fc;

double sin_e1(double x) {
  double result = 0;
  int sign = 1;
  for(int i = 1; i < 25; i += 2) {
    result += sign * pow(x, i) / bm::unchecked_factorial<double>(i);
    sign = -sign;
  }
  return result;
}

double sin_e2(double x) {
  double result = 0;
  int sign = 1;
  double xx = x * x;
  double pw = x;
  double fti = 1.0;
  for(int i = 1; i < 25; i += 2) {
    fti /= i;
    result += sign * pw * fti;
    fti /= ( i + 1 );
    sign = -sign;
    pw  *= xx;
  }
  return result;
}

double sin_e3(double x) {
  double result = 0;
  for(int i = 25; i >= 1; i -= 2) {
    result += (((i - 1) % 4 == 0) ? 1 : -1 ) * pow(x, i) / bm::unchecked_factorial<double>(i);
  }
  return result;
}

double sin_e4(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 1; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x * res;
}

double sin_e5(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 3; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x + x * xx * res;
}

inline
double fsin(double x) {
  double result;
  asm ("fsin" :"=t" (result) : "0" (x));
  return result;
}

#define SIN(a) fsin(a)
//#define SIN(a) sin(a)
//#define SIN(a) sin_e5(a)
// ^^     . ^^

int main() {

  __uint128_t ft = 1;
  fc[1] = 1.0; //3 * 5;
  for(int i = 2; i < fc.size(); i++) {
    ft *= i;
    // factorial with sign for Taylor series
    fc[i] = (1.0 / ft) * (( (i - 2) % 4 < 2) ? -1 : 1);
  }
  vector<double> xv;
  xv.resize(8 * 2000000);
  //      0  M_PI/2
  for (int i = 0; i < xv.size(); i++) {
    //      .
    xv[i] = (M_PI / 2) * i / double(xv.size());
  }

  double res = 0;
  {
    boost::timer::auto_cpu_timer at;
    for(int i = 0; i < xv.size(); i++) {
      res += SIN(xv[i]);
    }
  }

  int co = 0, cn = 0;
  //      .
  __float128 avg = 0.0, div = 0.0;
  for(int i = 0; i < xv.size(); i++) {
    mynumber dold, dnew;
    dold.x = sin(xv[i]);
    dnew.x = SIN(xv[i]);
    __float128 q = sinq(xv[i]); // <= sinq  .
    __float128 dd = __float128(dnew.x) - q;
    //     .
    div += dd * dd;
    avg += dd;
    //  ,         .
    //   ,         .
    if( dold.i[0] != dnew.i[0] || dold.i[1] != dnew.i[1] ) {
      if( fabsq(q - dold.x) <= fabsq(q - dnew.x) )
        co++;
      else
        cn++;
    }
  }
  avg /= xv.size();
  div /= xv.size();

  cout << res << endl;

  //  ,           .
  cout << "Better libm: " <<  co << " / " << xv.size() << "(" << 100.0 * co / xv.size() << "%)" << endl;

  //  ,  ""         .
  cout << "Better new: " <<  cn << " / " << xv.size() << "(" << 100.0 * cn / xv.size() << "%)" << endl;

  //         .
  cout << "  Avg / std new: " << double(avg) << " / " << double(sqrtq( div - avg * avg )) << endl;
  return 0;
}

double sin_e1(double x) {
  double result = 0;
  int sign = 1;
  for(int i = 1; i < 25; i += 2) {
    result += sign * pow(x, i) / bm::factorial<double>(i);
    sign = -sign;
  }
  return result;
}

Bagaimana cara mempercepat program, saya rasa banyak orang yang langsung tahu. Menurut Anda, berapa kali perubahan Anda dapat mempercepat program? Versi yang dioptimalkan dan jawaban atas pertanyaan di bawah spoiler.

Versi program yang dioptimalkan.

double sin_e2(double x) {
  double result = 0;
  int sign = 1;
  double xx = x * x;
  double pw = x;
  double fti = 1.0;
  for(int i = 1; i < 25; i += 2) {
    fti /= i;
    result += sign * pw * fti;
    fti /= ( i + 1 );
    sign = -sign;
    pw  *= xx;
  }
  return result;
}

10000 (GNU C++ v10; -O2)

Meningkatkan akurasi

Metodologi

Akurasi perhitungan fungsi akan ditentukan oleh 2 parameter standar.

Simpangan baku dari dosa sejati (float128) dan rata-rata deviasi ini. Parameter terakhir dapat memberikan informasi penting tentang bagaimana fungsi kita berperilaku. Dia dapat secara sistematis meremehkan atau melebih-lebihkan hasilnya.

Selain parameter ini, kami akan memperkenalkan dua lagi. Bersama dengan fungsi kami, kami memanggil fungsi sin (ganda) yang dibangun ke dalam perpustakaan. Jika hasil dari dua fungsi: milik kita dan yang built-in tidak cocok (sedikit demi sedikit), maka kita tambahkan ke statistik mana dari kedua fungsi tersebut yang lebih jauh dari nilai sebenarnya.

Urutan penjumlahan

Mari kembali ke contoh aslinya. Bagaimana Anda dapat meningkatkan akurasinya "dengan cepat"? Mereka yang telah membaca artikel dengan cermat Apakah mungkin menambahkan nomor N tipe ganda dengan paling akurat? kemungkinan besar akan langsung memberikan jawaban. Anda perlu memutar siklus ke arah yang berlawanan. Untuk menambahkan dari modulo terkecil ke terbesar.

double sin_e3(double x) {
  double result = 0;
  for(int i = 25; i >= 1; i -= 2) {
    result += (((i - 1) % 4 == 0) ? 1 : -1 ) * pow(x, i) / bm::unchecked_factorial<double>(i);
  }
  return result;
}

Hasilnya ditunjukkan di tabel.

Fungsi	Kesalahan rata-rata	STD	Milik kita lebih baik	Libm yang lebih baik
sin_e1	-1.28562e-18	8.25717e-17	0,0588438%	53,5466%
sin_e3	-3.4074e-21	3.39727e-17	0,0423%	10,8049%
sin_e4	8.79046e-18	4.77326e-17	0,0686%	27,6594%
sin_e5	8.78307e-18	3.69995e-17	0,0477062%	13,5105%

Tampaknya menggunakan algoritme penjumlahan cerdas akan menghapus kesalahan hingga hampir 0, tetapi tidak demikian. Tentu saja algoritma tersebut akan memberikan peningkatan akurasi, namun untuk menghilangkan error secara tuntas diperlukan juga algoritma perkalian cerdas. Mereka ada, tetapi sangat mahal: ada banyak operasi yang tidak perlu. Penggunaannya tidak dibenarkan di sini. Namun, kami akan membahasnya nanti dalam konteks yang berbeda.

Hanya ada sedikit yang tersisa. Gabungkan algoritme yang cepat dan akurat. Untuk melakukan ini, mari kembali ke seri Taylor lagi. Mari batasi menjadi 4 anggota misalnya dan buat transformasi berikut.

\sin (x) \approx x (1 + x^{2} (- 1 / 3! + x^{2} (1 / 5! + x^{2} (- 1 / 7! + x^{2} \cdot 1 / 9!))))

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

Anda dapat memperluas tanda kurung dan memeriksa bahwa ekspresi asli diperoleh. Tampilan ini sangat mudah untuk dimasukkan ke dalam satu lingkaran.

double sin_e4(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 1; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x * res;
}

Bekerja cepat, tetapi akurasinya kalah dibandingkan e3. Sekali lagi, pembulatan adalah masalah. Mari kita lihat langkah terakhir dari loop dan ubah sedikit ekspresi aslinya.

\sin (x) \approx x + x \cdot x^{2} (- 1 / 3! + \dots))

$\sin(x)\approx x+x \cdot x^2(-1/3!+ \cdots))$

Dan kode yang sesuai.

double sin_e5(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 3; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x + x * xx * res;
}

Akurasinya menjadi dua kali lipat dibandingkan dengan libm. Jika Anda bisa menebak mengapa keakuratannya meningkat, tulis di komentar. Selain itu, ada satu hal lagi yang tidak menyenangkan tentang sin_e4, yang hilang dari sin_e5, terkait dengan akurasi. Coba tebak apa masalahnya. Di bagian selanjutnya, saya pasti akan membicarakannya secara detail.

Jika Anda menyukai artikel ini, maka selanjutnya saya akan memberi tahu Anda bagaimana GNU libc menghitung sinus dengan ULP maksimum 0,548.

Penghitungan yang akurat dan cepat untuk bilangan floating point menggunakan contoh fungsi sinus. Pendahuluan dan bagian 1

Dekomposisi sinus dalam deret Taylor.

Solusi paling sederhana

Meningkatkan akurasi

Metodologi

Urutan penjumlahan

More articles: