d6r9.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  gammln(double xx)
{
	int j;
	double temp; 
    double cof[6],stp,half,one,fpf,x,tmp,ser;
    cof[1] = 76.18009173;
    cof[2] = -86.50532033;
    cof[3] = 24.01409822;
    cof[4] = -1.231739516;
    cof[5] = 0.00120858003;
    cof[6] = -0.00000536382;
    stp = 2.50662827465;
    half = 0.5;
    one = 1.0;
    fpf = 5.5;
    x = xx - one;
    tmp = x + fpf;
    tmp = (x + half) * log(tmp) - tmp;
    ser = one;
    for (j = 1;j<=6;j++)
	{
        x = x + one;
        ser = ser + cof[j] / x;
    }
    temp = tmp + log(stp * ser);
	return temp;
}

double bnldev(double pp, int n, int &idum)
{
	int j,en,nold,bnl;
	double p,am,g,t,oldg,pc,plog,pclog,pold,sq,y,em,pi;
    pi = 3.141592654;
	nold=0;
	pold=0.0;
    if (pp <= 0.5)
	    p = pp;
    else
        p = 1.0 - pp;
    am = n * p;
    if (n < 25)
	{
        bnl = 0;
        for (j = 1; j<=n; j++)
            if (ran1(idum) < p) bnl = bnl + 1;
    }
    else
	{
		if (am < 1.0)
		{
        g = exp(-am);
        t = 1.0;
        for (j = 0;j<=n;j++)
			{
				t = t * ran1(idum);
				if (t < g) exit(1);
			}
        if (t >= g) j = n;
        bnl = j;
		}
        else
		{
             if (n!= nold)
					{
					en = n;
					oldg = gammln(en + 1.0);
					nold = n;
					}
			if (p != pold)
			{
	        pc = 1.0 - p;
            plog = log(p);
            pclog = log(pc);
            pold = p;
			}
			sq = sqrt(2.0 * am * pc);
			do
			{
				do
				{
				y = tan(pi * ran1(idum));
				em = sq * y + am;
				}
				while ((em < 0) || (em >= en + 1.0));
				em = int(em);
				t = en - em;
				t = exp(oldg-gammln(em+1.0)-gammln(t+1.0)+em*plog+t*pclog);
				t = 1.2 * sq * (1.0 + y *y) * t;
			}
			while (ran1(idum) > t);
        bnl = int(em);
		}
	}
    if (p != pp) bnl = n - bnl;
    t=bnl;
	return t;
}

void main()
{
    //program d6r9
    //driver for routine bnldev
    int idum,n = 20, npts = 10000, iscal = 200, llen = 50,nn=100,xm,klim,i,j,k,m;
    double pp,dist[21];
	char text[50];
    idum = -133;
    for (j = 1; j<=20; j++)
        dist[j] = 0.0;
    cout<<"Mean of Binomial distribution [n=1..20] -1 to end."<<endl;
    xm = 5;
    cout<<xm<<endl;
    if (xm <= 0) exit(1);
	pp=double(xm)/double(nn);
	for (i = 1; i<=npts; i++)
	{
        j = int(bnldev(pp,nn, idum));
        if ((j >= 0) && (j <= 20)) dist[j+1] = dist[j+1] + 1;
    }
    cout<< "    x     p[x]     graph:"<<endl;
    for (j = 1; j<=20; j++)
	{
        dist[j] = dist[j] / npts;
        for (k = 0; k<=50; k++)
            text[k] =' ';
        klim = int(iscal * dist[j]);
        if (klim > llen) klim = llen;
        for (k = 1; k<=klim; k++)
            text[k] = '*';
        cout<< setw(5)<<float(j-1);
        cout<< setw(10)<<dist[j]<<"    ";
        
        for (m = 1; m<=50; m++)
            cout<< text[m];
		cout<<endl; 
     }
}