public static double standardCDF(double x,boolean log_p){
  boolean i_tail=false;
  if (Double.isNaN(x)) {
    return Double.NaN;
  }
  double xden, xnum, temp, del, eps, xsq, y;
  int i;
  double p=x, cp=Double.NaN;
  boolean lower, upper;
  eps=DBL_EPSILON * 0.5;
  lower=!i_tail;
  upper=i_tail;
  y=Math.abs(x);
  if (y <= 0.67448975) {
    if (y > eps) {
      xsq=x * x;
      xnum=a[4] * xsq;
      xden=xsq;
      for (i=0; i < 3; i++) {
        xnum=(xnum + a[i]) * xsq;
        xden=(xden + b[i]) * xsq;
      }
    }
 else {
      xnum=xden=0.0;
    }
    temp=x * (xnum + a[3]) / (xden + b[3]);
    if (lower) {
      p=0.5 + temp;
    }
    if (upper) {
      cp=0.5 - temp;
    }
    if (log_p) {
      if (lower) {
        p=Math.log(p);
      }
      if (upper) {
        cp=Math.log(cp);
      }
    }
  }
 else   if (y <= M_SQRT_32) {
    xnum=c[8] * y;
    xden=y;
    for (i=0; i < 7; i++) {
      xnum=(xnum + c[i]) * y;
      xden=(xden + d[i]) * y;
    }
    temp=(xnum + c[7]) / (xden + d[7]);
    xsq=((int)(y * CUTOFF)) * 1.0 / CUTOFF;
    del=(y - xsq) * (y + xsq);
    if (log_p) {
      p=(-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp);
      if ((lower && x > 0.0) || (upper && x <= 0.0)) {
        cp=Math.log(1.0 - Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp);
      }
    }
 else {
      p=Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp;
      cp=1.0 - p;
    }
    if (x > 0.0) {
      temp=p;
      if (lower) {
        p=cp;
      }
      cp=temp;
    }
  }
 else   if (log_p || (lower && -37.5193 < x && x < 8.2924) || (upper && -8.2924 < x && x < 37.5193)) {
    xsq=1.0 / (x * x);
    xnum=p_[5] * xsq;
    xden=xsq;
    for (i=0; i < 4; i++) {
      xnum=(xnum + p_[i]) * xsq;
      xden=(xden + q[i]) * xsq;
    }
    temp=xsq * (xnum + p_[4]) / (xden + q[4]);
    temp=(M_1_SQRT_2PI - temp) / y;
    xsq=((int)(x * CUTOFF)) * 1.0 / CUTOFF;
    del=(x - xsq) * (x + xsq);
    if (log_p) {
      p=(-xsq * xsq * 0.5) + (-del * 0.5) + Math.log(temp);
      if ((lower && x > 0.0) || (upper && x <= 0.0)) {
        cp=Math.log(1.0 - Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp);
      }
    }
 else {
      p=Math.exp(-xsq * xsq * 0.5) * Math.exp(-del * 0.5) * temp;
      cp=1.0 - p;
    }
    if (x > 0.0) {
      temp=p;
      if (lower) {
        p=cp;
      }
      cp=temp;
    }
  }
 else {
    if (x > 0) {
      p=1.0;
      cp=0.0;
    }
 else {
      p=0.0;
      cp=1.0;
    }
  }
  return p;
}

public static double standardCDF ( double x , boolean log_p ) { boolean i_tail = false ; if ( Double . isNaN ( x ) ) { return Double . NaN ; } double xden , xnum , temp , del , eps , xsq , y ; int i ; double p = x , cp = Double . NaN ; boolean lower , upper ; eps = DBL_EPSILON * 0.5 ; lower = ! i_tail ; upper = i_tail ; y = Math . abs ( x ) ; if ( y <= 0.67448975 ) { if ( y > eps ) { xsq = x * x ; xnum = a [ 4 ] * xsq ; xden = xsq ; for ( i = 0 ; i < 3 ; i ++ ) { xnum = ( xnum + a [ i ] ) * xsq ; xden = ( xden + b [ i ] ) * xsq ; } } else { xnum = xden = 0.0 ; } temp = x * ( xnum + a [ 3 ] ) / ( xden + b [ 3 ] ) ; if ( lower ) { p = 0.5 + temp ; } if ( upper ) { cp = 0.5 - temp ; } if ( log_p ) { if ( lower ) { p = Math . log ( p ) ; } if ( upper ) { cp = Math . log ( cp ) ; } } } else if ( y <= M_SQRT_32 ) { xnum = c [ 8 ] * y ; xden = y ; for ( i = 0 ; i < 7 ; i ++ ) { xnum = ( xnum + c [ i ] ) * y ; xden = ( xden + d [ i ] ) * y ; } temp = ( xnum + c [ 7 ] ) / ( xden + d [ 7 ] ) ; xsq = ( ( int ) ( y * CUTOFF ) ) * 1.0 / CUTOFF ; del = ( y - xsq ) * ( y + xsq ) ; if ( log_p ) { p = ( - xsq * xsq * 0.5 ) + ( - del * 0.5 ) + Math . log ( temp ) ; if ( ( lower && x > 0.0 ) || ( upper && x <= 0.0 ) ) { cp = Math . log ( 1.0 - Math . exp ( - xsq * xsq * 0.5 ) * Math . exp ( - del * 0.5 ) * temp ) ; } } else { p = Math . exp ( - xsq * xsq * 0.5 ) * Math . exp ( - del * 0.5 ) * temp ; cp = 1.0 - p ; } if ( x > 0.0 ) { temp = p ; if ( lower ) { p = cp ; } cp = temp ; } } else if ( log_p || ( lower && - 37.5193 < x && x < 8.2924 ) || ( upper && - 8.2924 < x && x < 37.5193 ) ) { xsq = 1.0 / ( x * x ) ; xnum = p_ [ 5 ] * xsq ; xden = xsq ; for ( i = 0 ; i < 4 ; i ++ ) { xnum = ( xnum + p_ [ i ] ) * xsq ; xden = ( xden + q [ i ] ) * xsq ; } temp = xsq * ( xnum + p_ [ 4 ] ) / ( xden + q [ 4 ] ) ; temp = ( M_1_SQRT_2PI - temp ) / y ; xsq = ( ( int ) ( x * CUTOFF ) ) * 1.0 / CUTOFF ; del = ( x - xsq ) * ( x + xsq ) ; if ( log_p ) { p = ( - xsq * xsq * 0.5 ) + ( - del * 0.5 ) + Math . log ( temp ) ; if ( ( lower && x > 0.0 ) || ( upper && x <= 0.0 ) ) { cp = Math . log ( 1.0 - Math . exp ( - xsq * xsq * 0.5 ) * Math . exp ( - del * 0.5 ) * temp ) ; } } else { p = Math . exp ( - xsq * xsq * 0.5 ) * Math . exp ( - del * 0.5 ) * temp ; cp = 1.0 - p ; } if ( x > 0.0 ) { temp = p ; if ( lower ) { p = cp ; } cp = temp ; } } else { if ( x > 0 ) { p = 1.0 ; cp = 0.0 ; } else { p = 0.0 ; cp = 1.0 ; } } return p ; }
a more accurate and faster implementation of the cdf ( taken from function pnorm in the r statistical language ) this implementation has discrepancies depending on the programming language and system architecture in java , returned values become zero once z reaches - 37 . 5193 exactly on the machine tested in the other implementation , the returned value 0 at about z = - 8 in c , this 0 value is reached approximately z = - 37 . 51938 < p / > will later need to be optimised for beast
