d6r5.cpp

C语言的编程合集！是我自己编的！很有用的！大家可以下载！有意见麻烦大家提出来以便本人及时修改！

#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
 
static double r[98];
static int iff,ix1,ix2,ix3;
double ran1(int& idum)
{
    double rm1,rm2,t;
	int m1,m2,m3,ia1,ia2,ia3,ic1,ic2,ic3,j;
	 
    m1 = 259200, ia1 = 7141, ic1 = 54773, rm1 = 0.0000038580247,
    m2 = 134456, ia2 = 8121, ic2 = 28411, rm2 = 0.0000074373773,
    m3 = 243000, ia3= 4561, ic3 = 51349;
	
    if ((idum < 0) || (iff == 0))
	{
        iff = 1;
        ix1 =  (ic1 - idum) % m1;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix2 = ix1 % m2;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix3 = ix1 % m3;
        for (j = 1; j<=97; j++)
		{
            ix1 = (ia1 * ix1 + ic1)% m1;
            ix2 = (ia2 * ix2 + ic2)% m2;
            r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
        }
        idum = 1;
	}
      ix1 = (ia1 * ix1 + ic1) % m1;
      ix2 = (ia2 * ix2 + ic2) % m2;
      ix3 = (ia3 * ix3 + ic3) % m3;
      j = 1 + int((97 * ix3) / m3);
      if ((j > 97) || (j < 1)) 
	  {
		  cout<< "abnormal exit"<<endl;
	      return 1;
	  }
      t = r[j];
      r[j] = (float(ix1) + float(ix2) * rm2) * rm1;
	  return t;
}


double expdev(int & idum)
{
	
	double temp;
    double dum;
	dum=0.0;
    while (dum == 0)
        dum = ran1(idum);
    temp = -log(dum);
	return temp;
}

void main()
{
    //program d6r5
    //driver for routine expdev
	int npts,i,j,idum;
	double y;
	double ee,trig[22], x[22],expect,total;
    npts = 10000; ee = 2.718281828;
    for (i = 1; i<=21; i++)
	{
        trig[i] = (i - 1) / 20.0;
        x[i] = 0.0;
    }
    idum = -1;
    for (i = 1; i<=npts; i++)
	{
        y = expdev(idum);
        for (j = 2; j<=21; j++)
		{
            if ((y < trig[j]) && (y > trig[j - 1])) 
			{
                x[j] = x[j] + 1.0;
			}
        }
    }
    total = 0.0;
    for (i = 2; i<=21; i++)
	{
        total = total + x[i];
    }
    cout<<"Exponential distribution with"<<npts<<"points;"<<endl;
    cout<<"      interval       observed    expected"<<endl;
    for (i = 2; i<=21; i++)
	{
        x[i] = x[i] / total;
        expect = exp(-(trig[i - 1] + trig[i]) / 2.0);
        expect = expect * 0.05 * ee / (ee - 1.0);  
	    cout<<setiosflags(ios::fixed);
		cout<<setprecision(2)<<setw(8)<<trig[i - 1];
        cout<<setprecision(2)<<setw(8)<<trig[i];
		cout<<setprecision(4)<<setw(12)<<x[i];
        cout<<setprecision(4)<<setw(12)<<expect<<endl;
	}
}