🤽🏽 🚏 👸🏾 使用正弦函数示例可以快速准确地计算浮点数。引言和第1部分 👩🏻 🦌 ☕️

仔细阅读来自阿特姆·卡拉瓦耶夫加上浮点数。该主题非常有趣，我想继续进行下去，并通过示例演示如何在实践中使用浮点数。让我们以GNU glibc库（libm）作为参考。为了使本文不会太无聊，我们将添加一个具有竞争力的组件：我们不仅将尝试重复，而且还将尝试改进库代码，从而使其更快/更准确。

例如，我选择了三角正弦函数。这是一项广泛的功能，其数学在学校和大学中是众所周知的。同时，随着其实现，将有许多生动的数字“正确”工作示例。我将使用double作为浮点数。

在本系列文章中，计划了很多内容，从数学到机器代码和编译器选项。撰写文章的语言是C ++，但没有“多余”。与C语言相反，即使对于不熟悉C语言的人来说，工作示例也将更具可读性并且占用更少的行。

文章将采用浸入法写作。将讨论子任务，然后将这些子任务组合在一起以解决问题。

将正弦分解为泰勒级数。

正弦函数扩展为无限泰勒级数。

\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$

显然，除了存在无限和的解析公式的情况以外，我们无法计算无限级数。但这不是我们的情况）））假设我们要计算区间中的正弦

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

$[0, \frac{\pi}{2}]$ ... 我们将在第3部分中更详细地讨论工作。

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

$\sin(\frac{\pi}{2})=1$ 估计，根据以下条件找到可以舍弃的第一项

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

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

e

$e$ 它是数字1和大于1的最小数字之间的差。大致来说，这是尾数（wiki）的最后一位。用蛮力解方程很容易。对于

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

$e \approx 2.22\times10^{-16}$ ... 我管理

n = 23

$n=23$ 可以被删除。术语数量的正确选择将在以下部分之一中讨论，因此，今天我们将“安全”并将这些术语应用于

n = 25

$n=25$ 包括的。

最后一个学期比大约少10,000倍

e

$e$ ...

最简单的解决方案

手已经发痒了，我们这样写：

测试程序的全文

#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;
}

如何加快程序速度，我认为很多人马上就想到了。您认为您的更改可以加速程序多少次？优化版本和扰流板下问题的答案。

程序的优化版本。

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)

提高准确性

方法

函数的计算精度将由2个标准参数确定。

与真实正弦的标准偏差（float128）和该偏差的平均值。最后一个参数可以提供有关函数功能的重要信息。她可以系统地低估或高估结果。

除了这些参数外，我们还将介绍另外两个。与我们的函数一起，我们调用内置在库中的sin（double）函数。如果两个函数（我们的函数和内置函数）的结果不匹配（逐位），则我们将两个函数中的哪个函数与真实值相去甚远，添加到统计信息中。

求和顺序

让我们回到原始示例。您如何“迅速”提高其准确性？那些认真阅读文章的人是否可以最准确地添加N个双精度类型的数字？很有可能会立即给出答案。必须以相反的方向转动循环。从最小模到最大模相加。

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;
}

结果显示在表中。

功能	平均误差	性病	更好的我们	更好的libm
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％

使用智能求和算法似乎可以将错误几乎消除为0，但事实并非如此。当然，这些算法将提高准确性，但是要完全消除错误，还需要智能乘法算法。它们存在，但是它们非常昂贵：有很多不必要的操作。在这里没有理由使用它们。但是，稍后我们将在不同的背景下再次讨论它们。

剩下的很少了。结合快速准确的算法。为此，让我们再次回到泰勒系列。让我们将其限制为例如4个成员，并进行以下转换。

\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!))))$

您可以展开括号并检查是否获得了原始表达式。此视图非常容易放入循环中。

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;
}

运行速度很快，但与e3相比，准确性有所下降。同样，舍入是一个问题。让我们看一下循环的最后一步，对原始表达式进行一些转换。

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

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

和相应的代码。

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;
}

与libm相比，准确性提高了一倍。如果您可以猜测为什么精度提高了，请在注释中写下。另外，关于sin_e4还有更多，更不愉快的事情，sin_e5中缺少与准确性有关的事情。尝试猜测是什么问题。在下一部分中，我一定会详细告诉您。

如果您喜欢这篇文章，那么在下一篇文章中，我将告诉您GNU libc如何计算最大ULP为0.548的正弦值。

使用正弦函数示例可以快速准确地计算浮点数。引言和第1部分

将正弦分解为泰勒级数。

最简单的解决方案

提高准确性

方法

求和顺序

More articles: