gt:             xnum = p_[5] * xsq;
context: /*
   * NormalDistribution.java
   *
   * Copyright (c) 2002-2011 Alexei Drummond, Andrew Rambaut and Marc Suchard
   *
   * This file is part of BEAST.
   * See the NOTICE file distributed with this work for additional
   * information regarding copyright ownership and licensing.
   *
   * BEAST is free software; you can redistribute it and/or modify
   * it under the terms of the GNU Lesser General Public License as
   * published by the Free Software Foundation; either version 2
   * of the License, or (at your option) any later version.
   *
   *  BEAST is distributed in the hope that it will be useful,
   *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   *  GNU Lesser General Public License for more details.
   *
   * You should have received a copy of the GNU Lesser General Public
   * License along with BEAST; if not, write to the
   * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
   * Boston, MA  02110-1301  USA
   */
  
  package dr.math.distributions;
  
  import dr.math.ErrorFunction;
  import dr.math.UnivariateFunction;
  
  /**
   * normal distribution (pdf, cdf, quantile)
   *
   * @author Korbinian Strimmer
   * @version $Id: NormalDistribution.java,v 1.7 2005/05/24 20:26:01 rambaut Exp $
   */
  public class NormalDistribution implements Distribution {
      //
      // Public stuff
      //
  
      /**
       * Constructor
       */
      public NormalDistribution(double mean, double sd) {
          this.m = mean;
          this.sd = sd;
      }
  
      public double getMean() {
          return m;
      }
  
      public void setMean(double value) {
          m = value;
      }
  
      public double getSD() {
          return sd;
      }
  
      public void setSD(double value) {
          sd = value;
      }
  
      public double pdf(double x) {
          return pdf(x, m, sd);
      }
  
      public double logPdf(double x) {
          return logPdf(x, m, sd);
      }
  
      public double cdf(double x) {
          return cdf(x, m, sd);
      }
  
      public double quantile(double y) {
          return quantile(y, m, sd);
      }
  
      public double mean() {
          return mean(m, sd);
      }
  
      public double variance() {
          return variance(m, sd);
      }
  
      public final UnivariateFunction getProbabilityDensityFunction() {
          return pdfFunction;
      }
  
      private final UnivariateFunction pdfFunction = new UnivariateFunction() {
          public final double evaluate(double x) {
              return pdf(x);
          }
  
          public final double getLowerBound() {
              return Double.NEGATIVE_INFINITY;
          }
  
          public final double getUpperBound() {
              return Double.POSITIVE_INFINITY;
          }
      };
  
      /**
       * probability density function
       *
       * @param x  argument
       * @param m  mean
       * @param sd standard deviation
       * @return pdf at x
       */
      public static double pdf(double x, double m, double sd) {
          double a = 1.0 / (Math.sqrt(2.0 * Math.PI) * sd);
          double b = -(x - m) * (x - m) / (2.0 * sd * sd);
  
          return a * Math.exp(b);
      }
  
      /**
       * the natural log of the probability density function of the distribution
       *
       * @param x  argument
       * @param m  mean
       * @param sd standard deviation
       * @return log pdf at x
       */
      public static double logPdf(double x, double m, double sd) {
          double a = 1.0 / (Math.sqrt(2.0 * Math.PI) * sd);
          double b = -(x - m) * (x - m) / (2.0 * sd * sd);
  
          return Math.log(a) + b;
      }
  
      /**
       * cumulative density function
       *
       * @param x  argument
       * @param m  mean
       * @param sd standard deviation
       * @return cdf at x
       */
      public static double cdf(double x, double m, double sd) {
          return cdf(x, m, sd, false);
  //        double a = (x - m) / (Math.sqrt(2.0) * sd);
  //
  //        return 0.5 * (1.0 + ErrorFunction.erf(a));
      }
  
      /**
       * quantiles (=inverse cumulative density function)
       *
       * @param z  argument
       * @param m  mean
       * @param sd standard deviation
       * @return icdf at z
       */
      public static double quantile(double z, double m, double sd) {
          return m + Math.sqrt(2.0) * sd * ErrorFunction.inverseErf(2.0 * z - 1.0);
      }
  
      /**
       * mean
       *
       * @param m  mean
       * @param sd standard deviation
       * @return mean
       */
      public static double mean(double m, double sd) {
          return m;
      }
  
      /**
       * variance
       *
       * @param m  mean
       * @param sd standard deviation
       * @return variance
       */
      public static double variance(double m, double sd) {
          return sd * sd;
      }
  
  
  
      /** 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
       *
       * Will later need to be optimised for BEAST
       *
       * @param x     argument
       * @param mu    mean
       * @param sigma standard deviation
       * @param log_p is p logged
       * @return cdf at x
       */
      public static double cdf(double x, double mu, double sigma, boolean log_p) {
  
          if(Double.isNaN(x) || Double.isNaN(mu) || Double.isNaN(sigma)) {
              return Double.NaN;
          }
          if(Double.isInfinite(x) && mu == x) { /* x-mu is NaN */
              return Double.NaN;
          }
          if (sigma <= 0) {
              if(sigma < 0) {
                  return Double.NaN;
              }
              return (x < mu) ? 0.0 : 1.0;
          }
          double p = (x - mu) / sigma;
          if(Double.isInfinite(p)) {
              return (x < mu) ? 0.0 : 1.0;
          }
          return standardCDF(p, log_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
       *
       * Will later need to be optimised for BEAST
       *
       * @param x     argument
       * @param log_p is p logged
       * @return cdf at x
       */
      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) { /* Normal.quantile(3/4, 1, 0) = 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) {
              /* Evaluate pnorm for 0.67448975 = Normal.quantile(3/4, 1, 0) < |x| <= sqrt(32) ~= 5.657 */
  
              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]);
  
              //do_del(y);
              //swap_tail;
              //#define do_del(X)							\
              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;
              }
              //#define swap_tail						\
              if (x > 0.0) {
                  temp = p;
                  if(lower) {
                      p = cp;
                  }
                  cp = temp;
              }
          }
          /* else	  |x| > sqrt(32) = 5.657 :
           * the next two case differentiations were really for lower=T, log=F
           * Particularly	 *not*	for  log_p !
           * Cody had (-37.5193 < x  &&  x < 8.2924) ; R originally had y < 50
           * Note that we do want symmetry(0), lower/upper -> hence use y
           */
          else if(log_p || (lower && -37.5193 < x  &&  x < 8.2924)
                  || (upper && -8.2924  < x  &&  x < 37.5193)) {
  
              /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
              xsq = 1.0 / (x * x);
  