Coding complex equations into C++
I have here equations from a paper by <PERSON>. They're for counting the events in a <PERSON>-distributed data set.
$$\begin{align} \Pr(N(t)=n) &= \sum_{j=n}^\infty{\frac{(-1)^{j+n}(\lambda t^c)^j \alpha_j^n}{\Gamma(cj+1)}} \ E(N(t)) &= \sum_{n=1}^\infty{\sum_{j=n}^\infty{\frac{n(-1)^{j+n}(\lambda t^c)^j \alpha_j^n}{\Gamma(cj+1)}}} \end{align}$$
Where $Pr(N(t) = n)$ is the probability that the count will be $n$, and $E(N(t))$ is the expected value. 346
Below are the other vectors.
$$\begin{align} \alpha_m^l &= \sum_{m=0}^{l-1}{\frac{\Gamma(cm+1)\Gamma(cl-cm+1)}{\Gamma(m+1)\Gamma(l-m+1)}} \ \alpha_j^0 &= \frac{\Gamma(cj+1)}{\Gamma(j+1)}, \quad j=0,1,2,\ldots \ \alpha_j^{n+1} &= \sum_{m=n}^{j-1}{\alpha_m^n \frac{\Gamma(cj-cm+1)}{\Gamma(j-m+1)}} \end{align}$$
I've programmed these into C++, but I'm not getting the right results. I know because when I get the expected value, I get an extremely negative number (too many zeroes to count). I'd like some input on where I got my code wrong?
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
using namespace std;
double alpha_m(int n, double shape){
double sum = 0;
for(int m = 0; m < n-1; m += 1){
//printf("%i < %i\n", m, n-1);
double num = tgamma((shape*m) + 1) * tgamma((shape*n) - (shape*m) + 1);
// printf("num %f , shape %f, m %i\n", num, shape, m);
// system("pause");
double denom = tgamma(m + 1) * tgamma(n - m + 1);
//printf("denum %f \n", denom);
sum += num/denom;
//printf("sum%f\n", sum);
}
// printf("sum%f\n", sum);
// printf("done");
return sum;
}
double alpha_j(int n, int j, double shape){
//printf("n: %i", n);
//system("pause");
double sum = 0;
if(n == 0){
double result = tgamma((shape*j) + 1)/tgamma(j + 1);
return result;
} else if(n >= 1) {
for(int m = n-1; m < j-1; m += 1){
//printf("%i < %i\n", m, j-1);
double result = tgamma((shape*j) - (shape*m) + 1)/tgamma(j - m + 1);
// printf("result %f\n", result);
sum += result;
// printf("sum %f\n", sum);
}
double alpha = alpha_m(n,shape);
// printf("sum %f, alpha %f\n", sum, alpha);
sum = sum*alpha;
// printf("sum %f\n", sum);
return sum;
}
}
int main()
{
/* double x;
printf("libm");
for(x=1.0; x <= 10.0; x+=1.0) {
printf("%15.8lf\n", tgamma(x/3. 65