Continued fraction expansion #1 for incomplete beta integral.

@param a Value.
@param b Value.
@param x Value.
@return Result.
771664
public static double Incbcf(double a, double b, double x) {
        double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
        double k1, k2, k3, k4, k5, k6, k7, k8;
        double r, t, ans, thresh;
        int n;
        double big = 4.503599627370496e15;
        double biginv = 2.22044604925031308085e-16;

        k1 = a;
        k2 = a + b;
        k3 = a;
        k4 = a + 1.0;
        k5 = 1.0;
        k6 = b - 1.0;
        k7 = k4;
        k8 = a + 2.0;

        pkm2 = 0.0;
        qkm2 = 1.0;
        pkm1 = 1.0;
        qkm1 = 1.0;
        ans = 1.0;
        r = 1.0;
        n = 0;
        thresh = 3.0 * Constants.DoubleEpsilon;

        do {
            xk = -(x * k1 * k2) / (k3 * k4);
            pk = pkm1 + pkm2 * xk;
            qk = qkm1 + qkm2 * xk;
            pkm2 = pkm1;
            pkm1 = pk;
            qkm2 = qkm1;
            qkm1 = qk;

            xk = (x * k5 * k6) / (k7 * k8);
            pk = pkm1 + pkm2 * xk;
            qk = qkm1 + qkm2 * xk;
            pkm2 = pkm1;
            pkm1 = pk;
            qkm2 = qkm1;
            qkm1 = qk;

            if (qk != 0) r = pk / qk;
            if (r != 0) {
                t = Math.abs((ans - r) / r);
                ans = r;
            } else
                t = 1.0;

            if (t < thresh) return ans;

            k1 += 1.0;
            k2 += 1.0;
            k3 += 2.0;
            k4 += 2.0;
            k5 += 1.0;
            k6 -= 1.0;
            k7 += 2.0;
            k8 += 2.0;

            if ((Math.abs(qk) + Math.abs(pk)) > big) {
                pkm2 *= biginv;
                pkm1 *= biginv;
                qkm2 *= biginv;
                qkm1 *= biginv;
            }
            if ((Math.abs(qk) < biginv) || (Math.abs(pk) < biginv)) {
                pkm2 *= big;
                pkm1 *= big;
                qkm2 *= big;
                qkm1 *= big;
            }
        } while (++n < 300);

        return ans;
    }
1275231
public static double ComplementedIncomplete(double a, double x) {
        final double big = 4.503599627370496e15;
        final double biginv = 2.22044604925031308085e-16;
        double ans, ax, c, yc, r, t, y, z;
        double pk, pkm1, pkm2, qk, qkm1, qkm2;

        if (x <= 0 || a <= 0) return 1.0;

        if (x < 1.0 || x < a) return 1.0 - Incomplete(a, x);

        ax = a * Math.log(x) - x - Log(a);
        if (ax < -Constants.LogMax) return 0.0;

        ax = Math.exp(ax);

        // continued fraction
        y = 1.0 - a;
        z = x + y + 1.0;
        c = 0.0;
        pkm2 = 1.0;
        qkm2 = x;
        pkm1 = x + 1.0;
        qkm1 = z * x;
        ans = pkm1 / qkm1;

        do {
            c += 1.0;
            y += 1.0;
            z += 2.0;
            yc = y * c;
            pk = pkm1 * z - pkm2 * yc;
            qk = qkm1 * z - qkm2 * yc;
            if (qk != 0) {
                r = pk / qk;
                t = Math.abs((ans - r) / r);
                ans = r;
            } else
                t = 1.0;

            pkm2 = pkm1;
            pkm1 = pk;
            qkm2 = qkm1;
            qkm1 = qk;
            if (Math.abs(pk) > big) {
                pkm2 *= biginv;
                pkm1 *= biginv;
                qkm2 *= biginv;
                qkm1 *= biginv;
            }
        } while (t > Constants.DoubleEpsilon);

        return ans * ax;
    }
90624
public static long subtractExact(long a, long b) {
        long result = a - b;
        checkNoOverflow((a ^ b) >= 0 | (a ^ result) >= 0);
        return result;
    }
1261182
public static double regularizedIncompleteGamma(double s, double x) {
        if (s < 0.0) {
            throw new IllegalArgumentException("Invalid s: " + s);
        }

        if (x < 0.0) {
            throw new IllegalArgumentException("Invalid x: " + x);
        }

        double igf = 0.0;

        if (x < s + 1.0) {
            // Series representation
            igf = regularizedIncompleteGammaSeries(s, x);
        } else {
            // Continued fraction representation
            igf = regularizedIncompleteGammaFraction(s, x);
        }

        return igf;
    }
100686
public static long safeMultiply(long a, int b) {
        switch (b) {
            case -1:
                if (a == Long.MIN_VALUE) {
                    throw new ArithmeticException("Multiplication overflows a long: " + a + " * " + b);
                }
                return -a;
            case 0:
                return 0L;
            case 1:
                return a;
        }
        long total = a * b;
        if (total / b != a) {
            throw new ArithmeticException("Multiplication overflows a long: " + a + " * " + b);
        }
        return total;
    }
