IEEE 浮点异常 - 为什么？
IEEE floating point exception - why?
我正在制作一个简单的函数来重现 Secant method，但我认为我在精度上遇到了一些问题。
首先，这是函数（以及调用它的主要方法，以及与之一起使用的测试函数）：
double secant_method(double(*f)(double), double a, double b){
double c;
for (int i = 0; i < 10; i++){
c = a - f(a) * (a - b) / (f(a) - f(b));
b = a; a = c;
}
return c;
}
typedef double (*func)(double);
//test function - x^3 + 4x - 10
double test(double x){
return (x*x*x) + (4*x) - 10;
}
int main(){
func f = &test;
double ans;
ans = secant_method(f, 0, 2);
printf("\nRoot found:\t%.*g\n", DECIMAL_DIG * 2, ans);
return 0;
}
注意：函数secant_method()中的for循环只循环了10次。这就是我的问题所在。
按原样打印时，一切正常。它给出了小数点后 16 位的正确输出：
但是，当我向 secant_method() 中的 for 循环添加迭代时，会发生这种情况：
为什么会这样？我是否达到了 C 可以处理的最大表示？
我通读了另一个 post 的 this great answer，但在其中，关于我收到的异常 (-1.#IND)，它只是说我的结果不是号码，或者我在进行某种非法操作。
编辑： 使用 (x*x*x) + (4*x) - 10 + sin(x) 作为我的测试函数给出了正确的答案 - 但只有当我在 i < 9 而不是 [=18= 时循环] (x*x*x) + (4*x) - 10
-1.#IND是微软输出不确定值的方式，具体来说就是NaN.
一种 可能发生这种情况的方式是 0 / 0 但我会检查 所有 操作以查看问题出在哪里谎言：
double secant_method(double(*f)(double), double a, double b){
double c;
printf("DBG =====\n");
for (int i = 0; i < 10; i++){
printf("\nDBG -----\n");
printf("DBG i: %d\n",i);
printf("DBG a: %30f\n",a);
printf("DBG b: %30f\n",b);
printf("DBG c: %30f\n",c);
printf("DBG f(a): %30f\n",f(a));
printf("DBG a-b: %30f\n",a-b);
printf("DBG f(b): %30f\n",f(b));
printf("DBG f(a)-f(b): %30f\n",f(a)-f(b));
printf("DBG f(a)*(a-b): %30f\n",f(a)*(a-b));
printf("DBG f(a)*(a-b)/(f(a)-f(b)): %30f\n",f(a)*(a-b)/(f(a)-f(b)));
c = a - f(a) * (a - b) / (f(a) - f(b));
b = a; a = c;
}
return c;
}
一旦你有了调试输出，然后你就可以弄清楚实际问题是什么，并采取策略来避免它。
当我这样做时，我看到（最后）：
DBG -----
DBG i: 8
DBG a: 1.556773264394211375716281509085
DBG b: 1.556773264393484179635152031551
DBG c: 1.556773264394211375716281509085
DBG f(a): -0.000000000000000987057657830803
DBG a-b: 0.000000000000727196081129477534
DBG f(b): -0.000000000008196943991622962500
DBG f(a)-f(b): 0.000000000008195956933965131697
DBG f(a)*(a-b): -0.000000000000000000000000000718
DBG f(a)*(a-b)/(f(a)-f(b)): -0.000000000000000087577871187781
DBG -----
DBG i: 9
DBG a: 1.556773264394211375716281509085
DBG b: 1.556773264394211375716281509085
DBG c: 1.556773264394211375716281509085
DBG f(a): -0.000000000000000987057657830803
DBG a-b: 0.000000000000000000000000000000
DBG f(b): -0.000000000000000987057657830803
DBG f(a)-f(b): 0.000000000000000000000000000000
DBG f(a)*(a-b): -0.000000000000000000000000000000
DBG f(a)*(a-b)/(f(a)-f(b)): nan
Root found: nan
所以你可以看到，在第十次迭代中，a 和 b 变得相等，因此 f(a) 和 f(b) 也相等。所以你得到了表达式：
something * 0 / 0
如前所述，这会给你 0 / 0 或 NaN。
关于如何修复它，您只需要避免除以零，因为这会给您带来 NaN 或无穷大。因此，您可以改用以下函数：
double secant_method(double(*f)(double), double a, double b){
double c;
for (int i = 0; i < 1000; i++) {
if (f(a) == f(b)) break;
c = a - f(a) * (a - b) / (f(a) - f(b));
b = a; a = c;
}
return c;
}
一千个循环应该足以得到一个像样的答案，如果你要除以零，它会提前退出。
如果您想要 更高 精度，您可以查看 long double 类型或切换到使用任意精度算术库之一，例如 GMP 或 MPIR。
这通常需要更多的工作，但您可以获得一些令人印象深刻的结果。这个程序，建立在 MPIR 上：
#include <stdio.h>
#include <mpir.h>
void secant_method(mpf_t result, void(*f)(mpf_t, mpf_t), mpf_t a, mpf_t b){
mpf_t c, fa, fb, temp1, temp2;
mpf_init (fa);
mpf_init (fb);
mpf_init (temp1);
mpf_init (temp2);
for (int i = 0; i < 1000; i++){
printf("DBG i: %d\n",i);
f (fa, a);
f (fb, b);
if (mpf_cmp (fa, fb) == 0) break;
mpf_set (temp1, a);
mpf_sub (temp1, temp1, b);
mpf_set (temp2, fa);
mpf_sub (temp2, temp2, fb);
mpf_set (result, fa);
mpf_mul (result, result, temp1);
mpf_div (result, result, temp2);
mpf_sub (result, result, a);
mpf_neg (result, result);
mpf_set (b, a);
mpf_set (a, result);
}
}
void test (mpf_t result, mpf_t x){
mpf_t temp;
mpf_set (result, x);
mpf_pow_ui (result, result, 3);
mpf_init_set (temp, x);
mpf_mul_ui (temp, temp, 4);
mpf_add (result, result, temp);
mpf_set_ui (temp, 10);
mpf_sub (result, result, temp);
mpf_clear (temp);
}
int main(){
mpf_t ans, a, b;
mpf_set_default_prec (8000);
mpf_init_set_ui (ans, 0);
mpf_init_set_ui (a, 0);
mpf_init_set_ui (b, 2);
secant_method (ans, &test, a, b);
mpf_out_str (stdout, 10, 0, ans);
return 0;
}
输出精度更高，大约二分半 千 位：
DBG i: 1
:
DBG i: 19
0.155677326439421146326886324730853302634853266143
22856485101283627988036767055520913212330822780959
93349183787687346999781239000417393618333668026011
02048595843228945228507966189601958673920851932189
20626590635658264390975889008832048255537650792123
54916373054888140164770654992918100928227714960414
65208113116379497717707745267800989233875981344305
90022883167106124203999713536673991376957068731244
91919087980169395013246250812213656324598765244218
15974098310512802880727074335472786858740154287363
31949470951650710072488856623955478366217474755111
76368234254761541647442609230138418167182918204711
66713459423756284737546964906061587903876515793884
14091165347411853670752820576131460960421137744435
73729141652832258144582021037373967987171478026002
48487515446248979731517957120705447608265161099693
33098235693813752370774508652788986557620510981156
19907950657355934071535840759135251701581523712307
00051674680667972152582339710574822560693109306285
91240827697915787078746087225027856691436076089912
35551789799825731841345891629028445554314717823386
07885164744100235567602875364878328805811271289098
87558119684442289199181352023304600847178256323082
57317198584882656089836229208443415369358460418542
84083408696290686178971039756668669303212658278679
39542421457300944206839268283788585029652481323614
65995074020560963212330914882733926627309382310653
39023265929195094492468196461296569155421718696631
73798097369621805062145075113127308161572398104766
37356504104570136778437926442139603916930640425421
15655156674699552536588332891562053247342008145504
44336211031437923307615880759201695011419324719812
46482293928341901673056596202744639074280785106031
90197472588293352508389295101867514582271001202777
85575614897203080940643669476500979934666490279524
88486176409290187337498631681392563044899541391612
88438904336237873504970887963071622208868799638373
42186338496601471274609131141920820263780493617795
89714798662834913192777810386631915415021934333441
01797098172897161215116673422762953435902633516501
73788202968876596925999628999004575114529754782488
59959395407324243559011982543407738505315960009874
36510513519775603567237051670918870105777288994910
85524037720122749091827520695838000086150188462000
63190624219373460624686216781527327604063990319908
56547016812115842640285111265677758613385414834511
69237199199725030839166586376374587900611430229333
87296847315023767826706323911923435564643861604120
017381909481e1
而且，如果你把那个数字传回 test() 函数，你会得到一个非常接近于零的数字，大约 -1.15 x 10<sup>- 2408</sup>。所以我想说这是对使用 double.
的一个相当大的改进
而且，就其价值而言，它只需要大约十分之一秒 CPU 的时间，所以至少用任意精度的算法来做到这一点是可行的。
甚至更精度，只需更改MPIR的默认精度设置，目前设置为：
mpf_set_default_prec (8000);
增加到 100,000 会给出超过 30,000 位有效数字的答案，最终 "close-to-zero" 答案约为 -5 x 10<sup>-30103</sup>.
#IND是由零除以零引起的。调试时在循环中插入一个简单的 isnan(c) 检查，你会发现 a 和 b 最终变得相等，这导致 a - b 和 f(a) - f(b) 为零。
我认为选角是个问题。编译器不会将您所有的常数解释为 double，而是将其解释为常规整数。添加一个“。” （或“.0”）在每个常量的末尾使它们加倍。因此 test() 应更改为：
return (x*x*x) + (4.*x) - 10.;
并且必须在 main 中更改为以下内容：
ans = secant_method(f, 0., 2.);
我在 Windows 上测试了这个，DECIMAL_DIG 定义为 9。我得到了没有点的 #INF，使用 11 次迭代得到了带点的 5e-315。
printf("\nRoot found:\t%.*g\n", DECIMAL_DIG * 2., ans);
我正在制作一个简单的函数来重现 Secant method，但我认为我在精度上遇到了一些问题。
首先，这是函数（以及调用它的主要方法，以及与之一起使用的测试函数）：
double secant_method(double(*f)(double), double a, double b){
double c;
for (int i = 0; i < 10; i++){
c = a - f(a) * (a - b) / (f(a) - f(b));
b = a; a = c;
}
return c;
}
typedef double (*func)(double);
//test function - x^3 + 4x - 10
double test(double x){
return (x*x*x) + (4*x) - 10;
}
int main(){
func f = &test;
double ans;
ans = secant_method(f, 0, 2);
printf("\nRoot found:\t%.*g\n", DECIMAL_DIG * 2, ans);
return 0;
}
注意：函数secant_method()中的for循环只循环了10次。这就是我的问题所在。
按原样打印时，一切正常。它给出了小数点后 16 位的正确输出：
但是，当我向 secant_method() 中的 for 循环添加迭代时，会发生这种情况：
为什么会这样？我是否达到了 C 可以处理的最大表示？
我通读了另一个 post 的 this great answer，但在其中，关于我收到的异常 (-1.#IND)，它只是说我的结果不是号码，或者我在进行某种非法操作。
编辑： 使用 (x*x*x) + (4*x) - 10 + sin(x) 作为我的测试函数给出了正确的答案 - 但只有当我在 i < 9 而不是 [=18= 时循环] (x*x*x) + (4*x) - 10
-1.#IND是微软输出不确定值的方式，具体来说就是NaN.
一种 可能发生这种情况的方式是 0 / 0 但我会检查 所有 操作以查看问题出在哪里谎言：
double secant_method(double(*f)(double), double a, double b){
double c;
printf("DBG =====\n");
for (int i = 0; i < 10; i++){
printf("\nDBG -----\n");
printf("DBG i: %d\n",i);
printf("DBG a: %30f\n",a);
printf("DBG b: %30f\n",b);
printf("DBG c: %30f\n",c);
printf("DBG f(a): %30f\n",f(a));
printf("DBG a-b: %30f\n",a-b);
printf("DBG f(b): %30f\n",f(b));
printf("DBG f(a)-f(b): %30f\n",f(a)-f(b));
printf("DBG f(a)*(a-b): %30f\n",f(a)*(a-b));
printf("DBG f(a)*(a-b)/(f(a)-f(b)): %30f\n",f(a)*(a-b)/(f(a)-f(b)));
c = a - f(a) * (a - b) / (f(a) - f(b));
b = a; a = c;
}
return c;
}
一旦你有了调试输出，然后你就可以弄清楚实际问题是什么，并采取策略来避免它。
当我这样做时，我看到（最后）：
DBG -----
DBG i: 8
DBG a: 1.556773264394211375716281509085
DBG b: 1.556773264393484179635152031551
DBG c: 1.556773264394211375716281509085
DBG f(a): -0.000000000000000987057657830803
DBG a-b: 0.000000000000727196081129477534
DBG f(b): -0.000000000008196943991622962500
DBG f(a)-f(b): 0.000000000008195956933965131697
DBG f(a)*(a-b): -0.000000000000000000000000000718
DBG f(a)*(a-b)/(f(a)-f(b)): -0.000000000000000087577871187781
DBG -----
DBG i: 9
DBG a: 1.556773264394211375716281509085
DBG b: 1.556773264394211375716281509085
DBG c: 1.556773264394211375716281509085
DBG f(a): -0.000000000000000987057657830803
DBG a-b: 0.000000000000000000000000000000
DBG f(b): -0.000000000000000987057657830803
DBG f(a)-f(b): 0.000000000000000000000000000000
DBG f(a)*(a-b): -0.000000000000000000000000000000
DBG f(a)*(a-b)/(f(a)-f(b)): nan
Root found: nan
所以你可以看到，在第十次迭代中，a 和 b 变得相等，因此 f(a) 和 f(b) 也相等。所以你得到了表达式：
something * 0 / 0
如前所述，这会给你 0 / 0 或 NaN。
关于如何修复它，您只需要避免除以零，因为这会给您带来 NaN 或无穷大。因此，您可以改用以下函数：
double secant_method(double(*f)(double), double a, double b){
double c;
for (int i = 0; i < 1000; i++) {
if (f(a) == f(b)) break;
c = a - f(a) * (a - b) / (f(a) - f(b));
b = a; a = c;
}
return c;
}
一千个循环应该足以得到一个像样的答案，如果你要除以零，它会提前退出。
如果您想要 更高 精度，您可以查看 long double 类型或切换到使用任意精度算术库之一，例如 GMP 或 MPIR。
这通常需要更多的工作，但您可以获得一些令人印象深刻的结果。这个程序，建立在 MPIR 上：
#include <stdio.h>
#include <mpir.h>
void secant_method(mpf_t result, void(*f)(mpf_t, mpf_t), mpf_t a, mpf_t b){
mpf_t c, fa, fb, temp1, temp2;
mpf_init (fa);
mpf_init (fb);
mpf_init (temp1);
mpf_init (temp2);
for (int i = 0; i < 1000; i++){
printf("DBG i: %d\n",i);
f (fa, a);
f (fb, b);
if (mpf_cmp (fa, fb) == 0) break;
mpf_set (temp1, a);
mpf_sub (temp1, temp1, b);
mpf_set (temp2, fa);
mpf_sub (temp2, temp2, fb);
mpf_set (result, fa);
mpf_mul (result, result, temp1);
mpf_div (result, result, temp2);
mpf_sub (result, result, a);
mpf_neg (result, result);
mpf_set (b, a);
mpf_set (a, result);
}
}
void test (mpf_t result, mpf_t x){
mpf_t temp;
mpf_set (result, x);
mpf_pow_ui (result, result, 3);
mpf_init_set (temp, x);
mpf_mul_ui (temp, temp, 4);
mpf_add (result, result, temp);
mpf_set_ui (temp, 10);
mpf_sub (result, result, temp);
mpf_clear (temp);
}
int main(){
mpf_t ans, a, b;
mpf_set_default_prec (8000);
mpf_init_set_ui (ans, 0);
mpf_init_set_ui (a, 0);
mpf_init_set_ui (b, 2);
secant_method (ans, &test, a, b);
mpf_out_str (stdout, 10, 0, ans);
return 0;
}
输出精度更高，大约二分半 千 位：
DBG i: 1
:
DBG i: 19
0.155677326439421146326886324730853302634853266143
22856485101283627988036767055520913212330822780959
93349183787687346999781239000417393618333668026011
02048595843228945228507966189601958673920851932189
20626590635658264390975889008832048255537650792123
54916373054888140164770654992918100928227714960414
65208113116379497717707745267800989233875981344305
90022883167106124203999713536673991376957068731244
91919087980169395013246250812213656324598765244218
15974098310512802880727074335472786858740154287363
31949470951650710072488856623955478366217474755111
76368234254761541647442609230138418167182918204711
66713459423756284737546964906061587903876515793884
14091165347411853670752820576131460960421137744435
73729141652832258144582021037373967987171478026002
48487515446248979731517957120705447608265161099693
33098235693813752370774508652788986557620510981156
19907950657355934071535840759135251701581523712307
00051674680667972152582339710574822560693109306285
91240827697915787078746087225027856691436076089912
35551789799825731841345891629028445554314717823386
07885164744100235567602875364878328805811271289098
87558119684442289199181352023304600847178256323082
57317198584882656089836229208443415369358460418542
84083408696290686178971039756668669303212658278679
39542421457300944206839268283788585029652481323614
65995074020560963212330914882733926627309382310653
39023265929195094492468196461296569155421718696631
73798097369621805062145075113127308161572398104766
37356504104570136778437926442139603916930640425421
15655156674699552536588332891562053247342008145504
44336211031437923307615880759201695011419324719812
46482293928341901673056596202744639074280785106031
90197472588293352508389295101867514582271001202777
85575614897203080940643669476500979934666490279524
88486176409290187337498631681392563044899541391612
88438904336237873504970887963071622208868799638373
42186338496601471274609131141920820263780493617795
89714798662834913192777810386631915415021934333441
01797098172897161215116673422762953435902633516501
73788202968876596925999628999004575114529754782488
59959395407324243559011982543407738505315960009874
36510513519775603567237051670918870105777288994910
85524037720122749091827520695838000086150188462000
63190624219373460624686216781527327604063990319908
56547016812115842640285111265677758613385414834511
69237199199725030839166586376374587900611430229333
87296847315023767826706323911923435564643861604120
017381909481e1
而且，如果你把那个数字传回 test() 函数，你会得到一个非常接近于零的数字，大约 -1.15 x 10<sup>- 2408</sup>。所以我想说这是对使用 double.
而且，就其价值而言，它只需要大约十分之一秒 CPU 的时间，所以至少用任意精度的算法来做到这一点是可行的。
甚至更精度，只需更改MPIR的默认精度设置，目前设置为：
mpf_set_default_prec (8000);
增加到 100,000 会给出超过 30,000 位有效数字的答案，最终 "close-to-zero" 答案约为 -5 x 10<sup>-30103</sup>.
#IND是由零除以零引起的。调试时在循环中插入一个简单的 isnan(c) 检查，你会发现 a 和 b 最终变得相等，这导致 a - b 和 f(a) - f(b) 为零。
我认为选角是个问题。编译器不会将您所有的常数解释为 double，而是将其解释为常规整数。添加一个“。” （或“.0”）在每个常量的末尾使它们加倍。因此 test() 应更改为：
return (x*x*x) + (4.*x) - 10.;
并且必须在 main 中更改为以下内容：
ans = secant_method(f, 0., 2.);
我在 Windows 上测试了这个，DECIMAL_DIG 定义为 9。我得到了没有点的 #INF，使用 11 次迭代得到了带点的 5e-315。
printf("\nRoot found:\t%.*g\n", DECIMAL_DIG * 2., ans);