ecm算法.txt

ECM算法在数据丢失的情况下

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

void main()
{
	long double simprl1(long double x0,long double x1,int m,long double n);
	long double partialfun(long double p);
	long double lagam(long double x);//完全gamma函数
	long double lbgam(long double a,long double x) ;//不完全gamma函数
	long double a=1;//初始化参数a
	long double b=4;//初始化参数b
	long double c=4;//门限值
	long double gi_sum=0,yi_sum=0;//对前r个gi求和,对前r个yi求和
	long double sum_gi,sum_yi,gsum,ysum;
	long double b_temp,a_temp;

	long double yi;
	long double yi_b[100000];
	long int n=0,r=0;//总数据个数；有效数据个数;
	long int m=12;
	float dis;

	
	FILE *Test_File,*Result_File,*Result_File1;

	
//	if((Test_File=fopen("E:\\VC\\ECM\\标准数据\\9-27\\data.txt","r"))==NULL)
//		printf("Can't open Test_File.\n");

	if((Test_File=fopen("E:\\VC\\ECM\\data.txt","r"))==NULL)
		printf("Can't open Test_File.\n");

	if((Result_File=fopen("E:\\VC\\ECM\\result.txt","w"))==NULL)
		printf("Can't open Result_File1.\n");

	if((Result_File1=fopen("C:\\MATLAB7\\work\\data.txt","w"))==NULL)
		printf("Can't open Result_File1.\n");

/*	while(!feof(Test_File))
	{
		fscanf(Test_File,"%f %lf",&dis,&yi);
		if(dis<30&&dis>=20)
		{
			if(yi!=0)
			{
				yi_b[r]=yi+85;
				fprintf(Result_File1,"%f\n",yi_b[r]);
				r++;				
			}
			n++;
		}
	}
	*/
	while(!feof(Test_File))
	{
		fscanf(Test_File,"%lf",&yi);
		if(yi!=0)
		{
			yi_b[r]=yi;
			r++;
		}
		n++;
	}
	double pi=double((n-r))/n;
	//fprintf(Result_File,"%f\n",pi);
	printf("%f\n%d\n%d\n",pi,n,r);
	for(long int i=0;i<r-1;i++)
	{
		gi_sum=gi_sum+log(yi_b[i]);//对前r个gi求和
		yi_sum=yi_sum+yi_b[i];//对前r个yi求和
	}

	long double Q1=0;
	long double Q2;
	long double Q_temp;
	long double bvar,avar;
	i=1;
	while (1)
	{
		sum_gi=(n-r)*(log(b)+simprl1(0.00000000000000001,c/b,m,a)/lbgam(a,c/b));
		sum_yi=(n-r)*(b*lbgam(a+1,c/b)/lbgam(a,c/b));
		gsum=gi_sum+sum_gi;
		ysum=yi_sum+sum_yi;
		long double	p=gsum/n-log(b);
		b_temp=ysum/(n*a);
		a_temp=partialfun(p);
		b=b_temp;
		a=a_temp;
		Q2=-n*log(lagam(a))-n*a*b+(a-1)*gsum-ysum/b;
		Q_temp=Q2-Q1;
		if(Q_temp<0)
			Q_temp=-Q_temp;
		if((Q_temp)<5.0e-9)
			break;
		i++;
		Q1=Q2;
		avar=1/(n*simprl1(0.00000000000000001,10000000000000000000,m,a));
		bvar=1/(2*ysum/(b*b*b)-n*a/(b*b));
		if((Q_temp)<3.0e-2)			
			fprintf(Result_File,"%18.10Lf\t%18.10Lf\t%18.10Lf\t%18.10Lf\n",a,b,avar,bvar);
	}
	
	fprintf(Result_File,"%18.10Lf\t%18.10Lf\t%18.10Lf\t%18.10Lf\n",a,b,Q_temp);
	fclose(Test_File);
	fclose(Result_File);
}

long double simprl1(long double x0,long double x1,int m,long double n)
{
	long double h=(x1-x0)/(2*m);
	long double s1=0;
	long double s2=0;
	long double x,s;
	long double myfun1(long double x,long double n);
	long double myfun2(long double x,long double n);
	for(int k=1;k<=m;k++)
	{
		x=x0+h*(2*k-1);
		s1=s1+myfun1(x,n);
	}
	for(k=1;k<m;k++)
	{
		x=x0+h*2*k;
		s2=s2+myfun1(x,n);
	}
	s=h*(myfun1(x0,n)+myfun1(x1,n)+4*s1+2*s2)/3;
	return s;
}

