#include<stdio.h>#include<math.h>#include<stdlib.h>
//连乘表达式double function(double x){ int i; double s=1; for(i=1;i<21;i++) s*=(x-i); return(s);}
double diff(double x){ int i; double s=0; for(i=1;i<21;i++) s+=function(x)/(x-i); return(s);}
int qr(long double *a,int n,long double *u,long double *v,long double eps,int itmax){ int i,j,k,ii,jj,kk; long double x,y,p,q,r; long double q00,q01,q02,q11,q12,q22; /* 用于求Q*/ int is1,is2,n1; long double *a1; long double b,c,s; if(itmax == 0) /* 已经不能再迭代*/ { printf("fail/n"); return(0); } if(n==1) /* 矩阵是1阶*/ { u[0] = a[0]; v[0] = 0.0; return(1); } if(n==2) { b = (a[0]+a[3]); c = a[0]*a[3]-a[1]*a[2]; s = b*b-4.0*c; y = sqrt(fabs(s)); if (s > 0.0) { if(b > 0.0) u[0] = (b+y)/2.0; else u[0]=(b-y)/2.0; v[0]=0.0; u[1]=c/u[0]; v[1]=0.0; } else { u[0] = b/2.0; v[0] = y/2.0; u[1] = u[0]; v[1] = -v[0]; } return(1); } is1 = 0; is2 = 0; while(is2 < n-1) { is2++; j = is2*n+is2; if(fabs(a[j-1]) < eps*(fabs(a[j-n-1])+fabs(a[j]))) { n1 = is2-is1; a1 = (long double*)malloc(n1*n1*sizeof(long double)); for(i=0; i<n1; i++) for(j=0; j<n1; j++) a1[i*n1+j] = a[(i+is1)*n+j+is1]; qr(a1,n1,u+is1,v+is1,eps,itmax); free(a1); is1 = is2; } } if(is1>0 && is1<n) { n1=n-is1; a1=(long double*)malloc(n1*n1*sizeof(long double)); for(i=0; i<n1; i++) for(j=0; j<n1; j++) a1[i*n1+j] = a[(i+is1)*n+j+is1]; qr(a1,n1,u+is1,v+is1,eps,itmax); free(a1); return(1); } else if(is1 == n) return(1); for (k=0; k<n-1; k++) { if(k==0) { x = a[(n-2)*n+n-2]+a[n*n-1]; y = a[(n-2)*n+n-2]*a[n*n-1]-a[(n-2)*n+n-1]*a[(n-1)*n+n-2]; p = a[0]*(a[0]-x)+a[1]*a[n]+y; q = a[n]*(a[0]+a[n+1]-x); r = a[n]*a[2*n+1]; } else { p = a[k*n+k-1]; q = a[(k+1)*n+k-1]; if(k != n-2) r= a[(k+2)*n+k-1]; else r= 0.0; } if ((fabs(p)+fabs(q)+fabs(r))!=0.0) { if(p<0.0) s = -sqrt(p*p+q*q+r*r); else s = sqrt(p*p+q*q+r*r); if(k!=0) a[k*n+k-1]=-s; q00 = -p/s; q01 = -q/s; q02 = -r/s; q11 = -q00-q02*r/(p+s); q12 = q01*r/(p+s); q22 = -q00-q01*q/(p+s); for(j=k; j<n; j++) { ii = k*n+j; jj = (k+1)*n+j; kk = (k+2)*n+j; p = q00*a[ii]+q01*a[jj]; q = q01*a[ii]+q11*a[jj]; r = q02*a[ii]+q12*a[jj]; if (k!=n-2) { p = p+q02*a[kk]; q = q+q12*a[kk]; r = r+q22*a[kk]; a[kk] = r; } a[ii] = p; a[jj] = q; } j=k+3; if(j>=n-1) j=n-1; for(i=0; i<=j; i++) { ii = i*n+k; jj = i*n+k+1; p = q00*a[ii]+q01*a[jj]; q = q01*a[ii]+q11*a[jj]; r = q02*a[ii]+q12*a[jj]; if(k!=n-2) { kk = i*n+k+2; p = p+q02*a[kk]; q = q+q12*a[kk]; r = r+q22*a[kk]; a[kk]=r; } a[ii]=p; a[jj]=q; } } if(k > 0) { a[(k+1)*n+k-1] = 0.0; if(k != n-2) a[(k+2)*n+k-1] = 0.0; } } i = qr(a,n,u,v,eps,itmax-1); return(i);}
void main(){ int i,j,n=20; double x0,x[20]; long double eps=1.0e-12; //eps 精度要求,用于判断元素是否为0 int itmax=50000; itmax 最大迭代次数 long double u[20],v[20]; u 保存特征值的实部. long double *H; long double b[21]={1,-2.09999999880791e2,20615,-1.256850e6,5.3327946e7, -1.672280820e9,4.0171771630e10,-7.56111184500e11,1.1310276995381e13, -1.35585182899530e14,1.307535010540400e15,-1.0142299865511500e16,6.3030812099294900e16, -3.11333643161391000e17,1.206647803780370000e18,-3.599979517947610000e18,8.037811822645050000e18, -1.2870931245151000000e19,1.3803759753640700000e19,-8.752948036761600000e18,2.432902008176640000e18}; for(i=0;i<20;i++) { x[i]=i+1.1; } printf("原方程的根如下:/n"); for(j=0;j<20;j++) { do { x0=x[j]; x[j]=x0-function(x0)/diff(x0); }while(fabs(x[j]-x0)>1e-8); printf("x[%d]=%lf/n",j,x[j]); } printf("方程的19次系数做微小扰动+pow(2,-23)后的解如下:/n"); H=(long double*)malloc(sizeof(long double)*n*n); /* 生成的上H矩阵*/ for(i=0;i<n;i++) H[i]=-1.0*b[i+1]/b[0]; /* 第一行*/ for(i=n;i<n*n;i++) H[i]=0.0; /* 下面为0*/ for(i=1;i<n;i++) H[i*n+i-1]=1.0; /* 次对角线为1*/ qr(H,n,u,v,eps,itmax); for(j=19;j>=0;j--) { if(v[j]>0) printf("x[%d]=%.12lf+%.12lfi/n",20-j,u[j],v[j]); else if(v[j]<0) printf("x[%d]=%.12lf%.12lfi/n",20-j,u[j],v[j]); else printf("x[%d]=%.12lf/n",20-j,u[j]); }}