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

double x,y,z,s,alpha,beta,h,hh,P,Pmax;

double Bin(double N, double k, double hh);


/*    We begin by setting two constants, hh<h.

      The three loops break the region given in Claim 9 into
      cubes of side length hh.  For each point (x,y,z) given by 
      the loops, we bound g(x,y,z) over the cube 
      
             C(x,y,z) = [x,x+h] X [y,y+h] X [z,z+h].  

      Note that cubes over which we bound g(x,y,z) are slightly
      expanded versions of the cubes in the partition given by 
      the loops.  So, we bound g(x,y,z) over an overalpping 
      collection of cubes.  This averts the possibility that some 
      small part of the region is not checked due to round-off 
      error.

      In order to bound g(x,y,z) we view the function as
      a product of eight terms.  We take the product of absolute 
      upper bounds over C(x,y,z) of each of the eight terms.
    
      We output the upper bound produced over a cube C(x,y,z)
      only when it exceeds the current maximum upper bound.  

*/


int main(void)
{
  hh = 0.00125;
  h  = 0.00126;
  Pmax=.1;

  for(x=.015;x<=0.415;x=x+hh) 
    { 
      printf("x=%8f\n",x);

      for(y=0;y<=1-2*x;y=y+hh)
	  { 

	    s= 1-x-y;
	    if( s> 0.415) s =0.415;

	    for(z=0;z<=3*y;z=z+hh)
	      {

		alpha = 1 + s - 3*x - y - z;
		beta = 2 - s - 2*y + z;

	        P = Bin(1,x,h);

		P = P*Bin(1-x,y,h);

		P = P*Bin(3*(1-x-y),alpha - 7*h,7*h);

		P = P* Bin(3*y+3*h, z,h);
                                                        
                P = P*pow(x + h, beta - 2*h);

                P = P*pow(1 - x, alpha - 5*h);

                P = P*pow(1 - x - y, z);

		P = P*pow(x + y + 2*h, 3*y - z -h);
		   
                if(P>Pmax)
                   {
		    Pmax=P;
                    printf("Pmax=%8f x=%8f y=%8f z=%8f \n",P,x,y,z);
       		   }
	      }
	  }
    }
}

/*   
      The procedure Bin provides an upper bound on the function
     
                        a^a/( b^b (a-b)^{a-b})

      where 0<b<a.  The variable a is taken to be fixed and b lies 
      in the range  (b,b+eps).  Of course, in some cases a is not
      fixed over the given cube, here we simply call Bin with the
      largest possible value of a over the cube in question.

*/


double Bin(double a, double b, double eps)
{
  if( b + eps < a/2) return pow(a,a)/(pow(b+eps,b+eps)*pow(a-b-eps,a-b-eps));
  if( b > a/2) return pow(a,a)/(pow(b,b)*pow(a-b,a-b));
  return pow(2,a);
}