float sinx(float x)
{
static const float a[] = {-.1666666664,.0083333315,-.0001984090,.0000027526,-.0000000239};
float xsq = x*x;
float temp = x*(1 + a[0]*xsq + a[1]*xsq*xsq + a[2]* xsq*xsq*xsq+a[3]*xsq*xsq*xsq*xsq+ a[4]*xsq*xsq*xsq*xsq*xsq);
return temp;
}
这些常数是如何计算的？如何计算cos和tan使用这种方法？我可以扩展它以获得更高的精度吗？我想我需要添加更多常量？
上面描述的“快速”正弦误差与等次泰勒多项式的关系图。
在撰写本文时，几乎所有答案都涉及函数 sin 的泰勒展开，但如果函数的作者是认真的，他不会使用泰勒系数。泰勒系数往往会产生一个多项式近似，它在接近零时比必要的要好，并且在远离零时越来越差。目标通常是在 -π/2…π/2 等范围内获得一致良好的近似值。对于多项式近似，这可以通过应用Remez 算法来获得。一个脚踏实地的解释是这篇文章。
通过该方法获得的多项式系数接近泰勒系数，因为两个多项式都试图逼近相同的函数，但是对于相同数量的操作，多项式可能更精确，或者对于相同（均匀）质量涉及更少的操作的近似值。
我不能仅通过查看它们来判断您问题中的系数是完全泰勒系数还是由 Remez 算法获得的略有不同的系数，但即使不是，它也可能是应该使用的。
最后，无论是谁写的都(1 + a[0]*xsq + a[1]*xsq*xsq + a[2]* xsq*xsq*xsq+a[3]*xsq*xsq*xsq*xsq+ a[4]*xsq*xsq*xsq*xsq*xsq)需要阅读更好的多项式评估方案，例如霍纳的：
1 + xsq*(a[0]+ xsq*(a[1] + xsq*(a[2] + xsq*(a[3] + xsq*a[4]))))使用 N 次乘法而不是 N 2 /2。
它们是-1/6, 1/120, -1/5040.. 等等。
或者更确切地说：-1/3!, 1/5!, -1/7!, 1/9!... 等
在这里查看 sin x 的泰勒级数：
它的正下方有 cos x：
对于 cos x，如上图所示，常数为-1/2!, 1/4!, -1/6!, 1/8!...
tan x 略有不同：
所以要为 cosx 调整它：
float cosx(float x)
{
static const float a[] = {-.5, .0416666667,-.0013888889,.0000248016,-.0000002756};
float xsq = x*x;
float temp = (1 + a[0]*xsq + a[1]*xsq*xsq + a[2]* xsq*xsq*xsq+a[3]*xsq*xsq*xsq*xsq+ a[4]*xsq*xsq*xsq*xsq*xsq);
return temp;
}
这些系数与《数学函数手册》编辑中给出的系数相同。Abramowitz 和 Stegan (1964)，第 76 页，归功于 Carlson 和 Goldstein，Rational approximations of functions，洛斯阿拉莫斯科学实验室 (1955)。
第一个可以在http://www.jonsson.eu/resources/hmf/pdfwrite_600dpi/hmf_600dpi_page_76.pdf中找到。
第二个在http://www.osti.gov/bridge/servlets/purl/4374577-0deJO9/4374577.pdf。第 37 页给出：
关于您的第三个问题，“我可以扩展它以获得更高的精度吗？”，http ://lol.zoy.org/wiki/doc/maths/remez 有一个 Remez 算法的可下载 C++ 实现；它提供（未经我检查）的 6 阶多项式的系数sin：
error: 3.9e-14
9.99999999999624e-1
-1.66666666660981e-1
8.33333330841468e-3
-1.98412650240363e-4
2.75568408741356e-6
-2.50266363478673e-8
1.53659375573646e-10
或者，当然，您需要从 float 更改为 double 才能实现任何改进。这也可以回答您关于cosand的第二个问题tan。
另外，我在评论中看到最后需要一个定点答案。大约 26 年前，我在 8031 汇编器中实现了 32 位定点版本；我会尝试挖掘它，看看它是否有任何有用的东西。
更新：如果您坚持使用 32 位双精度，那么我能看到的将精度提高“一两位”的唯一方法是忘记浮点并使用定点。令人惊讶的是，谷歌似乎没有发现任何东西。以下代码提供了概念验证，在标准 Linux 机器上运行：
#include <stdio.h>
#include <math.h>
#include <stdint.h>
// multiply two 32-bit fixed-point fractions (no rounding)
#define MUL32(a, b) ((uint64_t)(a) * (b) >> 32)
// sin32: Fixed-point sin calculation for first octant, coefficients from
// Handbook for Computing Elementary Functions, by Lyusternik et al, p. 89.
// input: 0 to 0xFFFFFFFF, giving fraction of octant 0 to PI/8, relative to 2**32
// output: 0 to 0.7071, relative to 2**32
static uint32_t sin32(uint32_t x) { // x in 1st octant, = radians/PI*8*2**32
uint32_t y, x2 = MUL32(x, x); // x2 = x * x
y = 0x000259EB; // a7 = 0.000 035 877 1
y = 0x00A32D1E - MUL32(x2, y); // a5 = 0.002 489 871 8
y = 0x14ABBA77 - MUL32(x2, y); // a3 = 0.080 745 367 2
y = 0xC90FDA73u - MUL32(x2, y); // a1 = 0.785 398 152 4
return MUL32(x, y);
}
int main(void) {
int i;
for (i = 0; i < 45; i += 2) { // 0 to 44 degrees
const double two32 = 1LL << 32;
const double radians = i * M_PI / 180;
const uint32_t octant = i / 45. * two32; // fraction of 1st octant
printf("%2d %+.10f %+.10f %+.10f %+.0f\n", i,
sin(radians) - sin32(octant) / two32,
sin(radians) - sinf(radians),
sin(radians) - (float)sin(radians),
sin(radians) * two32 - sin32(octant));
}
return 0;
}
这些系数来自于Lyusternik等人的《计算初等函数手册》 ，第 1 页。89，这里：
我选择这个特定函数的唯一原因是它比你的原始系列少一个术语。
结果是：
0 +0.0000000000 +0.0000000000 +0.0000000000 +0
2 +0.0000000007 +0.0000000003 +0.0000000012 +3
4 +0.0000000010 +0.0000000005 +0.0000000031 +4
6 +0.0000000012 -0.0000000029 -0.0000000011 +5
8 +0.0000000014 +0.0000000011 -0.0000000044 +6
10 +0.0000000014 +0.0000000050 -0.0000000009 +6
12 +0.0000000011 -0.0000000057 +0.0000000057 +5
14 +0.0000000006 -0.0000000018 -0.0000000061 +3
16 -0.0000000000 +0.0000000021 -0.0000000026 -0
18 -0.0000000005 -0.0000000083 -0.0000000082 -2
20 -0.0000000009 +0.0000000095 -0.0000000107 -4
22 -0.0000000010 -0.0000000007 +0.0000000139 -4
24 -0.0000000009 -0.0000000106 +0.0000000010 -4
26 -0.0000000005 +0.0000000065 -0.0000000049 -2
28 -0.0000000001 -0.0000000032 -0.0000000110 -0
30 +0.0000000005 -0.0000000126 -0.0000000000 +2
32 +0.0000000010 +0.0000000037 -0.0000000025 +4
34 +0.0000000015 +0.0000000193 +0.0000000076 +7
36 +0.0000000013 -0.0000000141 +0.0000000083 +6
38 +0.0000000007 +0.0000000011 -0.0000000266 +3
40 -0.0000000005 +0.0000000156 -0.0000000256 -2
42 -0.0000000009 -0.0000000152 -0.0000000170 -4
44 -0.0000000005 -0.0000000011 -0.0000000282 -2
因此我们看到这个定点计算比or精确十倍，并且对 29 位是正确的。使用舍入而不是截断只取得了边际改进。sinf()(float)sin()MUL32()
该函数正在计算sin使用其泰勒展开式的值：
这些常数是各种-1/3！，1/5！依此类推（参见例如这里的泰勒级数的其他函数）。
现在，如果您指定了系列的每个项，则 sin(x) 的泰勒展开对于每个x都是精确的，但是 AFAIK 有更快、更精确的方法来确定软件中的三角函数。
此外，许多处理器提供直接在处理器中实现的此类功能（例如，在 x86 上，它们有现成的操作码），因此通常无需为这些东西烦恼。
cos(x) = sqrt(1 - sin^2(x))
tan(x) = sin(x)/cos(x)
Sin(x) = x -x^3/3! + x^5/5! + (-1)^k*x^(2k+1)/(2k+1)! , k = 1, 2, ...
这是无限函数