Returns the probability of k of a specific number of tries n and probability p

@param k
@param p
@param n
@return
12168
public static double binomial(int k, double p, int n) {
        if(k<0 ||  p<0 || n<1) {
            throw new IllegalArgumentException("All the parameters must be positive and n larger than 1.");
        }
        
        k = Math.min(k, n); 
        
        /*
        //Slow and can't handle large numbers
        $probability=StatsUtilities::combination($n,$k)*pow($p,$k)*pow(1-$p,$n-$k);
        */

        //fast and can handle large numbers
        //Cdf(k)-Cdf(k-1)
        double probability = approxBinomialCdf(k,p,n); 
        if(k>0) {
            probability -= approxBinomialCdf(k-1,p,n);
        }
        
        return probability;
    }
111804
def EvalBinomialPmf(k, n, p):
    
    return scipy.stats.binom.pmf(k, n, p)
622485
static
	public double probability(ContinuousDistribution distribution, Number x){

		if(distribution instanceof GaussianDistribution){
			return probability((GaussianDistribution)distribution, x);
		} else

		if(distribution instanceof PoissonDistribution){
			return probability((PoissonDistribution)distribution, x);
		}

		throw new UnsupportedElementException(distribution);
	}
429648
public static double logChoose(int n, int k) {
        if (n < 0 || k < 0 || k > n) {
            throw new IllegalArgumentException(String.format("Invalid n = %d, k = %d", n, k));
        }

        return Math.logFactorial(n) - Math.logFactorial(k) - Math.logFactorial(n - k);
    }
158878
function randProbPrimeRounds(k,n) {
    var ans, i, divisible, B; 
    B=30000;  //B is largest prime to use in trial division
    ans=int2bigInt(0,k,0);
    
    //optimization: try larger and smaller B to find the best limit.
    
    if (primes.length==0)
      primes=findPrimes(30000);  //check for divisibility by primes <=30000

    if (rpprb.length!=ans.length)
      rpprb=dup(ans);

    for (;;) { //keep trying random values for ans until one appears to be prime
      //optimization: pick a random number times L=2*3*5*...*p, plus a 
      //   random element of the list of all numbers in [0,L) not divisible by any prime up to p.
      //   This can reduce the amount of random number generation.
      
      randBigInt_(ans,k,0); //ans = a random odd number to check
      ans[0] |= 1; 
      divisible=0;
    
      //check ans for divisibility by small primes up to B
      for (i=0; (i<primes.length) && (primes[i]<=B); i++)
        if (modInt(ans,primes[i])==0 && !equalsInt(ans,primes[i])) {
          divisible=1;
          break;
        }      
      
      //optimization: change millerRabin so the base can be bigger than the number being checked, then eliminate the while here.
      
      //do n rounds of Miller Rabin, with random bases less than ans
      for (i=0; i<n && !divisible; i++) {
        randBigInt_(rpprb,k,0);
        while(!greater(ans,rpprb)) //pick a random rpprb that's < ans
          randBigInt_(rpprb,k,0);
        if (!millerRabin(ans,rpprb))
          divisible=1;
      }
      
      if(!divisible)
        return ans;
    }  
  }
275749
public static int nextInt(int n)
  {
    Random random = getRandom();

    int value = random.nextInt(n);

    if (! _isTest)
      _freeRandomList.free(random);

    return value;
  }
1151715
@Reference(title = "Fast and accurate computation of binomial probabilities", //
      authors = "C. Loader", booktitle = "", //
      url = "http://projects.scipy.org/scipy/raw-attachment/ticket/620/loader2000Fast.pdf", //
      bibkey = "web/Loader00")
  private static double stirlingError(int n) {
    // Try to use a table value:
    if(n < 16) {
      return STIRLING_EXACT_ERROR[n << 1];
    }
    final double nn = n * n;
    // Use the appropriate number of terms
    if(n > 500) {
      return (S0 - S1 / nn) / n;
    }
    if(n > 80) {
      return ((S0 - (S1 - S2 / nn)) / nn) / n;
    }
    if(n > 35) {
      return ((S0 - (S1 - (S2 - S3 / nn) / nn) / nn) / n);
    }
    return ((S0 - (S1 - (S2 - (S3 - S4 / nn) / nn) / nn) / nn) / n);
  }
