How to calculate large numbers in c (large factorials and powers)-CodePudding

I want to calculate the probability of getting 40 or less 3's when throwing 234 dice. So I use binomial distribution, which uses large powers (up to 40 in this case) and also large factorials (up to 234! in this case). While calculating these large numbers I get this error:

"Floating point exception (core dumped)"

Is there any way to fix this error and calculate properly, since double and even long double is not enough to represent (1/6)^40 or (234!)?

My code:

int factorial(int n) 
{
    if(n == 0 or n == 1)
    {
        return 1;
    }
    else
    {
        return n * factorial(n - 1);
    }
}

double binomialDensity(int n, int k, double p)
{
    return (factorial(n) / (factorial(k) * factorial(n-k)) * pow(p, k) * pow((1 - p), (n-k)));
}

int main() 
{
    long double probability;
    for(int i = 0; i <= 40; i  )
    {
        probability  = binomialDensity(234, i, 1/6);
    }
    cout << probability << endl;

    return 0;
}

CodePudding user response：

You do not need intermediate results of such large magnitude. Consider that for example

factorial(n) / factorial(k) == (k 1) * (k 2) * .... * n

The right hand side is of much smaller magnitude than the individual terms on the left side.

The formula you used is good for maths because it is readable and convenient for analysis. To compute the result you should use something else. The following is using loops to accumulate the final result. It is certainly not the most numerically stable but it should explain the idea:

#include <cmath>
#include <iostream>

int factorial(int n) 
{
    if(n == 0 or n == 1)
    {
        return 1;
    }
    else
    {
        return n * factorial(n - 1);
    }
}


double binomialDensity(int n, int k, double p)
{
    return (factorial(n) / (factorial(k) * factorial(n-k)) * pow(p, k) * pow((1 - p), (n-k)));
}

double binomialDensity2(int n,int k, double p) {
    double result = 1;
    for (int i = k 1; i <= n;    i) result *= i;    // factorial(n) / factorial(k)
    for (int i = 1;   i <= n-k;  i) result /= i;    // factorial(n-k)
    for (int i = 0;   i < k;     i) result *=p;     // pow(p,k)
    for (int i = 0;   i < n-k;   i) result *=(1-p); // pow(p-1,n-k)
    return result;
}


int main() {
  int n = 10;
  double p = 0.5;
  for (int k=1; k < n;   k){
      std::cout << binomialDensity(n,k,p) - binomialDensity2(n,k,p) << "\n";
  }
}

Live Demo

main compares your solution with the one that uses loops in a range where your solution still works. binomialDensity2 is still using relatively large intermediate results. This can be mitigated by combining the loops that I kept seperated here for the sake of clarity.

CodePudding user response：

Numerical programming is all about avoiding techniques that produce large numbers, especially if they cancel during the calculation.

Using Pascal's triangle is the better approach here.

If you are worried about runtime performance with that approach, then program it with metaprogramming, constexpr programming &c. &c.