7
$\begingroup$

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?

5 Answers 5

9

You need Stirling's approximation. It is very accurate for large factorials.

  • 0
    Stirling's approximation would not help at all. The same restrictions apply for this formula too. What's good with having to calculate 100000 in the exponent?2015-10-17
  • 0
    @PanosK.: I don't see any $N!$ there. The next step would be to use the normal approximation for the $k \choose i$ term2015-10-17
  • 0
    Yes there are no factorials there indeed. I am talking about computers. It does look better on paper and it's better when i do the calculations. but the complexity, or even better the time and memory needed to compute this number (100000) had no practical implementation in my mobile phone calculator for example. I can calculate max 480! With Stirling's approximation it's exactly the same.2015-10-17
  • 0
    And don't forget that O(n!) < O(n^n)2015-10-17
5

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.

1

Here is a good writeup about implementing Stirling's approximation, along with a reference implementation:

http://threebrothers.org/brendan/blog/stirlings-approximation-formula-clojure/

1
 import java.util.HashMap;

/**
 * Convenience methods for deterministic combinations.
 * 
 * Usage:
 * 
 
 * 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
 * 
* * Notes: *
  • 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) *
  • For very big numbers, use Stirling's approximation, @see http://en.wikipedia.org/wiki/Stirling%27s_approximation *
  • 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 cache = new HashMap (); 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)!) *

    * Example usage: *

  • over(52,2) = 2598960 *
  • over(5,2) = 10 *
  • 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
    • 0
      Why would that be faster ?2014-10-09
    • 0
      because it reuses already computed factors.2014-10-09
    0

    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.