593370
public static double calculateScaledResidual(double[][] A, double[] x, double[] b){
		DoubleMatrix2D AMatrix = DoubleFactory2D.dense.make(A); 
		DoubleMatrix1D xVector = DoubleFactory1D.dense.make(x);
		DoubleMatrix1D bVector = DoubleFactory1D.dense.make(b);
		return calculateScaledResidual(AMatrix, xVector, bVector);
	}
257056
public static long saturatedMultiply(long a, long b) {
        // see checkedMultiply for explanation
        int leadingZeros = Long.numberOfLeadingZeros(a) + Long.numberOfLeadingZeros(~a) + Long.numberOfLeadingZeros(b) + Long.numberOfLeadingZeros(~b);
        if (leadingZeros > Long.SIZE + 1) {
            return a * b;
        }
        // the return value if we will overflow (which we calculate by overflowing a long :) )
        long limit = Long.MAX_VALUE + ((a ^ b) >>> (Long.SIZE - 1));
        if (leadingZeros < Long.SIZE | (a < 0 & b == Long.MIN_VALUE)) {
            // overflow
            return limit;
        }
        long result = a * b;
        if (a == 0 || result / a == b) {
            return result;
        }
        return limit;
    }
1189297
@Pure
	@Inline(value="($1 < $2)", constantExpression=true)
	public static boolean operator_lessThan(long a, double b) {
		return a < b;
	}
430510
public static double uniformCdf(double x, double a, double b) {
        if(a>=b) {
            throw new IllegalArgumentException("The a must be smaller than b.");
        }
        
        double probabilitySum;
        if(x<a) {
            probabilitySum=0.0;
        }
        else if(x<b) {
            probabilitySum=(x-a)/(b-a);
        }
        else {
            probabilitySum=1;
        }

        return probabilitySum;
    }
956384
public static int multiplyExact(int a, int b) {
        long result = (long) a * b;
        checkNoOverflow(result == (int) result);
        return (int) result;
    }
563134
@Pure
	@Inline(value="($1 ^ $2)", constantExpression=true)
	public static int bitwiseXor(int a, int b) {
		return a ^ b;
	}
390881
@Pure
	@Inline(value="($1 * $2)", constantExpression=true)
	public static long operator_multiply(short a, long b) {
		return a * b;
	}
797239
public static function continuedFraction($a, $b, $x)
    {
        $maxit = 100;
        $eps = 3e-16;
        $fpmin = 1e-30;
        
        $qab = $a + $b;
        $qap = $a + 1;
        $qam = $a - 1;

        $c = 1.0;
        $d = 1.0 - $qab*$x/$qap;

        if (abs($d)<$fpmin) {
            $d = $fpmin;
        }

        $d = 1.0/$d;

        $h = $d;

        for ($m = 1; $m < $maxit; $m++) {
            $m2 = 2*$m;
            $aa = $m*($b-$m)*$x/(($qam + $m2)*($a + $m2));
            $d = 1.0 + $aa*$d;

            if (abs($d)<$fpmin) {
                $d = $fpmin;
            }

            $c = 1.0 + $aa/$c;

            if (abs($c)<$fpmin) {
                $c = $fpmin;
            }

            $d = 1.0/$d;
            $h = $h*$d*$c;
            $aa = -($a + $m)*($qab + $m)*$x/(($a+$m2)*($qap+$m2));
            $d = 1.0 + $aa*$d;

            if (abs($d)<$fpmin) {
                $d = $fpmin;
            }

            $c = 1.0 + $aa/$c;

            if (abs($c)<$fpmin) {
                $c = $fpmin;
            }

            $d = 1.0/$d;
            $del = $d*$c;
            $h = $h*$del;

            if (abs($del-1.0)< $eps) {
                break;
            }
        }
        return $h;
    }
156394
public static double logpdf(double x, double mu, double beta) {
    if(x == Double.NEGATIVE_INFINITY) {
      return Double.NEGATIVE_INFINITY;
    }
    final double z = (x - mu) / beta;
    return -z - FastMath.exp(-z) - FastMath.log(beta);
  }
832474
public static BigInteger lcm(BigInteger a, BigInteger b) {
		return lcm(a, b, Calculator.BIG_INTEGER_CALCULATOR);
	}
931516
public static int safeSubtract(int a, int b) {
        int result = a - b;
        // check for a change of sign in the result when the inputs have the different signs
        if ((a ^ result) < 0 && (a ^ b) < 0) {
            throw new ArithmeticException("Subtraction overflows an int: " + a + " - " + b);
        }
        return result;
    }
1080422
@javax.annotation.Nonnull
  public static double[] sum(@javax.annotation.Nonnull final double[] a, final double b) {
    return com.simiacryptus.util.ArrayUtil.op(a, (x) -> x + b);
  }