系数与数学函数手册中给出的系数相同。 Abramowitz 和 Stegan (1964),第 76 页,归功于 Carlson 和 Goldstein,Rational approximations of functions,Los Alamos Scientific Laboratory (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++ 实现;它提供(未经我检查)sin 的 6 阶多项式的系数:
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 才能实现任何改进。这也可以回答您关于cos 和tan 的第二个问题。
另外,我在 cmets 中看到最终需要一个定点答案。大约 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;
}
系数来自 Handbook for Computing Elementary Functions,Lyusternik 等人,p。 89、here:
我选择这个特殊函数的唯一原因是它比你原来的系列少了一个词。
结果是:
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
因此,我们看到这个定点计算比sinf() 或(float)sin() 精确大约十倍,并且正确到29 位。在 MUL32() 中使用舍入而不是截断只取得了边际改进。