I am trying to calculate the following function in floating-point arithmetic.

$$f(c,z)=\frac{(c-1)z}{(z-1)^2}\left( \sum_{k=2}^{c-1}\frac{1}{c-k}\left(\frac{z-1}{z}\right)^k-\left(\frac{z-1}{z}\right)^c\log(1-z)\right)$$

where $z\in(0,1)$ and $c \in \mathbb{N}$, and $c>1$.

The following implementation, which I exemplary display as Matlab code, works for some inputs.

```
function res = hypergeo(c,z)
theSum = 0;
shared = (z-1)/z;
for k=2:(c-1)
theSum = theSum+shared^k/(c-k);
end
prefactor = (c-1)*z*(z-1)^(-2);
res = prefactor*(theSum-shared^c*log(1-z));
end
```

However, for example for $c=100$, and $z=0.1$, it returns -4.5288e+79, which is clearly wrong. I know that the correct answer for this case is $1.001002$.

The problems seem to occur if $z$ is small or $c$ is large. This leads to the terms $\left(\frac{z-1}{z}\right)^k$ and $\left(\frac{z-1}{z}\right)^c$ becoming quite large. For the example, $\left(\frac{0.1-1}{0.1}\right)^{100}=-99^{100}$. This leads me to believe that the reason for the function to return the wrong result is some kind of error accumulation due to the finite precision of floating-point arithmetic. Since I have many subtractions in the formula it might be a loss of significance (https://en.wikipedia.org/wiki/Loss_of_significance).

Does anybody see a way to transform the expression such that those numeric problems do not occur anymore? I tried the brute-force solution of increasing the number of bits used (in the R implementation) but this did not resolve the problem. My intuition is that I should somehow avoid those exponential terms but I do not know how.

**UPDATE:**
I updated the code according to the suggestions of @ManfredWeis. It now reads

```
function res = hypergeo(c,z)
theSum = 0;
for k=2:(c-1)
theSum = theSum + (z-1)^(k-2)/(z^(k-1)*(c-k));
end
res = (c-1)*(theSum-(z-1)^(c-2)/z^(c-1)*log(1-z));
end
```

Unfortunately, this did not help much. For $c=100,z=0.1$, I get -4.58e+79.