C 中正弦函数的实现不起作用
Implementation of sine function in C not working
我试图在 C 中实现正弦函数，但我得到了奇怪的结果。以下是我用来计算正弦的三个函数：
#define PI 3.14159265358979323846
#define DEPTH 16
double sine(long double);
long double pow(long double, unsigned int);
unsigned int fact(unsigned int);
double sine(long double x) {
long double i_x = x *= PI/180;
int n = 3, d = 0, sign = -1;
// fails past 67 degrees
for (; d < DEPTH; n += 2, d++, sign *= -1) {
x += pow(i_x, n) / fact(n) * sign;
}
return x;
}
long double pow(long double base, unsigned int exp) {
double answer = 1;
while (exp) {
answer *= base;
exp--;
}
return answer;
}
unsigned int fact(unsigned int n) {
unsigned int answer = 1;
while (n > 1) {
answer *= n--;
}
return answer;
}
为了测试它，我一直在针对内置正弦函数对其进行测试，如下所示：
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
main() {
for (int i = 0; i <= 180; i++) {
printf("sin(%i) = %lf, %lf\n", i, sine(i), sin(i*3.14159265358979323846/180));
}
exit(EXIT_SUCCESS);
}
向上 67 度，计算与内置函数相同。但是，随着它增加到超过 67，它通常会偏离实际值越来越远。
这是一个示例输出：
>> sin(100) = 0.987711, 0.984808
>> sin(101) = 0.986885, 0.981627
>> sin(102) = 0.987056, 0.978148
>> sin(103) = 0.988830, 0.974370
>> sin(104) = 0.993060, 0.970296
>> sin(105) = 1.000948, 0.965926
>> sin(106) = 1.014169, 0.961262
>> sin(107) = 1.035052, 0.956305
>> sin(108) = 1.066807, 0.951057
>> sin(109) = 1.113846, 0.945519
>> sin(110) = 1.182194, 0.939693
>> sin(111) = 1.280047, 0.933580
>> sin(112) = 1.418502, 0.927184
>> sin(113) = 1.612527, 0.920505
>> sin(114) = 1.882224, 0.913545
>> sin(115) = 2.254492, 0.906308
>> sin(116) = 2.765192, 0.898794
>> sin(117) = 3.461969, 0.891007
...
>> sin(180) = 8431648.192239, 0.000000
有人知道为什么会这样吗？我在 Windows 7 上使用 Visual Studio 2017，如果它提供了任何有用的信息。
你的多项式级数计算方式是numerically unstable. Try horner's method，比幂计算更稳定。
每次你的 for loop 进步，n 增加 2，因此对于 DEPTH = 16，在循环结束时你必须计算数字的阶乘为大到 30 而你使用的 unsigned int 只能存储大到 2^32 = 4294967296 ~= 12! 的值，这会导致你的阶乘函数溢出，进而给你错误的阶乘。
即使你使用 long double 并且我已经在我的评论中声明 MSCRT 中的 long double 映射到 double (Reference) 你仍然会可能在更大的角度看到一些异常，因为虽然 double 可以存储与 1.8E+308 一样大的值，但它在 2^53 = 9007199254740992 ~= 18! 处失去了粒度（即 2^53 + 1 存储为 double 等于 2^53)。因此，一旦你的角度上升，这种行为的影响就会变得越来越大，以至于你在 printf().
中使用的小数点后 6 位精度是显而易见的。
虽然您的方向是正确的，但您应该使用像 GMP 或 libcrypto 这样的 bignum 库。他们可以在不损失精度的情况下执行这些计算。
顺便说一句，由于您是在 Windows 7 上开发的，这意味着您使用的是 x86 或 x86-64。在这些平台上，x87 能够执行 80 位的扩展精度（根据 754 标准）操作，但我不知道编译器内部函数可以在不借助汇编代码的情况下为您提供这种能力。
我还想将您的注意力引向 范围缩小 技术。虽然我仍然建议使用 bignum 库，但如果你在 0 和 90 度之间表现良好（0 和 45，如果我要更严格的话），你可以计算所有其他角度的 sine() 仅通过简单的 三角恒等式.
更新：
实际上我要纠正自己在阶乘计算中使用 double 的问题。在编写了一个简单的程序后，我验证了当我使用 double 存储阶乘时，即使我高于 18，它们也是正确的。经过一番思考后，我意识到在 factorials 的情况下，double 粒度的情况稍微复杂一些。我举个例子说明一下：
19! = 19 * 18 * ... * 2 * 1
在这个数字 18, 16, 14, ... , 2 中都是 2 的倍数并且由于乘以 2 相当于二进制表示中的左移，所以 [= 中的所有低位40=] 已经是 0，因此当 double 对大于 2^53 的 整数 进行舍入时，这些阶乘不受影响。您可以通过计算 2 的数量 16 来计算 19! 的二进制表示中最低有效零的数量。 （对于 20!，它是 18）
我要去 1.8e+308 检查所有阶乘是否不受影响。我会告诉你最新的结果。
更新 2：
如果我们使用 double 来保存阶乘，它们会受到从 23! 开始的四舍五入的影响。它可以很容易地显示出来，因为 2^74 < 23! < 2^75 这意味着至少需要 75 位精度来表示它，但是由于 23! 有 19 个最低有效位，其值为 0, 所以它需要 75 - 19 = 56 比 double.
提供的 53 位大
对于22!来说是51位（可以自己计算）
你的问题在这里：
for (; d < DEPTH; n += 2, d++, sign *= -1) {
x += pow(i_x, n) / fact(n) * sign;
}
您错误地使用了 d < DEPTH，而它应该是 n < DEPTH，d 与您在循环中的计算无关。以下应该可以工作——尽管我还没有编译测试。
for (; n < DEPTH; n += 2, sign *= -1) {
x += pow(i_x, n) / fact(n) * sign;
}
注意： 12 的 DEPTH（例如，带项 1, 3, 5, ... 11 的泰勒级数展开）足以满足 3e-10错误 -- 60-degrees 处的 3 个十亿分之一。 （虽然误差随着 0-360 之间角度的增加而增加，20 的 DEPTH 将在整个范围内保持误差小于 1.0e-8。）
启用编译器警告会捕获sine.
中未使用的d
这里是一个代码变化的例子（注意：Gnu 为 PI 提供了一个常量 M_PI）：
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#define DEPTH 16
/* n factorial */
uint64_t nfact (int n)
{
if (n <= 0) return 1;
uint64_t s = n;
while (--n)
s *= n;
return s;
}
/* y ^ x */
double powerd (const double y, const int x)
{
if (!x) return 1;
double r = y;
for (int i = 1; i < x; i++)
r *= y;
return r;
}
double sine (double deg)
{
double rad = deg * M_PI / 180.0,
x = rad;
int sign = -1;
for (int n = 3; n < DEPTH; n += 2, sign *= -1)
x += sign * powerd (rad, n) / nfact (n);
return x;
}
int main (void) {
printf (" deg sin sine\n\n");
for (int i = 0; i < 180; i++)
printf ("%3d %11.8f %11.8f\n", i, sin (i * M_PI / 180.0), sine (i));
return 0;
}
示例Use/Output
$ ./bin/sine
deg sin sine
0 0.00000000 0.00000000
1 0.01745241 0.01745241
2 0.03489950 0.03489950
3 0.05233596 0.05233596
4 0.06975647 0.06975647
5 0.08715574 0.08715574
6 0.10452846 0.10452846
7 0.12186934 0.12186934
8 0.13917310 0.13917310
9 0.15643447 0.15643447
10 0.17364818 0.17364818
11 0.19080900 0.19080900
12 0.20791169 0.20791169
13 0.22495105 0.22495105
14 0.24192190 0.24192190
15 0.25881905 0.25881905
16 0.27563736 0.27563736
17 0.29237170 0.29237170
18 0.30901699 0.30901699
19 0.32556815 0.32556815
20 0.34202014 0.34202014
21 0.35836795 0.35836795
22 0.37460659 0.37460659
23 0.39073113 0.39073113
24 0.40673664 0.40673664
25 0.42261826 0.42261826
26 0.43837115 0.43837115
27 0.45399050 0.45399050
28 0.46947156 0.46947156
29 0.48480962 0.48480962
30 0.50000000 0.50000000
31 0.51503807 0.51503807
32 0.52991926 0.52991926
33 0.54463904 0.54463904
34 0.55919290 0.55919290
35 0.57357644 0.57357644
36 0.58778525 0.58778525
37 0.60181502 0.60181502
38 0.61566148 0.61566148
39 0.62932039 0.62932039
40 0.64278761 0.64278761
41 0.65605903 0.65605903
42 0.66913061 0.66913061
43 0.68199836 0.68199836
44 0.69465837 0.69465837
45 0.70710678 0.70710678
46 0.71933980 0.71933980
47 0.73135370 0.73135370
48 0.74314483 0.74314483
49 0.75470958 0.75470958
50 0.76604444 0.76604444
51 0.77714596 0.77714596
52 0.78801075 0.78801075
53 0.79863551 0.79863551
54 0.80901699 0.80901699
55 0.81915204 0.81915204
56 0.82903757 0.82903757
57 0.83867057 0.83867057
58 0.84804810 0.84804810
59 0.85716730 0.85716730
60 0.86602540 0.86602540
61 0.87461971 0.87461971
62 0.88294759 0.88294759
63 0.89100652 0.89100652
64 0.89879405 0.89879405
65 0.90630779 0.90630779
66 0.91354546 0.91354546
67 0.92050485 0.92050485
68 0.92718385 0.92718385
69 0.93358043 0.93358043
70 0.93969262 0.93969262
71 0.94551858 0.94551858
72 0.95105652 0.95105652
73 0.95630476 0.95630476
74 0.96126170 0.96126170
75 0.96592583 0.96592583
76 0.97029573 0.97029573
77 0.97437006 0.97437006
78 0.97814760 0.97814760
79 0.98162718 0.98162718
80 0.98480775 0.98480775
81 0.98768834 0.98768834
82 0.99026807 0.99026807
83 0.99254615 0.99254615
84 0.99452190 0.99452190
85 0.99619470 0.99619470
86 0.99756405 0.99756405
87 0.99862953 0.99862953
88 0.99939083 0.99939083
89 0.99984770 0.99984770
90 1.00000000 1.00000000
91 0.99984770 0.99984770
92 0.99939083 0.99939083
93 0.99862953 0.99862953
94 0.99756405 0.99756405
95 0.99619470 0.99619470
96 0.99452190 0.99452190
97 0.99254615 0.99254615
98 0.99026807 0.99026807
99 0.98768834 0.98768834
100 0.98480775 0.98480775
101 0.98162718 0.98162718
102 0.97814760 0.97814760
103 0.97437006 0.97437006
104 0.97029573 0.97029573
105 0.96592583 0.96592583
106 0.96126170 0.96126170
107 0.95630476 0.95630476
108 0.95105652 0.95105652
109 0.94551858 0.94551858
110 0.93969262 0.93969262
111 0.93358043 0.93358043
112 0.92718385 0.92718385
113 0.92050485 0.92050485
114 0.91354546 0.91354546
115 0.90630779 0.90630779
116 0.89879405 0.89879405
117 0.89100652 0.89100652
118 0.88294759 0.88294759
119 0.87461971 0.87461971
120 0.86602540 0.86602540
121 0.85716730 0.85716730
122 0.84804810 0.84804810
123 0.83867057 0.83867057
124 0.82903757 0.82903757
125 0.81915204 0.81915204
126 0.80901699 0.80901699
127 0.79863551 0.79863551
128 0.78801075 0.78801075
129 0.77714596 0.77714596
130 0.76604444 0.76604444
131 0.75470958 0.75470958
132 0.74314483 0.74314482
133 0.73135370 0.73135370
134 0.71933980 0.71933980
135 0.70710678 0.70710678
136 0.69465837 0.69465836
137 0.68199836 0.68199835
138 0.66913061 0.66913060
139 0.65605903 0.65605902
140 0.64278761 0.64278760
141 0.62932039 0.62932038
142 0.61566148 0.61566146
143 0.60181502 0.60181501
144 0.58778525 0.58778523
145 0.57357644 0.57357642
146 0.55919290 0.55919288
147 0.54463904 0.54463901
148 0.52991926 0.52991924
149 0.51503807 0.51503804
150 0.50000000 0.49999996
151 0.48480962 0.48480958
152 0.46947156 0.46947152
153 0.45399050 0.45399045
154 0.43837115 0.43837109
155 0.42261826 0.42261820
156 0.40673664 0.40673657
157 0.39073113 0.39073105
158 0.37460659 0.37460651
159 0.35836795 0.35836786
160 0.34202014 0.34202004
161 0.32556815 0.32556804
162 0.30901699 0.30901686
163 0.29237170 0.29237156
164 0.27563736 0.27563720
165 0.25881905 0.25881887
166 0.24192190 0.24192170
167 0.22495105 0.22495084
168 0.20791169 0.20791145
169 0.19080900 0.19080873
170 0.17364818 0.17364788
171 0.15643447 0.15643414
172 0.13917310 0.13917274
173 0.12186934 0.12186895
174 0.10452846 0.10452803
175 0.08715574 0.08715526
176 0.06975647 0.06975595
177 0.05233596 0.05233537
178 0.03489950 0.03489886
179 0.01745241 0.01745170
基于深度的错误检查
为了回应有关计算误差的评论，您调查了与 sin 和 cos 的泰勒级数展开相关的误差，方法是通过改变 DEPTH 并在 0-360（或 0-2PI）、
范围内使用类似于以下的内容设置最大误差 EMAX 1.0e-8
#define DEPTH 20
#define EMAX 1.0e-8
...
/* sine as above */
...
/* cos with taylor series expansion to n = DEPTH */
long double cose (const long double deg)
{
long double rad = deg * M_PI / 180.0,
x = 1.0;
int sign = -1;
for (int n = 2; n < DEPTH; n += 2, sign *= -1)
x += sign * powerd (rad, n) / nfact (n);
return x;
}
int main (void) {
for (int i = 0; i < 180; i++) {
long double sinlibc = sin (i * M_PI / 180.0),
coslibc = cos (i * M_PI / 180.0),
sints = sine (i),
costs = cose (i),
serr = fabs (sinlibc - sints),
cerr = fabs (coslibc - costs);
if (serr > EMAX)
fprintf (stderr, "sine error exceeds limit of %e\n"
"%3d %11.8Lf %11.8Lf %Le\n",
EMAX, i, sinlibc, sints, serr);
if (cerr > EMAX)
fprintf (stderr, "cose error exceeds limit of %e\n"
"%3d %11.8Lf %11.8Lf %Le\n",
EMAX, i, coslibc, costs, cerr);
}
return 0;
}
如果你检查一下，你会发现对于小于 DEPTH 20（每次展开 10 个项）的任何东西，对于更高的角度，误差将超过 1.0e-8。令人惊讶的是，对于低至 12（6 项）的 DEPTH 值，扩展在第一象限上非常准确。
Addemdum - 使用 0-90 和象限
提高了泰勒级数的准确性
在正常的泰勒级数展开中，误差随着角度的增大而增大。而且......因为有些人不能不修补，我想进一步比较 libc sin/cos 和泰勒级数之间的准确性，如果计算限于 0-90 度和 libc sin/cos 的剩余时间段 90-360 由 0-90 结果的象限 (2, 3 & 4) 镜像处理。它非常有效。
例如，仅处理角度 0-90 和 90 - 180、180 - 270 和 270 - 360 之间的包围角以及初始 angle % 360 的结果会产生结果与 libc 数学库函数相当。 libc 和 8 & 10 术语泰勒级数展开之间的最大误差分别为：
来自 libc 的最大错误 sin/cos
与TSLIM 16
sine_ts max err at : 90.00 deg -- 6.023182e-12
cose_ts max err at : 270.00 deg -- 6.513370e-11
与TSLIM 20
sine_ts max err at : 357.00 deg -- 5.342948e-16
cose_ts max err at : 270.00 deg -- 3.557149e-15
（大量角度显示完全没有区别）
sine 和 cose 与泰勒级数的调整版本如下：
double sine (const double deg)
{
double fp = deg - (int64_t)deg, /* save fractional part of deg */
qdeg = (int64_t)deg % 360, /* get equivalent 0-359 deg angle */
rad, sine_deg; /* radians, sine_deg */
int pos_quad = 1, /* positive quadrant flag 1,2 */
sign = -1; /* taylor series term sign */
qdeg += fp; /* add fractional part back to angle */
/* get equivalent 0-90 degree angle, set pos_quad flag */
if (90 < qdeg && qdeg <= 180) /* in 2nd quadrant */
qdeg = 180 - qdeg;
else if (180 < qdeg && qdeg <= 270) { /* in 3rd quadrant */
qdeg = qdeg - 180;
pos_quad = 0;
}
else if (270 < qdeg && qdeg <= 360) { /* in 4th quadrant */
qdeg = 360 - qdeg;
pos_quad = 0;
}
rad = qdeg * M_PI / 180.0; /* convert to radians */
sine_deg = rad; /* save copy for computation */
/* compute Taylor-Series expansion for sine for TSLIM / 2 terms */
for (int n = 3; n < TSLIM; n += 2, sign *= -1) {
double p = rad;
uint64_t f = n;
for (int i = 1; i < n; i++) /* pow */
p *= rad;
for (int i = 1; i < n; i++) /* nfact */
f *= i;
sine_deg += sign * p / f; /* Taylor-series term */
}
return pos_quad ? sine_deg : -sine_deg;
}
和 cos
double cose (const double deg)
{
double fp = deg - (int64_t)deg, /* save fractional part of deg */
qdeg = (int64_t)deg % 360, /* get equivalent 0-359 deg angle */
rad, cose_deg = 1.0; /* radians, cose_deg */
int pos_quad = 1, /* positive quadrant flag 1,4 */
sign = -1; /* taylor series term sign */
qdeg += fp; /* add fractional part back to angle */
/* get equivalent 0-90 degree angle, set pos_quad flag */
if (90 < qdeg && qdeg <= 180) { /* in 2nd quadrant */
qdeg = 180 - qdeg;
pos_quad = 0;
}
else if (180 < qdeg && qdeg <= 270) { /* in 3rd quadrant */
qdeg = qdeg - 180;
pos_quad = 0;
}
else if (270 < qdeg && qdeg <= 360) /* in 4th quadrant */
qdeg = 360 - qdeg;
rad = qdeg * M_PI / 180.0; /* convert to radians */
/* compute Taylor-Series expansion for sine for TSLIM / 2 terms */
for (int n = 2; n < TSLIM; n += 2, sign *= -1) {
double p = rad;
uint64_t f = n;
for (int i = 1; i < n; i++) /* pow */
p *= rad;
for (int i = 1; i < n; i++) /* nfact */
f *= i;
cose_deg += sign * p / f; /* Taylor-series term */
}
return pos_quad ? cose_deg : -cose_deg;
}
发现兔子尾迹...
您的代码中存在多个问题：
您用不同的原型重新定义了标准函数 pow()。当您 link 将程序作为 en 可执行文件时，这可能会导致问题。使用不同的 anme，例如 pow_int.
您应该在 sine 函数之前将 pow_int 和 fact 函数定义为 static。它可能允许在编译时进行更好的优化。
确实fact受限于类型unsigned int的范围，远小于类型long double的精度。超出 12 的阶乘具有不正确的值，导致精度损失。
您实际上可以增量计算项，节省大量计算并避免潜在的精度损失。
没有参数的main()的原型是int main(void)
PI/180 的计算方式为 double，不如 long double 精确。您应该将表达式写为 x = x * PI / 180;
DEPTH以提高精度。至少还有4个术语带来实质性的改进。
您应该缩小范围：利用正弦函数的对称性和周期性，可以用更少的项对 x 模 90 甚至 45 度进行计算。
这是修改后的版本：
#include <stdio.h>
#include <math.h>
#define PI_L 3.14159265358979323846264338327950288L
#define PI 3.14159265358979323846264338327950288
#define DEPTH 24
double sine(long double x) {
long double res, term, x2, t1;
int phase;
x = remquol(x, 90, &phase);
if (phase & 1)
x = 90 - x;
x = x * PI_L / 180; // convert x to radians
x2 = x * x; // pre-compute x^2
// compute the sine series: x - x^3/3! + x^5/5! ...
res = term = x; // the first term is x
for (int n = 1; n < DEPTH; n += 4) {
// to reduce precision loss, compute 2 terms for each iteration
t1 = term * x2 / ((n + 1) * (n + 2));
term = t1 * x2 / ((n + 3) * (n + 4));
// update the result with the difference of the terms
res += term - t1;
}
if (phase & 2)
res = -res;
return (double)res;
}
int main(void) {
printf("deg sin sine delta\n\n");
for (int i = 0; i <= 360; i += 10) {
double s1 = sin(i * PI / 180);
double s2 = sine(i);
printf("%3i %20.17f %20.17f %g\n", i, s1, s2, s2 - s1);
}
return 0;
}
输出为：
deg sin sine delta
0 0.00000000000000000 0.00000000000000000 0
10 0.17364817766693033 0.17364817766693036 2.77556e-17
20 0.34202014332566871 0.34202014332566871 0
30 0.49999999999999994 0.50000000000000000 5.55112e-17
40 0.64278760968653925 0.64278760968653936 1.11022e-16
50 0.76604444311897801 0.76604444311897801 0
60 0.86602540378443860 0.86602540378443860 0
70 0.93969262078590832 0.93969262078590843 1.11022e-16
80 0.98480775301220802 0.98480775301220802 0
90 1.00000000000000000 1.00000000000000000 0
100 0.98480775301220802 0.98480775301220802 0
110 0.93969262078590843 0.93969262078590843 0
120 0.86602540378443882 0.86602540378443860 -2.22045e-16
130 0.76604444311897812 0.76604444311897801 -1.11022e-16
140 0.64278760968653947 0.64278760968653936 -1.11022e-16
150 0.49999999999999994 0.50000000000000000 5.55112e-17
160 0.34202014332566888 0.34202014332566871 -1.66533e-16
170 0.17364817766693025 0.17364817766693036 1.11022e-16
180 0.00000000000000012 -0.00000000000000000 -1.22465e-16
190 -0.17364817766693047 -0.17364817766693036 1.11022e-16
200 -0.34202014332566866 -0.34202014332566871 -5.55112e-17
210 -0.50000000000000011 -0.50000000000000000 1.11022e-16
220 -0.64278760968653925 -0.64278760968653936 -1.11022e-16
230 -0.76604444311897790 -0.76604444311897801 -1.11022e-16
240 -0.86602540378443837 -0.86602540378443860 -2.22045e-16
250 -0.93969262078590821 -0.93969262078590843 -2.22045e-16
260 -0.98480775301220802 -0.98480775301220802 0
270 -1.00000000000000000 -1.00000000000000000 0
280 -0.98480775301220813 -0.98480775301220802 1.11022e-16
290 -0.93969262078590854 -0.93969262078590843 1.11022e-16
300 -0.86602540378443860 -0.86602540378443860 0
310 -0.76604444311897812 -0.76604444311897801 1.11022e-16
320 -0.64278760968653958 -0.64278760968653936 2.22045e-16
330 -0.50000000000000044 -0.50000000000000000 4.44089e-16
340 -0.34202014332566855 -0.34202014332566871 -1.66533e-16
350 -0.17364817766693127 -0.17364817766693036 9.15934e-16
360 -0.00000000000000024 0.00000000000000000 2.44929e-16
如上所示，sine() 函数似乎比我系统上的标准 sin 函数更精确：sin(180 * M_PI / 128) 应该是 0。同样，sin(150 * M_PI / 128) 应该是 0.5。
将 main 中的角度范围更改为 -90 到 90 仍会覆盖整个正弦范围。但是由于泰勒级数是从零开始的，因此 DEPTH 值可以减少到 7。如前所述，使事实函数为 64 位无符号
将解决 67 度问题。
我试图在 C 中实现正弦函数，但我得到了奇怪的结果。以下是我用来计算正弦的三个函数：
#define PI 3.14159265358979323846
#define DEPTH 16
double sine(long double);
long double pow(long double, unsigned int);
unsigned int fact(unsigned int);
double sine(long double x) {
long double i_x = x *= PI/180;
int n = 3, d = 0, sign = -1;
// fails past 67 degrees
for (; d < DEPTH; n += 2, d++, sign *= -1) {
x += pow(i_x, n) / fact(n) * sign;
}
return x;
}
long double pow(long double base, unsigned int exp) {
double answer = 1;
while (exp) {
answer *= base;
exp--;
}
return answer;
}
unsigned int fact(unsigned int n) {
unsigned int answer = 1;
while (n > 1) {
answer *= n--;
}
return answer;
}
为了测试它，我一直在针对内置正弦函数对其进行测试，如下所示：
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
main() {
for (int i = 0; i <= 180; i++) {
printf("sin(%i) = %lf, %lf\n", i, sine(i), sin(i*3.14159265358979323846/180));
}
exit(EXIT_SUCCESS);
}
向上 67 度，计算与内置函数相同。但是，随着它增加到超过 67，它通常会偏离实际值越来越远。
这是一个示例输出：
>> sin(100) = 0.987711, 0.984808
>> sin(101) = 0.986885, 0.981627
>> sin(102) = 0.987056, 0.978148
>> sin(103) = 0.988830, 0.974370
>> sin(104) = 0.993060, 0.970296
>> sin(105) = 1.000948, 0.965926
>> sin(106) = 1.014169, 0.961262
>> sin(107) = 1.035052, 0.956305
>> sin(108) = 1.066807, 0.951057
>> sin(109) = 1.113846, 0.945519
>> sin(110) = 1.182194, 0.939693
>> sin(111) = 1.280047, 0.933580
>> sin(112) = 1.418502, 0.927184
>> sin(113) = 1.612527, 0.920505
>> sin(114) = 1.882224, 0.913545
>> sin(115) = 2.254492, 0.906308
>> sin(116) = 2.765192, 0.898794
>> sin(117) = 3.461969, 0.891007
...
>> sin(180) = 8431648.192239, 0.000000
有人知道为什么会这样吗？我在 Windows 7 上使用 Visual Studio 2017，如果它提供了任何有用的信息。
你的多项式级数计算方式是numerically unstable. Try horner's method，比幂计算更稳定。
每次你的 for loop 进步，n 增加 2，因此对于 DEPTH = 16，在循环结束时你必须计算数字的阶乘为大到 30 而你使用的 unsigned int 只能存储大到 2^32 = 4294967296 ~= 12! 的值，这会导致你的阶乘函数溢出，进而给你错误的阶乘。
即使你使用 long double 并且我已经在我的评论中声明 MSCRT 中的 long double 映射到 double (Reference) 你仍然会可能在更大的角度看到一些异常，因为虽然 double 可以存储与 1.8E+308 一样大的值，但它在 2^53 = 9007199254740992 ~= 18! 处失去了粒度（即 2^53 + 1 存储为 double 等于 2^53)。因此，一旦你的角度上升，这种行为的影响就会变得越来越大，以至于你在 printf().
虽然您的方向是正确的，但您应该使用像 GMP 或 libcrypto 这样的 bignum 库。他们可以在不损失精度的情况下执行这些计算。
顺便说一句，由于您是在 Windows 7 上开发的，这意味着您使用的是 x86 或 x86-64。在这些平台上，x87 能够执行 80 位的扩展精度（根据 754 标准）操作，但我不知道编译器内部函数可以在不借助汇编代码的情况下为您提供这种能力。
我还想将您的注意力引向 范围缩小 技术。虽然我仍然建议使用 bignum 库，但如果你在 0 和 90 度之间表现良好（0 和 45，如果我要更严格的话），你可以计算所有其他角度的 sine() 仅通过简单的 三角恒等式.
更新：
实际上我要纠正自己在阶乘计算中使用 double 的问题。在编写了一个简单的程序后，我验证了当我使用 double 存储阶乘时，即使我高于 18，它们也是正确的。经过一番思考后，我意识到在 factorials 的情况下，double 粒度的情况稍微复杂一些。我举个例子说明一下：
19! = 19 * 18 * ... * 2 * 1
在这个数字 18, 16, 14, ... , 2 中都是 2 的倍数并且由于乘以 2 相当于二进制表示中的左移，所以 [= 中的所有低位40=] 已经是 0，因此当 double 对大于 2^53 的 整数 进行舍入时，这些阶乘不受影响。您可以通过计算 2 的数量 16 来计算 19! 的二进制表示中最低有效零的数量。 （对于 20!，它是 18）
我要去 1.8e+308 检查所有阶乘是否不受影响。我会告诉你最新的结果。
更新 2：
如果我们使用 double 来保存阶乘，它们会受到从 23! 开始的四舍五入的影响。它可以很容易地显示出来，因为 2^74 < 23! < 2^75 这意味着至少需要 75 位精度来表示它，但是由于 23! 有 19 个最低有效位，其值为 0, 所以它需要 75 - 19 = 56 比 double.
53 位大
对于22!来说是51位（可以自己计算）
你的问题在这里：
for (; d < DEPTH; n += 2, d++, sign *= -1) {
x += pow(i_x, n) / fact(n) * sign;
}
您错误地使用了 d < DEPTH，而它应该是 n < DEPTH，d 与您在循环中的计算无关。以下应该可以工作——尽管我还没有编译测试。
for (; n < DEPTH; n += 2, sign *= -1) {
x += pow(i_x, n) / fact(n) * sign;
}
注意： 12 的 DEPTH（例如，带项 1, 3, 5, ... 11 的泰勒级数展开）足以满足 3e-10错误 -- 60-degrees 处的 3 个十亿分之一。 （虽然误差随着 0-360 之间角度的增加而增加，20 的 DEPTH 将在整个范围内保持误差小于 1.0e-8。）
启用编译器警告会捕获sine.
d
这里是一个代码变化的例子（注意：Gnu 为 PI 提供了一个常量 M_PI）：
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#define DEPTH 16
/* n factorial */
uint64_t nfact (int n)
{
if (n <= 0) return 1;
uint64_t s = n;
while (--n)
s *= n;
return s;
}
/* y ^ x */
double powerd (const double y, const int x)
{
if (!x) return 1;
double r = y;
for (int i = 1; i < x; i++)
r *= y;
return r;
}
double sine (double deg)
{
double rad = deg * M_PI / 180.0,
x = rad;
int sign = -1;
for (int n = 3; n < DEPTH; n += 2, sign *= -1)
x += sign * powerd (rad, n) / nfact (n);
return x;
}
int main (void) {
printf (" deg sin sine\n\n");
for (int i = 0; i < 180; i++)
printf ("%3d %11.8f %11.8f\n", i, sin (i * M_PI / 180.0), sine (i));
return 0;
}
示例Use/Output
$ ./bin/sine
deg sin sine
0 0.00000000 0.00000000
1 0.01745241 0.01745241
2 0.03489950 0.03489950
3 0.05233596 0.05233596
4 0.06975647 0.06975647
5 0.08715574 0.08715574
6 0.10452846 0.10452846
7 0.12186934 0.12186934
8 0.13917310 0.13917310
9 0.15643447 0.15643447
10 0.17364818 0.17364818
11 0.19080900 0.19080900
12 0.20791169 0.20791169
13 0.22495105 0.22495105
14 0.24192190 0.24192190
15 0.25881905 0.25881905
16 0.27563736 0.27563736
17 0.29237170 0.29237170
18 0.30901699 0.30901699
19 0.32556815 0.32556815
20 0.34202014 0.34202014
21 0.35836795 0.35836795
22 0.37460659 0.37460659
23 0.39073113 0.39073113
24 0.40673664 0.40673664
25 0.42261826 0.42261826
26 0.43837115 0.43837115
27 0.45399050 0.45399050
28 0.46947156 0.46947156
29 0.48480962 0.48480962
30 0.50000000 0.50000000
31 0.51503807 0.51503807
32 0.52991926 0.52991926
33 0.54463904 0.54463904
34 0.55919290 0.55919290
35 0.57357644 0.57357644
36 0.58778525 0.58778525
37 0.60181502 0.60181502
38 0.61566148 0.61566148
39 0.62932039 0.62932039
40 0.64278761 0.64278761
41 0.65605903 0.65605903
42 0.66913061 0.66913061
43 0.68199836 0.68199836
44 0.69465837 0.69465837
45 0.70710678 0.70710678
46 0.71933980 0.71933980
47 0.73135370 0.73135370
48 0.74314483 0.74314483
49 0.75470958 0.75470958
50 0.76604444 0.76604444
51 0.77714596 0.77714596
52 0.78801075 0.78801075
53 0.79863551 0.79863551
54 0.80901699 0.80901699
55 0.81915204 0.81915204
56 0.82903757 0.82903757
57 0.83867057 0.83867057
58 0.84804810 0.84804810
59 0.85716730 0.85716730
60 0.86602540 0.86602540
61 0.87461971 0.87461971
62 0.88294759 0.88294759
63 0.89100652 0.89100652
64 0.89879405 0.89879405
65 0.90630779 0.90630779
66 0.91354546 0.91354546
67 0.92050485 0.92050485
68 0.92718385 0.92718385
69 0.93358043 0.93358043
70 0.93969262 0.93969262
71 0.94551858 0.94551858
72 0.95105652 0.95105652
73 0.95630476 0.95630476
74 0.96126170 0.96126170
75 0.96592583 0.96592583
76 0.97029573 0.97029573
77 0.97437006 0.97437006
78 0.97814760 0.97814760
79 0.98162718 0.98162718
80 0.98480775 0.98480775
81 0.98768834 0.98768834
82 0.99026807 0.99026807
83 0.99254615 0.99254615
84 0.99452190 0.99452190
85 0.99619470 0.99619470
86 0.99756405 0.99756405
87 0.99862953 0.99862953
88 0.99939083 0.99939083
89 0.99984770 0.99984770
90 1.00000000 1.00000000
91 0.99984770 0.99984770
92 0.99939083 0.99939083
93 0.99862953 0.99862953
94 0.99756405 0.99756405
95 0.99619470 0.99619470
96 0.99452190 0.99452190
97 0.99254615 0.99254615
98 0.99026807 0.99026807
99 0.98768834 0.98768834
100 0.98480775 0.98480775
101 0.98162718 0.98162718
102 0.97814760 0.97814760
103 0.97437006 0.97437006
104 0.97029573 0.97029573
105 0.96592583 0.96592583
106 0.96126170 0.96126170
107 0.95630476 0.95630476
108 0.95105652 0.95105652
109 0.94551858 0.94551858
110 0.93969262 0.93969262
111 0.93358043 0.93358043
112 0.92718385 0.92718385
113 0.92050485 0.92050485
114 0.91354546 0.91354546
115 0.90630779 0.90630779
116 0.89879405 0.89879405
117 0.89100652 0.89100652
118 0.88294759 0.88294759
119 0.87461971 0.87461971
120 0.86602540 0.86602540
121 0.85716730 0.85716730
122 0.84804810 0.84804810
123 0.83867057 0.83867057
124 0.82903757 0.82903757
125 0.81915204 0.81915204
126 0.80901699 0.80901699
127 0.79863551 0.79863551
128 0.78801075 0.78801075
129 0.77714596 0.77714596
130 0.76604444 0.76604444
131 0.75470958 0.75470958
132 0.74314483 0.74314482
133 0.73135370 0.73135370
134 0.71933980 0.71933980
135 0.70710678 0.70710678
136 0.69465837 0.69465836
137 0.68199836 0.68199835
138 0.66913061 0.66913060
139 0.65605903 0.65605902
140 0.64278761 0.64278760
141 0.62932039 0.62932038
142 0.61566148 0.61566146
143 0.60181502 0.60181501
144 0.58778525 0.58778523
145 0.57357644 0.57357642
146 0.55919290 0.55919288
147 0.54463904 0.54463901
148 0.52991926 0.52991924
149 0.51503807 0.51503804
150 0.50000000 0.49999996
151 0.48480962 0.48480958
152 0.46947156 0.46947152
153 0.45399050 0.45399045
154 0.43837115 0.43837109
155 0.42261826 0.42261820
156 0.40673664 0.40673657
157 0.39073113 0.39073105
158 0.37460659 0.37460651
159 0.35836795 0.35836786
160 0.34202014 0.34202004
161 0.32556815 0.32556804
162 0.30901699 0.30901686
163 0.29237170 0.29237156
164 0.27563736 0.27563720
165 0.25881905 0.25881887
166 0.24192190 0.24192170
167 0.22495105 0.22495084
168 0.20791169 0.20791145
169 0.19080900 0.19080873
170 0.17364818 0.17364788
171 0.15643447 0.15643414
172 0.13917310 0.13917274
173 0.12186934 0.12186895
174 0.10452846 0.10452803
175 0.08715574 0.08715526
176 0.06975647 0.06975595
177 0.05233596 0.05233537
178 0.03489950 0.03489886
179 0.01745241 0.01745170
基于深度的错误检查
为了回应有关计算误差的评论，您调查了与 sin 和 cos 的泰勒级数展开相关的误差，方法是通过改变 DEPTH 并在 0-360（或 0-2PI）、
EMAX 1.0e-8
#define DEPTH 20
#define EMAX 1.0e-8
...
/* sine as above */
...
/* cos with taylor series expansion to n = DEPTH */
long double cose (const long double deg)
{
long double rad = deg * M_PI / 180.0,
x = 1.0;
int sign = -1;
for (int n = 2; n < DEPTH; n += 2, sign *= -1)
x += sign * powerd (rad, n) / nfact (n);
return x;
}
int main (void) {
for (int i = 0; i < 180; i++) {
long double sinlibc = sin (i * M_PI / 180.0),
coslibc = cos (i * M_PI / 180.0),
sints = sine (i),
costs = cose (i),
serr = fabs (sinlibc - sints),
cerr = fabs (coslibc - costs);
if (serr > EMAX)
fprintf (stderr, "sine error exceeds limit of %e\n"
"%3d %11.8Lf %11.8Lf %Le\n",
EMAX, i, sinlibc, sints, serr);
if (cerr > EMAX)
fprintf (stderr, "cose error exceeds limit of %e\n"
"%3d %11.8Lf %11.8Lf %Le\n",
EMAX, i, coslibc, costs, cerr);
}
return 0;
}
如果你检查一下，你会发现对于小于 DEPTH 20（每次展开 10 个项）的任何东西，对于更高的角度，误差将超过 1.0e-8。令人惊讶的是，对于低至 12（6 项）的 DEPTH 值，扩展在第一象限上非常准确。
Addemdum - 使用 0-90 和象限
在正常的泰勒级数展开中，误差随着角度的增大而增大。而且......因为有些人不能不修补，我想进一步比较 libc sin/cos 和泰勒级数之间的准确性，如果计算限于 0-90 度和 libc sin/cos 的剩余时间段 90-360 由 0-90 结果的象限 (2, 3 & 4) 镜像处理。它非常有效。
例如，仅处理角度 0-90 和 90 - 180、180 - 270 和 270 - 360 之间的包围角以及初始 angle % 360 的结果会产生结果与 libc 数学库函数相当。 libc 和 8 & 10 术语泰勒级数展开之间的最大误差分别为：
来自 libc 的最大错误 sin/cos
与TSLIM 16
sine_ts max err at : 90.00 deg -- 6.023182e-12
cose_ts max err at : 270.00 deg -- 6.513370e-11
与TSLIM 20
sine_ts max err at : 357.00 deg -- 5.342948e-16
cose_ts max err at : 270.00 deg -- 3.557149e-15
（大量角度显示完全没有区别）
sine 和 cose 与泰勒级数的调整版本如下：
double sine (const double deg)
{
double fp = deg - (int64_t)deg, /* save fractional part of deg */
qdeg = (int64_t)deg % 360, /* get equivalent 0-359 deg angle */
rad, sine_deg; /* radians, sine_deg */
int pos_quad = 1, /* positive quadrant flag 1,2 */
sign = -1; /* taylor series term sign */
qdeg += fp; /* add fractional part back to angle */
/* get equivalent 0-90 degree angle, set pos_quad flag */
if (90 < qdeg && qdeg <= 180) /* in 2nd quadrant */
qdeg = 180 - qdeg;
else if (180 < qdeg && qdeg <= 270) { /* in 3rd quadrant */
qdeg = qdeg - 180;
pos_quad = 0;
}
else if (270 < qdeg && qdeg <= 360) { /* in 4th quadrant */
qdeg = 360 - qdeg;
pos_quad = 0;
}
rad = qdeg * M_PI / 180.0; /* convert to radians */
sine_deg = rad; /* save copy for computation */
/* compute Taylor-Series expansion for sine for TSLIM / 2 terms */
for (int n = 3; n < TSLIM; n += 2, sign *= -1) {
double p = rad;
uint64_t f = n;
for (int i = 1; i < n; i++) /* pow */
p *= rad;
for (int i = 1; i < n; i++) /* nfact */
f *= i;
sine_deg += sign * p / f; /* Taylor-series term */
}
return pos_quad ? sine_deg : -sine_deg;
}
和 cos
double cose (const double deg)
{
double fp = deg - (int64_t)deg, /* save fractional part of deg */
qdeg = (int64_t)deg % 360, /* get equivalent 0-359 deg angle */
rad, cose_deg = 1.0; /* radians, cose_deg */
int pos_quad = 1, /* positive quadrant flag 1,4 */
sign = -1; /* taylor series term sign */
qdeg += fp; /* add fractional part back to angle */
/* get equivalent 0-90 degree angle, set pos_quad flag */
if (90 < qdeg && qdeg <= 180) { /* in 2nd quadrant */
qdeg = 180 - qdeg;
pos_quad = 0;
}
else if (180 < qdeg && qdeg <= 270) { /* in 3rd quadrant */
qdeg = qdeg - 180;
pos_quad = 0;
}
else if (270 < qdeg && qdeg <= 360) /* in 4th quadrant */
qdeg = 360 - qdeg;
rad = qdeg * M_PI / 180.0; /* convert to radians */
/* compute Taylor-Series expansion for sine for TSLIM / 2 terms */
for (int n = 2; n < TSLIM; n += 2, sign *= -1) {
double p = rad;
uint64_t f = n;
for (int i = 1; i < n; i++) /* pow */
p *= rad;
for (int i = 1; i < n; i++) /* nfact */
f *= i;
cose_deg += sign * p / f; /* Taylor-series term */
}
return pos_quad ? cose_deg : -cose_deg;
}
发现兔子尾迹...
您的代码中存在多个问题：
您用不同的原型重新定义了标准函数
pow()。当您 link 将程序作为 en 可执行文件时，这可能会导致问题。使用不同的 anme，例如pow_int.您应该在
sine函数之前将pow_int和fact函数定义为static。它可能允许在编译时进行更好的优化。确实
fact受限于类型unsigned int的范围，远小于类型long double的精度。超出12的阶乘具有不正确的值，导致精度损失。您实际上可以增量计算项，节省大量计算并避免潜在的精度损失。
没有参数的
main()的原型是int main(void)PI/180的计算方式为double，不如long double精确。您应该将表达式写为x = x * PI / 180;
应增加DEPTH以提高精度。至少还有4个术语带来实质性的改进。您应该缩小范围：利用正弦函数的对称性和周期性，可以用更少的项对
x模 90 甚至 45 度进行计算。
这是修改后的版本：
#include <stdio.h>
#include <math.h>
#define PI_L 3.14159265358979323846264338327950288L
#define PI 3.14159265358979323846264338327950288
#define DEPTH 24
double sine(long double x) {
long double res, term, x2, t1;
int phase;
x = remquol(x, 90, &phase);
if (phase & 1)
x = 90 - x;
x = x * PI_L / 180; // convert x to radians
x2 = x * x; // pre-compute x^2
// compute the sine series: x - x^3/3! + x^5/5! ...
res = term = x; // the first term is x
for (int n = 1; n < DEPTH; n += 4) {
// to reduce precision loss, compute 2 terms for each iteration
t1 = term * x2 / ((n + 1) * (n + 2));
term = t1 * x2 / ((n + 3) * (n + 4));
// update the result with the difference of the terms
res += term - t1;
}
if (phase & 2)
res = -res;
return (double)res;
}
int main(void) {
printf("deg sin sine delta\n\n");
for (int i = 0; i <= 360; i += 10) {
double s1 = sin(i * PI / 180);
double s2 = sine(i);
printf("%3i %20.17f %20.17f %g\n", i, s1, s2, s2 - s1);
}
return 0;
}
输出为：
deg sin sine delta
0 0.00000000000000000 0.00000000000000000 0
10 0.17364817766693033 0.17364817766693036 2.77556e-17
20 0.34202014332566871 0.34202014332566871 0
30 0.49999999999999994 0.50000000000000000 5.55112e-17
40 0.64278760968653925 0.64278760968653936 1.11022e-16
50 0.76604444311897801 0.76604444311897801 0
60 0.86602540378443860 0.86602540378443860 0
70 0.93969262078590832 0.93969262078590843 1.11022e-16
80 0.98480775301220802 0.98480775301220802 0
90 1.00000000000000000 1.00000000000000000 0
100 0.98480775301220802 0.98480775301220802 0
110 0.93969262078590843 0.93969262078590843 0
120 0.86602540378443882 0.86602540378443860 -2.22045e-16
130 0.76604444311897812 0.76604444311897801 -1.11022e-16
140 0.64278760968653947 0.64278760968653936 -1.11022e-16
150 0.49999999999999994 0.50000000000000000 5.55112e-17
160 0.34202014332566888 0.34202014332566871 -1.66533e-16
170 0.17364817766693025 0.17364817766693036 1.11022e-16
180 0.00000000000000012 -0.00000000000000000 -1.22465e-16
190 -0.17364817766693047 -0.17364817766693036 1.11022e-16
200 -0.34202014332566866 -0.34202014332566871 -5.55112e-17
210 -0.50000000000000011 -0.50000000000000000 1.11022e-16
220 -0.64278760968653925 -0.64278760968653936 -1.11022e-16
230 -0.76604444311897790 -0.76604444311897801 -1.11022e-16
240 -0.86602540378443837 -0.86602540378443860 -2.22045e-16
250 -0.93969262078590821 -0.93969262078590843 -2.22045e-16
260 -0.98480775301220802 -0.98480775301220802 0
270 -1.00000000000000000 -1.00000000000000000 0
280 -0.98480775301220813 -0.98480775301220802 1.11022e-16
290 -0.93969262078590854 -0.93969262078590843 1.11022e-16
300 -0.86602540378443860 -0.86602540378443860 0
310 -0.76604444311897812 -0.76604444311897801 1.11022e-16
320 -0.64278760968653958 -0.64278760968653936 2.22045e-16
330 -0.50000000000000044 -0.50000000000000000 4.44089e-16
340 -0.34202014332566855 -0.34202014332566871 -1.66533e-16
350 -0.17364817766693127 -0.17364817766693036 9.15934e-16
360 -0.00000000000000024 0.00000000000000000 2.44929e-16
如上所示，sine() 函数似乎比我系统上的标准 sin 函数更精确：sin(180 * M_PI / 128) 应该是 0。同样，sin(150 * M_PI / 128) 应该是 0.5。
将 main 中的角度范围更改为 -90 到 90 仍会覆盖整个正弦范围。但是由于泰勒级数是从零开始的，因此 DEPTH 值可以减少到 7。如前所述，使事实函数为 64 位无符号 将解决 67 度问题。