Strangely low numerical precision

[jaskij]

The program is supposed to compute the values of Euler's Gamma Function using an infinite product (formula #15 in the link), and it does so decently for low values, but the error gets too big for bigger values.
Using both Mathematica and Maxima for reference this is what I get for Gamma(19.9):
Gamma(19.9)~9.0406e16
Mathematica and Maxima difference: 9.6e2
Difference between my functions and Mathematica: 1.79e11
Difference between both of my functions: 1.21e5
What surprises me as well is similarity of results for both of my functions, because the second one gets to multiply values like 1e66 and 1e-61, while the first doesn't.
Compiled under Win7/Cygwin using g++ version 4.5.3 with -O2 flag.
Just now run it without -O2, will see if it makes a difference.

My question is: is this error due to my mistake, or is it due to inner workings of C++?
Code below.

long int const MAX_ITERS=1e8

long double eulerGamma(long double xVal){
	long double base=exp(-GAMMA*xVal)/xVal;
	long double product=1;
	for(int i=1;i<MAX_ITERS;i++){
		product*=exp(xVal/i)/(1+xVal/i);
	}
	return base*product;
}

long double eulerGamma2(long double xVal){
	long double base=exp(-GAMMA*xVal)/xVal;
	long double product=1;
	long double expon=0;
	for(int i=1;i<MAX_ITERS;i++){
		product*=(long double)i/((long double)i+xVal);
		expon+=(long double)1/(long double)i;
	}
	return base*exp(xVal*expon)*product;
}

Thanks in advance.

[nezachem]

The problem is that a Weierstrass product converges very poorly (my gut feeling is that it is worse than linear; correct me if I am wrong). That is, even 1e8 iterations may not be enough.

[jaskij]

For 19.9, the current product element was 1+1.8e-10, don't really know what to think about this.

Seems it really was due to the Weierstrass product, and since the formula is the part of the assignment, I'm gonna leave it as is.

Thanks nezachem!