511248
def miller_rabin(n, k, rnd=default_pseudo_random):
    '''
       
   '''
    s = 0
    d = n-1
    # Find nearest power of 2
    s = primitives.integer_bit_size(n)
    # Find greatest factor which is a power of 2
    s = fractions.gcd(2**s, n-1)
    d = (n-1) // s
    s = primitives.integer_bit_size(s) - 1
    while k:
        k = k - 1
        a = rnd.randint(2, n-2)
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for r in range(1, s-1):
            x = pow(x, 2, n)
            if x == 1:
                return False
            if x == n - 1:
                break
        else:
            return False
    return True
403206
@Override
  public int compareTo(Object o) {
    SortedClassProbability other = (SortedClassProbability) o;
    if (this.probability < other.probability) {
      return -1;
    }
    else if (this.probability > other.probability) {
      return 1;
    }
    else {
      return 0;
    }
  }
1213400
@Override protected double[] score0(double data[/*ncols*/], double preds[/*nclasses+1*/], double offset, int ntrees) {
    super.score0(data, preds, offset, ntrees);    // These are f_k(x) in Algorithm 10.4
    return score0Probabilities(preds, offset);
  }
460403
private final int m() {
		int n = 0;
		int i = k0;
		while (true) {
			if (i > j)
				return n;
			if (!cons(i))
				break;
			i++;
		}
		i++;
		while (true) {
			while (true) {
				if (i > j)
					return n;
				if (cons(i))
					break;
				i++;
			}
			i++;
			n++;
			while (true) {
				if (i > j)
					return n;
				if (!cons(i))
					break;
				i++;
			}
			i++;
		}
	}
517035
public static double bernoulliCdf(int k, double p) {
        if(p<0) {
            throw new IllegalArgumentException("The probability p can't be negative.");
        }
        
        double probabilitySum=0.0;
        if(k<0) {
            
        }
        else if(k<1) { //aka k==0
            probabilitySum=(1-p);
        }
        else { //k>=1 aka k==1
            probabilitySum=1.0;
        }
        
        return probabilitySum;
    }
1082698
public static double hypergeometricCdf(int k, int n, int Kp, int Np) {
        if(k<0 || n<0 || Kp<0 || Np<0) {
            throw new IllegalArgumentException("All the parameters must be positive.");
        }
        Kp = Math.max(k, Kp);
        Np = Math.max(n, Np);
        
        /*
        //slow!
        $probabilitySum=0;
        for($i=0;$i<=$k;++$i) {
            $probabilitySum+=self::hypergeometric($i,$n,$Kp,$Np);
        }
        */

        //fast and can handle large numbers
        //Cdf(k)-Cdf(k-1)
        double probabilitySum = approxHypergeometricCdf(k,n,Kp,Np);
        
        return probabilitySum;
    }
1159322
public static long count(int n, int m) {
		if(n == m) {
			return NumberUtil.factorial(n);
		}
		return (n > m) ? NumberUtil.factorial(n, n - m) : 0;
	}
1143523
public boolean nextBoolean (final double probability)
        {
        if (probability < 0.0 || probability > 1.0)
            throw new IllegalArgumentException ("probability must be between 0.0 and 1.0 inclusive.");
        if (probability==0.0) return false;             // fix half-open issues
        else if (probability==1.0) return true; // fix half-open issues
        return nextDouble() < probability; 
        }
1175866
public static function comb(int $n, int $k = 2) : int
    {
        return $k === 0 ? 1 : (int) (($n * self::comb($n - 1, $k - 1)) / $k);
    }
312812
def recomb_probability(cM, method="kosambi"):
    
    assert method in ("kosambi", "haldane")
    d = cM / 100.
    if method == "kosambi":
        e4d = exp(4 * d)
        return (e4d - 1) / (e4d + 1) / 2
    elif method == "haldane":
        return (1 - exp(-2 * d)) / 2