#include #include #include #include //--- 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) s=s1=s2=0; t1=exp(-b); t2=exp(-c); for(i=1;i<=k;i++){ s1+=t1; t1*=b/i; // s1=P(Xb=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=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*rp)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.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;j2*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 }