suma szeregu e^-x

Cześć, muszę napisać prosty programik do sumowania szeregu potęgowego e^-x z zadaną dokładnością. Jest to szereg naprzemienny, dla x>0. Chodzi o to, że rożnica od pewnego momentu miedzy funkcja biblioteczna a mają jest dość znacząca i nie wiem gdzie mam błąd.

long double minusexpo( double x,double epsilon) {

     int k=0;
     long  double S[20000] = {1}; 
     long double y =1;

     while(abs(S[k] *(-x/(k+1))) > epsilon) {

	  S[k+1]= S[k]*(-x/(k+1));
	  y = y+S[k+1];
	  k++;
    };
  return y;
};

output (cmath i moja funkcja,epsilon 1e-15):

 x: 3
cmath: 0.0497870683678639
moja: 0.0497870683678647
x: 6
cmath: 0.00247875217666636
moja: 0.00247875217666945
x: 9
cmath: 0.00012340980408668
moja: 0.000123409803898621
x: 12
cmath: 6.14421235332821e-006
moja: 6.14421051805556e-006
x: 15
cmath: 3.05902320501826e-007
moja: 3.05910001953705e-007
x: 18
cmath: 1.52299797447126e-008
moja: 1.59837207894053e-008
x: 21
cmath: 7.58256042791191e-010
moja: 1.50843165493256e-008
x: 24
cmath: 3.7751345442791e-011
moja: 3.44305353678885e-007
x: 27
cmath: 1.87952881653908e-012
moja: -5.864302722932e-006
x: 30
cmath: 9.35762296884017e-014
moja: 6.10304247944361e-006
x: 33
cmath: 4.6588861451034e-015
moja: 0.00185404651646217
x: 36
cmath: 2.31952283024357e-016
moja: -0.0331838429246613
x: 39
cmath: 1.15482241730158e-017
moja: -0.072867218513009

1. Nie trzeba ci w tym kodzie z tablicy, skoro nie używasz starych elementów ponownie.
2. Powinno być fabs a nie abs.
3. Algorytm jest jak widać numerycznie niestabilny. Poszukaj czegoś bardziej stabilnego numerycznie.