How to simplify or calculate a formula with very big factorials
I'm facing a practical problem where I've calculated a formula that, with the help of some programming, can bring me to my final answer. However, the numbers involved are so big that it takes ages to compute, and I think it might not be necessary. I have the following formula,
$\sum\limits^{k}_{i=m}(N-i)^{k-i}(\frac{1}{N})^k\frac{k!}{(k-i)!i!} \leq a$
from which I need to calculate N. The rest of the values are constants, but in the range of the 100,000s. The factorial there is giving me a headache, since the values involved are too large; what simplifications could I make that will loosen the bounds slightly and thereby simplify the calculation? Are there any standard tricks? Or perhaps a way to calculate this in matlab / octave?
You need Stirling's approximation. It is very accurate for large factorials.
You don't need to compute the individual factorials in order to compute $k!/(k-i)!i!$, since that's the binomial coefficient $\binom{k}{i}$. A simple algorithm for computing binomial coefficients can be found on Wikipedia. A more sophisticated algorithm is due to Goetgheluck (JSTOR); implementations can be found here and here.
Of course, with numbers of the size that you have, this might still not be feasible, and in this case I also recommend Stirling's formula.
Here is a good writeup about implementing Stirling's approximation, along with a reference implementation:
http://threebrothers.org/brendan/blog/stirlings-approximation-formula-clojure/
import java.util.HashMap;
/**
* Convenience methods for deterministic combinations.
*
* Usage:
* <pre>
* Factorial f = new Factorial();
* long ten_over_two = f.over(10,2); //takes a while, because no cache yet
* System.out.println( f.factorial( 8)); //instantaneous
* </pre>
*
* Notes:
* <li>instantiate and use, it will reuse its own cache, so after a while it is supposed to run in O(1).
* Code is super fast (0.0ms)
* <li>For very big numbers, use Stirling's approximation, @see http://en.wikipedia.org/wiki/Stirling%27s_approximation
* <li>Overflow is a big problem.
* the max value is 9223372036854775807 < 21!,
* 1.7976931348623157E308 <171!
* Either it is super fast or it breaks, the combinations grow too fast, and for 2k~n it is biggest;
* it cannot even handle (29 14) but for k/n small it can function barely.
*
*/
public class Factorial {
static private HashMap<Integer, Double> cache = new HashMap<Integer, Double> ();
public Double factorial(int n)
{
Integer N = Integer.valueOf( n );
if(cache.containsKey( N) ) return cache.get( N );
if(n<1) {
cache.put( N, 0D );
return 0D;
}
if(n==1) {
cache.put( N, 1D );
return 1D;
}
Double temp = 1D;
for (int i = 1; i < n; i++) {
temp*=i;
cache.put( i, temp );
}
Double prev = factorial(n-1);
if(prev>Double.MAX_VALUE/n) throw new IllegalArgumentException("factorial("+n+") overflow");
return n*prev;
}
// 52!/47! = 52*51*50*49*48 so extra(52,47) = 52*(51,47) ... and (47,47)=1
private Double extra(int n, int minus)
{
if(n<=minus) return 1D;
Double result = extra(n-1,minus);
if(result>Double.MAX_VALUE/n) throw new IllegalArgumentException( "overflow");
return n* result;
}
/**
*
* This method computes (n k) = n!/(k! * (n-k)!)
* <p/>
* Example usage:
* <li>over(52,2) = 2598960
* <li>over(5,2) = 10
* <li>over(2,2) = 1 //ignoring swaps
*
* @param n the big set
* @param k how much to draw out of n
* @return all combinations of drawing k out of n, ignoring swaps. Simply double the amount for regarding swaps.
*
*/
public Double over(int n, int k){
if(k==n) return 1D;
if(k>n) return 0D;
if(k<=0) return 0D;
//k>=1 and n>k
if(k<n-k) return over(n,n-k);
Double over = extra(n,k);
Double divisor = factorial(k);
return over/divisor;
}
/**
* basic test and usage
*/
public static void main(String[] args)
{
Double fuc= new Factorial().factorial( 170); //takes a while, because no cache yet
System.out.println(fuc);
{
int n = 100;
int k = 50;
double start = System.currentTimeMillis();
//Factorial f = new Factorial();
Double ten_over_two = new Factorial().over(n,k); //takes a while, because no cache yet
double delta = System.currentTimeMillis()-start;
System.out.println( "ten_over_two=" +ten_over_two + ": in " + delta + "ms");
}
{
int n = 8;
double start = System.currentTimeMillis();
Factorial f = new Factorial();
System.out.println( f.factorial( n)); //instantaneous
double delta = System.currentTimeMillis()-start;
System.out.println( "factorial("+n+") in " + delta + "ms");
}
}
}//~class}
Not sure what the question is but here goes. Let's assume you are trying to have a computer perform the calculation. 1M! would have an exponent of < 10^6M. Not sure about software out there for that large of an exponent. If you take the log to calculate then the number is just < 6M. Logs also have the advantage of being addition when multiplying, so doing a M of addition is not that big of a deal. You could also set the function to be between i and k which would speed things up, assuming a fairly large size of i. Additionally you could look up the value for 100K! (2.824229408 × 10 456573 the log of which is 456573.450899970914458) and then just iterate on values larger than that.