d6r7.cpp

这个是常用的数值算法用VC编写的。相信对大家有用哦。

#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 gamdev(int ia, int &idum)
{
	double x,v1,v2,y,s,e,t;
	int j,am;
    if (ia < 1)
	{
		cout<< "abnormal exit";
		exit(1);
	}
    if (ia < 6)
	{
        x = 1.0;
        for (j = 1; j<=ia; j++)
            x = x * ran1(idum);
        x = -log(x);
	}
    else
	{
      do
	  {
        do
		{
          do
		  {
            v1 = 2.0 * ran1(idum) - 1.0;
            v2 = 2.0 * ran1(idum) - 1.0;
		  }
          while (v1*v1 + v2*v2 > 1.0);
          y = v2 / v1;
          am = ia - 1;
          s = sqrt(2.0* am + 1.0);
          x = s * y + am;
		}
        while (x <= 0.0);
        e = (1.0 + y *y) * exp(am * log(x / am) - s * y);
	  }
      while (ran1(idum) > e);
	}
    t = x;
	return t;
}

void main()
{
    //program d6r7
    //driver for routine gamdev
    int idum,n = 20, npts = 10000, iscal = 200, llen = 50,ia,klim,i,j,k,m;
    double dist[22];
	char text[51];
    idum = -13;
    do
	{
        for (j = 0; j<=21; j++)
            dist[j] = 0.0;
        cout<<"Order of gamma distribution [n=1..20] -1 to end."<<endl;
        ia = 4;
        cout<<ia<<endl;
        if (ia <= 0) exit(1);
	}
    while (ia > 20);
    for (i = 1; i<=npts; i++)
	{
        j = int(gamdev(ia, idum)) + 1;
        if ((j >= 1) && (j <= 21)) dist[j] = dist[j] + 1;
    }
    cout<< "Gamma-distribution deviate, order "<<ia<<"of "<< npts<< " points"<<endl;
    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(4)<<float(j);
        cout<< setw(10)<<dist[j]<<"    ";
        for (m = 1; m<=50; m++)
            cout<< text[m];
		cout<<endl;
     }
}