public static double lbeta(double a, double b) {

        double corr, p, q;

        p = q = a;
        if (b < p) p = b;/* := min(a,b) */
        if (b > q) q = b;/* := max(a,b) */

        /*!* #ifdef IEEE_754 /*4!*/
        if (Double.isNaN(a) || Double.isNaN(b)) return a + b;
        /*!* #endif /*4!*/

        /* both arguments must be >= 0 */

        if (p < 0) {
            throw new java.lang.ArithmeticException("Math Error: DOMAIN");
            //    	return Double.NaN;
        }
        else if (p == 0) {
            return Double.POSITIVE_INFINITY;
        }
        /*!* #ifdef IEEE_754 /*4!*/
        else if (Double.isInfinite(q)) { return Double.NEGATIVE_INFINITY; }
        /*!* #endif /*4!*/

        if (p >= 10) {
            /* p and q are big. */
            corr = lgammacor(p) + lgammacor(q) - lgammacor(p + q);
            /*!* 	return log(q) * -0.5 + Constants.M_LN_SQRT_2PI + corr *!*/
            return java.lang.Math.log(q) * -0.5 + Constants.M_LN_SQRT_2PI + corr
            /*!* 		+ (p - 0.5) * log(p / (p + q)) + q * logrelerr(-p / (p + q)); *!*/
            + (p - 0.5) * java.lang.Math.log(p / (p + q)) + q * logrelerr(-p / (p + q));
        }
        else if (q >= 10) {
            /* p is small, but q is big. */
            corr = lgammacor(q) - lgammacor(p + q);
            /*!* 	return lgammafn(p) + corr + p - p * log(p + q) *!*/
            return lgammafn(p) + corr + p - p * java.lang.Math.log(p + q) + (q - 0.5) * logrelerr(-p / (p + q));
        }
        else
        /* p and q are small: p <= q > 10. */
        /*!* 	return log(gammafn(p) * (gammafn(q) / gammafn(p + q))); *!*/
        return java.lang.Math.log(gammafn(p) * (gammafn(q) / gammafn(p + q)));
    }