/*
 * estimatepi: produce an estimate of pi every display
 */

#include <glut.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>

typedef GLfloat point2[2];
point2 *p,*reds,*greens;
int np,npmax;
#define NPINC 100000

void display() {
 int i,nin=0,nout=0;
 for(i=0;i<np;i++) {
  p[i][0] = drand48()*2.0 - 1.0;
  p[i][1] = drand48()*2.0 - 1.0;
  }
 for(i=0;i<np;i++) {
  float r=p[i][0]*p[i][0]+p[i][1]*p[i][1];
  if (r>=1.0)  {
   reds[nout][0]=p[i][0];
   reds[nout][1]=p[i][1];
   nout++; 
   } else {
   greens[nin][0]=p[i][0];
   greens[nin][1]=p[i][1];
   nin++; 
   }
  }
 glClear(GL_COLOR_BUFFER_BIT);
 glColor3f(0.7,0.0,0.0);
 glBegin(GL_POINTS);
 for(i=0;i<nout;i++) 
  glVertex2fv(reds[i]);
 glEnd();
 glColor3f(0.0,0.7,0.0);
 glBegin(GL_POINTS);
 for(i=0;i<nin;i++) 
  glVertex2fv(greens[i]);
 glEnd();
 glFlush();
 fprintf(stderr,"Estimate: %f\n",(4.0*nin)/(1.0*(nin+nout)));
 }

void init() {
 int i;
 glClearColor(0.0,0.0,0.0,0.0);
 glMatrixMode(GL_PROJECTION);
 glLoadIdentity();
 gluOrtho2D(-1.0,1.0,-1.0,1.0);
 glMatrixMode(GL_MODELVIEW);
 npmax = NPINC;
 p = (point2 *)malloc(npmax*sizeof(point2));
 if (!p) { perror("malloc"); abort(); };
 np=npmax;
 reds = (point2 *)malloc(npmax*sizeof(point2));
 if (!reds) { perror("malloc"); abort(); };
 greens = (point2 *)malloc(npmax*sizeof(point2));
 if (!greens) { perror("malloc"); abort(); };
 }

int main(int argc,char *argv[]) {
 glutInit(&argc,argv);
 srand48(time(0));
 glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
 glutInitWindowSize(500,500);
 glutInitWindowPosition(0,0);
 glutCreateWindow("simple");
 glutDisplayFunc(display);
 init();
 glutMainLoop();
 }