/* Generating random numbers on the [0,1] interval,
	according to the g(x)=3 x^2 distribution function  */

#include <dos.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>


void plot(int x, int y, int color);
void setmode (int mode);
void cls(int color);


int xres=640;
int yres=480;

float lambdar, lambdai;

long int A[102];

main()
{ float x;
  long int i;
  int is, js;
  long int Amax;
  Amax=0;
  for(is=1; is<=101; is++) A[is]=0;
  for(i=1; i<=1000000; i++)
  {x=pow((float)(rand())/(RAND_MAX+1.0),0.333);
	A[(int)(x/0.01)]++;
	if (A[(int)(x/0.01)]>Amax) Amax=A[(int)(x/0.01)];
  }
  setmode(18);
  cls(7);

  for (is=0 ;is<=580; is++)
	  for (js=480; js>=480-(float)(A[is/6])/(float)(Amax)*400.0; js--)
		{
		  plot(is, js , 12);
		}

getchar();
}



void setmode(int mode)
{
  union REGS reg;
  reg.x.ax=mode;
  int86 (0x10,&reg,&reg);
}



void cls(int color)
{
    union REGS reg;

    reg.x.ax=0x0600;
    reg.x.cx=0;
    reg.x.dx=0x1E4F;
    reg.h.bh=color;
    int86(0x10,&reg,&reg);
}

void plot(int x, int y , int color)
{

   #define graph_out(index,val) {outp(0x3CE,index); outp(0x3CF,val);}
   int dummy,mask;
   char far * address;
   address = (char far *) 0xA0000000L+(long)y*xres/8L+((long)x/8L);
   mask=0x80 >> (x%8);
   graph_out(8,mask);
   graph_out(5,2);
   dummy= *address;
   *address=color;
   graph_out(5,0);
   graph_out(8,0xFF);
 }