📄 remes.htm
字号:
<html><head><meta http-equiv="Content-Type" content="text/html; charset=gb2312"><title></title></head><body><pre id="content" style="width:90%;white-space:normal;word-wrap:break-word;word-break:keep-all;">
#include<stdio.h>
#include<malloc.h>
#include<stdlib.h>
#include<conio.h>
#include<math.h>
void input_begin_end(float*,float*);
void input_n(int*);
void input_e(float*);
void Remes(float,float,int,float);
void compute_polynomial(float*,float*,int);
int compute_newx(float*,float*,float,float,float,int);
int sup(float*,float*,float*,float,float,float,int);
void Bisection(float*,float*,float,float,float,int);
float dkx(float*,float,int,int);
float bisect_dk(float*,float,float,float,int,int);
float difference(float*,float*,int);
float*my_malloc(int);
void main()
{
int n;
float e,begin,end;
clrscr();
printf("\nWelcome to use Remes Algorithm to compute the best uniformly-approximating \n\npolynomial of Exp(x).\n");
input_begin_end(&begin,&end);
input_n(&n);
input_e(&e);
Remes(begin,end,n,e);
printf("\nLots of thanks for your using.\n\nPress any key to quit.\n");
getch();
}
void Remes(float begin,float end,int n,float e)
{
int i;
float error;
float*p1=my_malloc(n+1);
float*p2=my_malloc(n+1);
float*x=my_malloc(n+2);
for(i=0;i<=n+1;i++)
{
*(x+i)=begin+(end-begin)/(n+1)*i;
}
compute_polynomial(p1,x,n);
do{
if(compute_newx(p1,x,begin,end,e,n))
{
compute_polynomial(p2,x,n);
error=difference(p1,p2,n+1);
for(i=0;i<=n;i++)
{
*(p1+i)=*(p2+i);
}
}
else
{
error=0;
}
}while(error>e);
printf("\nThe coefficients of the best %d_degree uniformly-approximating polynomial are as follows: \n",n);
for(i=0;i<=n;i++)
{
printf("\n%d-degree term: %f\n",n-i,*(p1+i));
}
free(p1);
free(p2);
free(x);
}
void compute_polynomial(float*p,float*x,int n)
{
int i,j,k,r;
float y,a,tempa,tempb;
float*A=my_malloc((n+2)*(n+2));
float*b=my_malloc(n+2);
for(i=0;i<=n+1;i++)
{
if(i%2)
{
*(A+(n+2)*i)=-1;
}
else
{
*(A+(n+2)*i)=1;
}
*(A+(n+2)*i+1)=1;
for(j=2;j<=n+1;j++)
{
*(A+(n+2)*i+j)=pow(*(x+i),j-1);
}
}
for(i=0;i<=n+1;i++)
{
*(b+i)=exp(*(x+i));
}
for(k=1;k<=n+1;k++)
{
r=k;
a=fabs(*(A+(n+2)*(k-1)+(k-1)));
for(i=k+1;i<=n+2;i++)
{
if(a<fabs(*(A+(n+2)*(i-1)+(k-1))))
{
r=i;
a=fabs(*(A+(n+2)*(i-1)+(k-1)));
}
}
if(r!=k)
{
tempb=*(b+(k-1));
*(b+(k-1))=*(b+(r-1));
*(b+(r-1))=tempb;
for(i=1;i<=n+2;i++)
{
tempa=*(A+(n+2)*(k-1)+(i-1));
*(A+(n+2)*(k-1)+(i-1))=*(A+(n+2)*(r-1)+(i-1));
*(A+(n+2)*(r-1)+(i-1))=tempa;
}
}
for(i=k+1;i<=n+2;i++)
{
*(A+(n+2)*(i-1)+(k-1))=*(A+(n+2)*(i-1)+(k-1))/(*(A+(n+2)*(k-1)+(k-1)));
*(b+(i-1))=*(b+(i-1))-*(A+(n+2)*(i-1)+(k-1))*(*(b+(k-1)));
for(j=k+1;j<=n+2;j++)
{
*(A+(n+2)*(i-1)+(j-1))=*(A+(n+2)*(i-1)+(j-1))-*(A+(n+2)*(i-1)+(k-1))*(*(A+(n+2)*(k-1)+(j-1)));
}
*(A+(n+2)*(i-1)+(k-1))=0;
}
}
for(i=n+2;i>=1;i--)
{
y=0;
for(j=n+2;j>i;j--)
{
y=y+*(b+(j-1))*(*(A+(n+2)*(i-1)+(j-1)));
}
*(b+(i-1))=(*(b+(i-1))-y)/(*(A+(n+2)*(i-1)+(i-1)));
}
for(i=0;i<=n;i++)
{
*(p+i)=*(b+(n+1-i));
}
free(A);
free(b);
}
int compute_newx(float*p,float*x,float begin,float end,float e,int n)
{
int i;
float xm;
if(sup(p,x,&xm,begin,end,e,n)==0)
{
return 0;
}
else
{
for(i=0;i<=n+1;i++)
{
if(xm<*(x+i))
{
break;
}
}
if(i==0)
{
if((dkx(p,xm,n,0)*dkx(p,*x,n,0))>0)
{
*x=xm;
}
else
{
for(i=n+1;i>=1;i--)
{
*(x+i)=*(x+i-1);
}
*x=xm;
}
}
else if(i==n+2)
{
if((dkx(p,xm,n,0)*dkx(p,*(x+n+1),n,0))>0)
{
*(x+n+1)=xm;
}
else
{
for(i=0;i<=n;i++)
{
*(x+i)=*(x+i+1);
}
*(x+n+1)=xm;
}
}
else
{
if((dkx(p,xm,n,0)*dkx(p,*(x+i),n,0))>0)
{
*(x+i)=xm;
}
else
{
*(x+i-1)=xm;
}
}
return 1;
}
}
int sup(float*p,float*x,float*pm,float begin,float end,float e,int n)
{
int i,j,num;
float*partition=my_malloc(n+2);
float*abs_value=my_malloc(n+2);
Bisection(p,partition,begin,end,e,n);
for(i=0;i<=n+1;i++)
{
*(abs_value+i)=fabs(dkx(p,*(partition+i),n,0));
}
for(i=1;i<=n+1;i++)
{
if(*(abs_value+i)>*(abs_value))
{
num=1;
*partition=*(partition+i);
*abs_value=*(abs_value+i);
}
else if(*(abs_value+i)==*(abs_value))
{
num++;
*partition=*(partition+i);
*abs_value=*(abs_value+i);
}
}
for(i=0;i<num;i++)
{
for(j=0;j<=n+1;j++)
{
if(*(x+j)==*(partition+i))
{
break;
}
}
if(j==n+2)
{
*pm=*(partition+i);
free(partition);
free(abs_value);
return 1;
}
}
free(partition);
free(abs_value);
return 0;
}
void Bisection(float*p,float*partition,float begin, float end,float e,int n)
{
int i,j,k;
float a,b,dka,dkb;
float*r=my_malloc(n+1);
for(i=0;i<=n+1;i++)
{
if(i==0)
{
*partition=begin;
}
else
{
*(partition+i)=end;
}
}
for(k=n;k>=1;k--)
{
j=0;
for(i=1;i<=n+1;i++)
{
a=*(partition+i-1);
b=*(partition+i);
dka=dkx(p,a,n,k);
dkb=dkx(p,b,n,k);
if(fabs(dka)<=e)
{
*(r+j)=a;
j++;
}
else if(dka*dkb<0)
{
*(r+j)=bisect_dk(p,a,b,e,n,k);
j++;
}
}
for(i=1;i<=j;i++)
{
*(partition+i)=*(r+i-1);
}
}
free(r);
}
float dkx(float*p,float x,int n,int k)
{
float dkx;
int i,j;
float*q=my_malloc(n+1);
for(i=0;i<=n;i++)
{
*(q+i)=*(p+i);
}
for(i=0;i<=k;i++)
{
for(j=1;j<=n-i;j++)
{
*(q+j)=*(q+j)+*(q+j-1)*x;
}
}
dkx=*(q+n-k);
free(q);
for(i=2;i<=k;i++)
{
dkx=dkx*i;
}
dkx=exp(x)-dkx;
return dkx;
}
float bisect_dk(float*p,float a,float b,float e,int n,int k)
{
float dka,dkb,dkc,c;
while(1)
{
c=(a+b)/2;
if(b-a<=e)
{
return c;
}
else
{
dka=dkx(p,a,n,k);
dkb=dkx(p,b,n,k);
dkc=dkx(p,c,n,k);
if(fabs(dka)<=e)
{
return a;
}
else if(fabs(dkb)<=e)
{
return b;
}
else if(dka*dkc<0)
{
b=c;
}
else
{
a=c;
}
}
}
}
float difference(float*vector1,float*vector2,int k)
{
int i;
float difference=0;
for(i=2;i<=k-1;i++)
{
difference=difference+fabs(*(vector1+i)-*(vector2+i));
}
return difference;
}
void input_n(int*pn)
{
printf("\nPlease input the degree.\"Enter\"to confirm.\n\nn=");
scanf("%d",pn);
while((*pn)<0)
{
printf("\nSorry,your inputting is illeagle.Please input again.\n\nA positive integer is expected.\n\nn=");
scanf("%d",pn);
}
}
void input_e(float*pe)
{
printf("\nPlease input the lagest error you would allow.\"Enter\"to confirm.\n\ne=");
scanf("%f",pe);
while((*pe)<=0)
{
printf("\nSorry,your inputting is illeagle.Please input again.\n\nA positive number is expected.\n\ne=");
scanf("%f",pe);
}
}
void input_begin_end(float*pbegin,float*pend)
{
printf("\nPlease input the beginning point of the interval you require.\"Enter\"to confirm.\n");
printf("\nBeginning point=");
scanf("%f",pbegin);
printf("\nPlease input the ending point of the interval you require.\"Enter\"to confirm.\n");
printf("\nEnding point=");
scanf("%f",pend);
while(*pbegin>=*pend)
{
printf("\nSorry,your inputting is illeagle.Pplease input again.\n\nThe beginning point shoud be strictly smller than the endding one.\n");
printf("\nPlease input the beginning point of the interval you require.\"Enter\"to confirm.\n");
printf("\nBeginning point=");
scanf("%f",pbegin);
printf("\nPlease input the ending point of the interval you require.\"Enter\"to confirm.\n");
printf("\nEnding point=");
scanf("%f",pend);
}
}
float* my_malloc(int k)
{
float*q;
q=(float*)malloc(4*k);
if(q)
{
return q;
}
else
{
printf("\nSorry,the program cannot run because Memory Request Is Rejected.\n");
printf("\nLots of thanks for your using.\n\nPress any key to quit.\n");
getch();
exit(-1);
}
}
</pre></body></html>⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -