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

//--- The following parameters may be changed --------------------------

#define	X	'U'		// Case X = 'U'/'O'/'I'/'B'
#define	KMIN	9		// Minimum value of k
#define	KMAX	15		// Maximum value of k 
#define ACC	100		// Increase to get better bound

//----------------------------------------------------------------------

#define PI 3.14159265358979323846 // pi

// Calculate area L(a,b,t) of lune inside disk radius a,
// outside disk radius b, angle of at the endpoints of the lune = t,
// Requires |t|<=PI, a,b>=0. For t<0, L(a,b,t) defined as L(b,a,-t).
//
double lune(double a,double b,double t){
 double c,p,q;
 if(t<0)return -lune(b,a,-t);	// L(a,b,t)=-L(b,a,-t) for t<0
 if(a<0||b<0||t>PI){printf("Error in argument to lune()\n");exit(1);}
 if(b==0)return PI*a*a;		// b=0 causes undefined value of q
 if(t==0)return 0;		// t=0 may cause problems when a=b
 c=sqrt(a*a+b*b-2*a*b*cos(t));	// distance between centers of disks
 q=acos((b*b+c*c-a*a)/(2*b*c));	// half angle of lune from center of D_b 
 p=q+t;				// half angle of lune from center of D_a
 return (p-sin(2*p)/2)*a*a-(q-sin(2*q)/2)*b*b;
}

// Calculate probability P(Po(a)>=k)
//
double pois_tail(int k, double a){
 double s,t;
 int i;
 s=0; t=1;
 for(i=0;i<k;t*=a/(++i))s+=t;
 return 1-exp(-a)*s;
}

// Calculate p_X(u)e^{|B|} given areas
// a=|A\cap C|, b=|A\setminus C|, c=|C\setminus (A\cup B)|.
//
double split(int k,double a,double b,double c){
 double s1,t1,s2,t2,s,t;
 int i;
// Write Xa, Xb, Xc for independent Po(a), Po(b), Po(c) RVs
 switch(X){
 case 'U':			// P(Xa+Xb and Xa+Xc >=k)
  s=s1=s2=0; t1=exp(-b); t2=exp(-c);
  for(i=1;i<=k;i++){
   s1+=t1; t1*=b/i;		// s1=P(Xb<i)
   s2+=t2; t2*=c/i;		// s2=P(Xc<i)
   s=s1+s2-s1*s2+a*s/(k-i+1);	// P(Xb or Xc<k-Xa & Xa>=k-i)/P(Xa=k-i)
  }
  return 1-exp(-a)*s;
 case 'B':			// P(Xa+Xb or Xa+Xc >=k)
  s=s1=s2=0; t1=exp(-b); t2=exp(-c);
  for(i=1;i<=k;i++){
   s1+=t1; t1*=b/i;		// s1=P(Xb<i)
   s2+=t2; t2*=c/i;		// s2=P(Xc<i)
   s=s1*s2+a*s/(k-i+1);		// P(Xb & Xc<k-Xa & Xa>=k-i)/P(Xa=k-i)
  }
  return 1-exp(-a)*s;
 case 'O':			// P(Xa+Xb >=k)
  return pois_tail(k,a+b);
 case 'I':			// P(Xa+Xc >=k)
  return pois_tail(k,a+c);
 default:
  printf("Option %c not allowed in split()\n",X); exit(1);
 }
}

// Calculate p_X(u) with |u-v|=a, angle Ovu=t
//
double prob(int k,double r,double a,double t){
 double p,A,C_AB,AC,B;
 if(t<0)t=-t;			// wlog t>=0 by symmetry
 if(2*r<a)return 1;		// u not covered by D_v
 p=acos(a/(2*r));		// critical angle \phi
 if(t>p)return 1;		// u not covered by D_v
 A=lune(a,r,p);			// |A|
 C_AB=(PI/3+sqrt(.75))*a*a;	// |C\(AuB)|=|(AuB)\C|
 if(t<=p-PI/3 && t>PI/3-p)	// Calculate |A\cap C| (see Appendix A)
  AC=lune(a,r,t);
 else if(t>p-PI/3 && t>PI/3-p)
  AC=lune(a,r,p-PI/3)+(t-p+PI/3)*a*a/2; // may use lune() with angle<0
 else
  AC=A-C_AB+lune(r,a,t);
 B=PI*a*a-A;			// |B|
 return split(k,AC,A-AC,C_AB)*exp(-B);
}

// P(failure for given k,r,s) minimized over s
// Estimate integrals using Simplson's rule,
// fast and accurate but not rigorous.
// This is used to find parameters r,s.
// Finding s is easy, we just stop integrating
// when we hit a local minimum (<1).
// n=accuracy parameter.
//
double fail_est(int k,double r,double *s,int n){
 double f,a,p,t,m,v;
 int i,j;
 n*=2;				// double n, so even
 f=0; *s=0; m=1;
 for(i=1;;i++){			// integrate over a
  a=(double)i/n;
  if(a>2*r)return m;		// return m anyway (s=2r)
  p=acos(a/(2*r))-1e-15;	// -1e-15 avoids errors
  for(t=0,j=0;j<=n;j++){	// integrate over angle
   t+=prob(k,r,a,j*p/n)*((j&1)*2+(j>0)+(j<n))/3.0; // Simpson
  }
  t*=2*a*(p/n);			// 2 to include negative angles
  if(i&1){
   f+=4*t/(3*n);		// Simpson
  }
  else {
   v=f+t/(3*n);			// finish integral
   v+=exp(lune(r,a,p)-PI*r*r);
   v=v*4*r*(r+a)+exp(-PI*r*r);	// estimate when s=a
   if(v<=m){m=v;*s=a;}
   else if(m<1)return m;
   f+=2*t/(3*n);		// Simpson
  }
 }
}

// min_{r,s} P(failure for given k,r,s).
// Minimizes fail_est() over r in [0,10] using golden
// section search. Fast and accurate, but neither the
// calculation of the failure probability, nor the
// method of minimization, is rigorously justified.
// This does not matter, since min_est() is just used
// to find suitable parameters r,s to feed into fail_rig(),
// and any value of r and s will then give a rigorous bound.
// (The worst that can happen is that the bound isn't very good.)
// n=accuracy parameter.
//
double min_est(int k,double *r,double *s, int n){
 double t,h,u,v;
 t=(1+sqrt(5))/2;		// golden ratio
 h=10;				// interval length
 *r=h/t;			// current point
 v=fail_est(k,*r,s,n);		// current value
 while(fabs(h)*n>.5){
  u=fail_est(k,*r-h/(t*t*t),s,n); // test new point
  if(u<=v){v=u;*r-=h/(t*t*t); h/=t;} // change interval
  else h=-h/t;
 }
 return v;			// return best found
}

// Rigorous upper bound for average p_X(u) value
// over all angles theta, for fixed value of a.
// Uses monotonicity in angle.
// Routine used in fail_rig().
// n=accuracy parameter
//
double tint_rig(int k,double r,double a,int n){
 double f,p,t;
 int j,m;
 p=acos(a/(2*r));		// angle limit
 f=prob(k,r,a,p);		// max of p_X(u)
 if(X=='O')return f;		// p_O(u) independent of angle!
 m=ceil((f-prob(k,r,a,0.))*n*200); // number of subdivisions
 if(m<1)m=1;
 for(j=1;j<m;j++)f+=prob(k,r,a,j*p/m);
 return f/m;			// upper Riemann sum
}

// P(failure for given k,r,s)
// Bounds probability rigorously
// Uses adaptive approach, changing step size.
//
double fail_rig(int k,double r,double s,int n){
 double f,a,am,b,bm,p,pm,t;
 int i,j,m;
 if(s>2*r)s=2*r;		// no point in s>2r
 f=exp(lune(r,s,acos(s/(2*r)))-PI*r*r);
 a=s;				// integrate from a=s to a=0
// we will integrate by bounding on intervals [am,a] of the radius
 while(a>0){
  t=tint_rig(k,r,a,n);		// average p_X(u) over larger radius
  am=a-1/(n*500*(t>.001?t:.001)); // choose smaller radius
  if(am<.0001)am=0;		// but not below 0
// also numerical problems close to a=0, so don't allow am too small
  p=acos(a/(2*r));		// max angle at radius a
  pm=acos(am/(2*r));		// max angle at radius am
  b=PI*a*a-lune(a,r,p);		// |B| at radius a
  bm=PI*am*am-lune(am,r,pm);	// |B| at radius am
  f+=exp(b-bm)*t*p*(a*a-am*am);	// pessimistic integral from am to a
// The p*(a*a-am*am) is area of annulus segment from angle
// -p to p, radius am to a. The e^{b-bm} adjusts to max e^{-|B|}.
// But p changes from am to a! so we have to deal with small
// triangular regions with angle between p and pm (or -p to -pm) and
// radius between am and a.
// But we always have p_X(u)\le (1 or 2)P(Po(|B|)=0)P(Po(\pi a^2)>=k)
// since B is empty and one (of two possible in case 'B') region
// of area \pi a^2 must have at least k points.
  t=(1+(X=='B'))*pois_tail(k,PI*a*a)*exp(-bm);
  f+=t*(pm-p)*(a*a-am*am);	// very crude bound for small triangles
  a=am;				// new value of max radius
 }
 return f*4*r*(r+s)+exp(-PI*r*r);
}

// Main routine. Prints all estimates and upper bounds.
//
main(){
 time_t start,end;		// just to time program, may be removed
 double r,s,t;
 int k;
 time(&start);			// start time
 for(k=KMIN;k<=KMAX;k++){
  printf("X=%c k=%d Accuracy=%d\n",X,k,ACC);
  t=min_est(k,&r,&s,ACC);
  printf("Estimate    = %f at r=%f s=%f\n",t,r,s);
// Use ACC*ACC for rigorous bound since estimate has O(1/n^2) error
// while rigorous bound has O(1/n) error (with much worse constant).
// 6e-7 adjusts for rounding errors in calculation (1e-7)
// and printing to 6 decimal places (5e-7)
  t=fail_rig(k,r,s,ACC*ACC)+6e-7;
  printf("Upper bound = %f\n\n",t);
  fflush(NULL);			// ensure data is displayed immediately
 }
 time(&end);			// end time
 printf("%d seconds\n",(int)(end-start)); // time taken
}