#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>

void main()
{
    //program d6r7
    //driver for routine gamdev
    int idum,n = 20, npts = 10000, iscal = 200, llen = 50, ia;
	int 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;
     }
}