/* Compute probabilities for passing random parameter selection phases in
   fractal.c
                                                                          */
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<time.h>

#define SAMPLE_COUNT  10000000

/* Compute cross product of 2 vectors. */
static double CrossProduct(double ax, double ay, double bx, double by)
{
   return ax * by - ay * bx;
}

/* Compute area of parallelogram from 3 corner vertices. */
static double GetParallelogramArea(double ax, double ay,
                                   double bx, double by,
                                   double cx, double cy)
{
   return fabs(CrossProduct(bx - ax, by - ay, cx - ax, cy - ay));
}

/* Get the rate of random parallelograms that will fall within the
   accepted range.                                                 */
static double GetAcceptRate(double scale, double min_area, double max_area)
{
   int i, j, samples, count;
   double vertex[3][2], area;

   count = 0;
   for(samples = 0; samples < SAMPLE_COUNT; samples++)
   {
      for(i = 0; i < 3; i++)
      {
         for(j = 0; j < 2; j++)
            vertex[i][j] = ((double)rand() / RAND_MAX) * scale;
      }
      area = GetParallelogramArea(vertex[0][0], vertex[0][1],
                                  vertex[1][0], vertex[1][1],
                                  vertex[2][0], vertex[2][1]);
      if( area >= min_area && area <= max_area )
         count++;
   }
   return (double)count / SAMPLE_COUNT;
}

/* Run test and output results. */
static void TestAcceptRate(double scale, double min_area, double max_area)
{
   printf("Max rand = %g, accept area = %g .. %g, accept rate = %.5f\n",
          scale, min_area, max_area, GetAcceptRate(scale, min_area, max_area));
}

/* Get a distribution of random parallelogram areas. */
int main(int argc, char **argv)
{
   srand(time(NULL));

   TestAcceptRate(0.8, 0.4, 1.0);
   TestAcceptRate(1.0, 0.3, 0.6);
   TestAcceptRate(1.0, 0.3, 0.7);
   TestAcceptRate(1.0, 0.4, 0.7);
   TestAcceptRate(1.0, 0.5, 0.8);
   return 0;
}