我正在使用 Maehly 程式來完善多項式的根,并偶然發現了一些有趣的東西:完全相同的代碼給了我兩個非常不同的輸出,具體取決于編譯它的機器。
代碼
#include <stdio.h>
#define MAX_ITERATION 1000
double poly(double x){
double coeff[9]={-61.688, 72.5235, 72.822, -108.519, -5.12949, 39.9139,-7.07373, -3.91823, 1.0};
double result=coeff[0];
double buffer;
for(int i=1; i<9;i ){
buffer=coeff[i];
for(int j=1;j<=i;j ){
buffer*=x;
}
result =buffer;
}
return result;
}
double poly_der(double x){
double coeff[8]={ 72.5235, 72.822, -108.519, -5.12949, 39.9139,-7.07373, -3.91823, 1.0};
double result=coeff[0];
double buffer;
for(int i=1; i<8;i ){
buffer=coeff[i]*(i 1);
for(int j=1;j<=i;j ){
buffer*=x;
}
result =buffer;
}
return result;
}
int main(){
double roots[8]={0.9, -1.1, 1.4, 1.4, -2.0, -2.0, 2.2, 2.2};
double factor;
double pol_eval;
//Implement Maehly-procedure
for(int i=0; i<MAX_ITERATION;i ){
for(int k=0;k<8;k ){
factor=0;
for(int j=0;j<k;j ){
factor =1/(roots[k]-roots[j]);
}
pol_eval=poly(roots[k]);
roots[k]-=pol_eval/(poly_der(roots[k])-(pol_eval*factor));
}
}
for(int i=0;i<8;i ){
printf("\n%d: x:%0.16f poly:%e \n",i,roots[i],poly(roots[i]));
}
}
Windows 輸出 (Windows10):
0: x:1.0072928773885637 poly:-8.437695e-015
1: x:-1.0004044550991309 poly:-2.375877e-014
2: x:1.3770602924650244 poly:-3.552714e-015
3: x:-2.5000428878301499 poly:0.000000e 000
4: x:-1.7318124315476966 poly:-1.136868e-013
5: x:3.0001628929552053 poly:9.094947e-013
6: x:2.2341265341600458 poly:-2.273737e-013
7: x:3.0001628929552049 poly:0.000000e 000
Linux 輸出(Debian GNU/Linux 10):
0: x:1.0072928773885637 poly:-8.437695e-15
1: x:-1.0004044550991309 poly:-2.375877e-14
2: x:1.3770602924650244 poly:-3.552714e-15
3: x:-2.5000428878301499 poly:0.000000e 00
4: x:-1.7318124315476959 poly:2.842171e-14
5: x:3.0001628929552093 poly:-1.818989e-12
6: x:2.2341265341600458 poly:-2.273737e-13
7: x:1.5318471775081237 poly:0.000000e 00
x 是多項式的拋光根,起始值保存在陣列中roots[8]。
你能幫我解釋一下這種行為嗎,最重要的是,幫助我了解如何避免將來發生類似的事情?
uj5u.com熱心網友回復:
我們有2個問題,為什么不同?哪個更好 - 或者可能是第三種方式?
OP報告FLT_EVAL_METHOD2 和 0 的不同值。這意味著 Windows 版本使用long double數學進行中間計算,而 Linux 版本僅使用double. 通常FLT_EVAL_METHOD==2是更正確的。
#include <float.h>
printf("%d\n", FLT_EVAL_METHOD);
FP 在減去幾乎相同的值方面有一個弱點。許多公共位的取消導致輕微的計算錯誤變得更加嚴重。由于各種原因,不同的編譯可能會有略有不同的舍入結果,盡管我 _expect 相同。 poly()通過重復計算來找到 的冪,x然后添加抵消的項。
在我的系統上FLT_EVAL_METHOD是 0。
修復poly()和poly_dev()
static const double coeff[9] = {-61.688, //
72.5235, 72.822, -108.519, -5.12949, //
39.9139, -7.07373, -3.91823, 1.0};
double poly(double x) {
double result = 0.0;
for (int i = 9; i-- > 0; ) {
result = result * x coeff[i];
}
return result;
}
double poly_der(double x) {
double result = 0.0;
for (int i = 9; i-- > 1; ) {
result = result * x coeff[i]*i;
}
return result;
}
只有double數學的結果可以與使用數學的較弱的計算代碼相媲美long double。
printf("%d: x:% 0.16f poly:% e \n", i, roots[i], poly(roots[i]));
0: x: 1.0072928773885645 poly: 0.000000e 00
1: x:-1.0004044550991309 poly:-7.105427e-15
2: x: 1.3770602924650324 poly:-4.973799e-14
3: x:-2.5000428878301495 poly:-6.394885e-13
4: x:-1.7318124315476964 poly:-2.842171e-14
5: x: 1.5318471775081377 poly: 0.000000e 00
6: x: 2.2341265341600520 poly: 0.000000e 00
7: x: 3.0001628929552009 poly:-2.771117e-13
轉載請註明出處,本文鏈接:https://www.uj5u.com/ruanti/383278.html
上一篇:如何縮短此代碼?-C編程
下一篇:使用scanf掃描表C
