lib/bivnor.c

   1 #include "xmath.h"
   2
   3 #define twopi 6.283195307179587
   4 #define con (twopi / 2.0) * 10.0e-10
   5
   6 double bivnor(ah, ak, r)
   7      double ah, ak, r;
   8 {
   9   /*
  10     based on alg 4628 comm. acm oct 73
  11     gives the probability that a bivariate normal exceeds (ah,ak).
  12     gh and gk are .5 times the right tail areas of ah, ak under a n(0,1)
  13
  14     Tranlated from FORTRAN to ratfor using struct; from ratfor to C by hand.
  15   */
  16   double a2, ap, b, cn, conex, ex, g2, gh, gk, gw, h2, h4, rr, s1, s2,
  17     sgn, sn, sp, sqr, t, temp, w2, wh, wk;
  18   int is;
  19
  20   temp = -ah;
  21   normbase(&temp, &gh);
  22   gh = gh / 2.0;
  23   temp = -ak;
  24   normbase(&temp, &gk);
  25   gk = gk / 2.0;
  26
  27   b = 0;
  28
  29   if (r==0)
  30     b = 4*gh*gk;
  31   else {
  32     rr = 1-r*r;
  33     if (rr<0)
  34       return(-1.0);
  35     if (rr!=0) {
  36       sqr = sqrt(rr);
  37       if (ah!=0) {
  38         b = gh;
  39         if (ah*ak<0)
  40           b = b-.5;
  41         else if (ah*ak==0)
  42           goto label10;
  43       }
  44       else if (ak==0) {
  45         b = atan(r/sqr)/twopi+.25;
  46         goto label50;
  47       }
  48       b = b+gk;
  49       if (ah==0)
  50         goto label20;
  51     label10:
  52       wh = -ah;
  53       wk = (ak/ah-r)/sqr;
  54       gw = 2*gh;
  55       is = -1;
  56       goto label30;
  57     label20:
  58       do {
  59         wh = -ak;
  60         wk = (ah/ak-r)/sqr;
  61         gw = 2*gk;
  62         is = 1;
  63       label30:
  64         sgn = -1;
  65         t = 0;
  66         if (wk!=0) {
  67           if (fabs(wk)>=1)
  68             if (fabs(wk)==1) {
  69               t = wk*gw*(1-gw)/2;
  70               goto label40;
  71             }
  72             else {
  73               sgn = -sgn;
  74               wh = wh*wk;
  75               normbase(&wh, &g2);
  76               wk = 1/wk;
  77               if (wk<0)
  78                 b = b+.5;
  79               b = b-(gw+g2)/2+gw*g2;
  80             }
  81           h2 = wh*wh;
  82           a2 = wk*wk;
  83           h4 = h2*.5;
  84           ex = 0;
  85           if (h4<150.0)
  86             ex = exp(-h4);
  87           w2 = h4*ex;
  88           ap = 1;
  89           s2 = ap-ex;
  90           sp = ap;
  91           s1 = 0;
  92           sn = s1;
  93           conex = fabs(con/wk);
  94           do {
  95             cn = ap*s2/(sn+sp);
  96             s1 = s1+cn;
  97             if (fabs(cn)<=conex)
  98               break;
  99             sn = sp;
 100             sp = sp+1;
 101             s2 = s2-w2;
 102             w2 = w2*h4/sp;
 103             ap = -ap*a2;
 104           } while (1);
 105           t = (atan(wk)-wk*s1)/twopi;
 106         label40:
 107           b = b+sgn*t;
 108         }
 109         if (is>=0)
 110           break;
 111       } while(ak!=0);
 112     }
 113     else if (r>=0)
 114       if (ah>=ak)
 115         b = 2*gh;
 116       else
 117         b = 2*gk;
 118     else if (ah+ak<0)
 119       b = 2*(gh+gk)-1;
 120   }
 121  label50:
 122   if (b<0)
 123     b = 0;
 124   if (b>1)
 125     b = 1;
 126
 127   return(b);
 128 }