Trying to calculate with the Bernoulli formula

[yabi2943993]

I am willing to calculate with the Bernoulli formula.
For that I have made the following code:

#include <iostream>
#include <math.h>

using namespace std;

int berechneWK();

int main()
{
    cout<<berechneWK()<<endl;
}
   
int berechneWK()
{
    int n,k,p,n2,k2,n3,k3,gemischt;
    
    cin>>n;
    cin>>k;
    cin>>p;
    
    long nfakultaet=n,kfakultaet=k,gemischtefakultaet=n-k;
    
    n2=n;
    k2=k;
    k3=k;
    n3=n;
    gemischt=n2-k2;
    while(n>0)
    {
        nfakultaet*=n-1;
        n--;
    }
    while(k>0)
    {
        kfakultaet*=k-1;
        k--;
    }
    while(gemischt>0)
    {
        gemischtefakultaet*=gemischt-1;
        gemischt--;
    }
    
    return ((nfakultaet)/(kfakultaet*gemischtefakultaet)) * pow(p,k3) * pow(1-p,n3-k3);
}

As I run this, I receive the error floating point exception (core dumped). I think the problem are the loops but I am not sure. Any idea how to solve that or tipp on how to calculate the probability more easy?

[jlb]

What exactly are you inputting into this program?

[helios]

First, you can't use int to calculate factorials. the largest factorial you can compute is 12!. You will probably need to find a way to simplify away factors from the multiplication.

Second, the way you're computing factorials is incorrect anyway. Look at what you're doing. The last value of n on lines 28-32 is n == 1, and you're multiplying nfakultaet by n - 1. The result is just zero. This is the immediate cause of the exception, since the divisor in (nfakultaet)/(kfakultaet*gemischtefakultaet) is zero.

[lastchance]

Work in doubles and take logs. log(n!) is perfectly tractable as log(n) + log(n-1) + ... + log(1).

At the end of the final calculation you can reverse with exp().

[lastchance]

Something like

#include <iostream>
#include <cmath>
using namespace std;

int main()
{
   int n, k;
   double p;
   cout << "Input the number of trials (n): ";   cin >> n;
   cout << "Input the probability of success on each trial (p): ";   cin >> p;
   cout << "Input the number of successes (k): ";   cin >> k;
   
   double log_answer = k * log( p ) + ( n - k ) * log( 1.0 - p );
   for ( int nn = n, kk = k; kk; nn--, kk-- ) log_answer += log( (double)nn / kk );
   cout << "Probability is " << exp( log_answer );
}

Input the number of trials (n): 100
Input the probability of success on each trial (p): 0.3
Input the number of successes (k): 35
Probability is 0.0467797

[dhayden]

You can compute large, exact values of the combinations function nCk by alternating the multiplication and division. The beauty is that each division results in an integer result. Here's an example, with debug code to show what it does:

#include <iostream>
unsigned long long comb(unsigned long long N, unsigned long long K)
{
    if (K > N || N==0 || K==0) return -1; // invalid inputs
    if (K > N-K) {
	K = N-K;
    }
    unsigned long long result = N;
    std::cout << N << 'C' << K << " = " << result;
    for (unsigned long long i=1; i<K;) {
	std::cout << '*' << N-i << '/' << i+1;
	result *= N-i;
	result /= ++i;
    }
    std::cout << '\n';
    return result;
}

int main()
{   
    unsigned long long n,k;
    while (std::cin >> n >> k) {
	std::cout << comb(n,k) << '\n';
    }
}

50000 40
50000C40 = 50000*49999/2*49998/3*49997/4*49996/5*49995/6*49994/7*49993/8*49992/9*49991/10*49990/11*49989/12*49988/13*49987/14*49986/15*49985/16*49984/17*49983/18*49982/19*49981/20*49980/21*49979/22*49978/23*49977/24*49976/25*49975/26*49974/27*49973/28*49972/29*49971/30*49970/31*49969/32*49968/33*49967/34*49966/35*49965/36*49964/37*49963/38*49962/39*49961/40
38691986793473693

[helios]

Did I miss something? How is 49999/2 an integer?

[dhayden]

The order of operatations is left to right, so it's (50000*49999)/2, which is an integer.

[helios]

Ah, duh.