root_618.c

用黄金数0.618作为3折求根逐步缩小数据范围的完整程序

/*黄金分割,找可跃函数f(x)的极低点i,值f(i);两分法求根;求值
初目标段,长1始0,插t1,t2点,[0,t1]长L1,[0,t2]长L2=(根号5-1)/2=.618
若f(t1)<f(t2),则i囿[0,t2],否则囿[t2,1]:弃左[0,t1]或右[t2,1],余2段
作新目标:目标长=L1(=(3-根号5)/2)+L2,迭代至<=FLT_EPSILON
新目标长L2:弃[t2,1],始0,长*L2+始=老L1;弃[0,t1],始L1,长*L1+始=老L2*/

#include <math.h>
#include <float.h>
#include <stdio.h>
#include <string.h>
#include <process.h>

#define	P1	31
//[0,4]x^3-2x-5:y(sqrt(6)/3=.8165)=-6.088667,y(.816496)=-6.088662;y(2.094552)=1e-6
#define	L2	.618034
#define	L1	(1-L2)
#define t2f2	f2=DEP(t2=L+W*L2)
#define	t1f1	f1=DEP(t1=L+W*L1);
#define	L_s		L=t2;
#define	R_s		W=t2;

unsigned char b;
double C[P1+1],L,W,t1=0,t2,f1,f2;

double DEP(double t){
double y=C[0];
for(b=1;b!=P1+1;b++)
y=y*t+C[b];
return y;
}

void prt(char*p,double z){
L+=z;
printf("\nFLT_EPSILON=%e;%s:%e=%e",FLT_EPSILON,p,L,DEP(L));
if('('==*p)printf(",(_R)=%e",DEP(W));
exit(0);
}

void main(int ac,char *av[]){
printf("[]_solo_min{Ci=dbl,poke +0:root}y=av1*x^%02d+..+C%02d;rng_L _R;when _L=_R:wit y(_L)",P1,P1);

  	if(1+P1+1+2==ac){
		for(b=0;b!=P1+1;b++){
			if('+'==*av[1+b])
				ac=0;
			C[b]=atof(av[1+b]);
		}

		L=atof(av[1+P1+1]);
		if(L==(W=atof(av[1+P1+1+1])))
			prt("val",0);

		if(ac)
			for(W=W-L,t2f2,t1f1;){
				if(W<=FLT_EPSILON)
					prt("min",FLT_EPSILON/2);

				if(f1<f2){//弃[t2,1],新t2叠老t1
					W=t2-L;
					t2=t1;
					f2=f1;
					t1f1
				}else{//弃[0,t1],新t1叠老t2
					W=L+W-t1;
					L=t1;
					t1=t2;
					f1=f2;
					t2f2;
				}
			}
		else
			if((f1=DEP(L))*DEP(W)<0){
				if(f1>0)
					ac=1;
				for(;;){
					if(W-L<=FLT_EPSILON)
						prt("root",(W-L)/2);

					if(f2=DEP(t2=(L+W)/2))
						if(f2>0)
							if(ac)//同符号
								L_s
							else
								R_s
						else
							if(ac)
								R_s
							else
								L_s
					else{
						L=t2;
						prt("root",0);
					}
				}
			}else
				prt("(_L)",0);
	}
}