long double myfun1(long double x,long double n)
{
	return log(x)*pow(x,n-1)*exp(-x);
}

long double myfun2(long double x,long double n)
{
	return pow(x,n-2)*exp(-x)+(n-1)*(n-2)*pow(x,n-3)*exp(-x);
}
/*解方程log(gamma(z))-p(z-1)=0*/
long double partialfun(long double p)
{
	long double x0=1.00000000000000000001,x1=15;
	long double lagam(long double x);
	long double result1=lagam(x1)-p*(x1-1),result2=lagam(x0)-p*(x0-1);
	long double x2=x1-result1*(x1-x0)/(result1-result2);  
	long int n=1;
	long double x3=x1-x0;
	while ((x3>=1.0e-4)&&(n<=1000000))     //判断两个条件截止
    {
		x0=x1;                                       //将x1赋给x0
		x1=x2;  //将x2赋给x1
		x3=x1-x0;
		if(x3<0)
			x3=-x3;
		result1=lagam(x1)-p*(x1-1);
		result2=lagam(x0)-p*(x0-1);
		x2=x1-result1*(x1-x0)/(result1-result2);    //迭代运算
		n=n+1;
	}
    return x2;
}

/*计算伽马(gamma)函数值,gamma函数积分区间为0到正无穷*/
long double lagam(long double x)
{ 
int i;
long double y,t,s,u;
static   long double a[11]={ 0.0000677106,-0.0003442342,
  0.0015397681,-0.0024467480,0.0109736958,
  -0.0002109075,0.0742379071,0.0815782188,
  0.4118402518,0.4227843370,1.0};
if (x<=0.0)
{ 
  printf("err**x<=0!\n"); 
  return(-1.0);
}
y=x;
if (y<=1.0)
{ 
  t=1.0/(y*(y+1.0)); 
  y=y+2.0;
}
else if (y<=2.0)
{ 
  t=1.0/y; 
  y=y+1.0;
}
else if (y<=3.0) 
{
  t=1.0;
}
else
{ 
  t=1.0;
  while (y>3.0)
  { 
   y=y-1.0; 
   t=t*y;
  }
}
s=a[0]; 
u=y-2.0;
for (i=1; i<=10; i++)
{
  s=s*u+a[i];
}
s=s*t;
return(s);
}

/*计算不完全伽马(gamma)函数值,gamma函数积分区间为0到a*/
long double lbgam(long double a,long double x) 
{ 
int n;
    long double p,q,d,s,s1,p0,q0,p1,q1,qq;
    if ((a<=0.0)||(x<0.0))
{ 
  if (a<=0.0) 
  {
   printf("err**a<=0!\n");
  }
        if (x<0.0) 
  {
   printf("err**x<0!\n");
  }
        return(-1.0);
}
    if (x+1.0==1.0) 
{
  return(0.0);
}
    if (x>1.0e+35) 
{
  return(1.0);
}
    q=log(x); 
q=a*q; 
qq=exp(q);
    if (x<1.0+a)
{ 
  p=a; 
  d=1.0/a; 
  s=d;
        for (n=1; n<=100; n++)
  { 
   p=1.0+p;
   d=d*x/p; 
   s=s+d;
   if (fabs(d)<fabs(s)*1.0e-07)
   { 
    s=s*exp(-x)*qq/lagam(a);
                return(s);
   }
  }
}
    else
{ 
  s=1.0/x; 
  p0=0.0; 
  p1=1.0; 
  q0=1.0; 
  q1=x;
        for (n=1; n<=100; n++)
  { 
   p0=p1+(n-a)*p0; 
   q0=q1+(n-a)*q0;
            p=x*p0+n*p1; 
   q=x*q0+n*q1;
            if (fabs(q)+1.0!=1.0)
   { 
    s1=p/q; 
    p1=p; 
    q1=q;
                if (fabs((s1-s)/s1)<1.0e-07)
    { 
     s=s1*exp(-x)*qq/lagam(a);
                    return(1.0-s);
    }
                s=s1;
   }
            p1=p; 
   q1=q;
  }
}
    printf("a too large !\n");
    s=1.0-s*exp(-x)*qq/lagam(a);
    return(s);
}