Numerical problems floating-point arithmetic

Question

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.

Answer 1

It is actually pretty simple if you are comfortable with Taylor series (definitely not MO level, so ask on MSE next time). Let $w=\frac{z-1}{z}$. If $|w|<1$, you are in no trouble computing the expression as it is. So let's consider the case $|w|>1$. Then $\frac1{1-z}=1-\frac 1w$, so your expression in parentheses (the one you really have trouble with) becomes $$ w^c(\frac 1w+\frac 1{2w^2}+\dots+\frac 1{(c-2)w^{c-2}}+\log(1-\frac 1w)) \\\ =-w^c\sum_{k=c-1}^\infty \frac 1{kw^k} =-w\sum_{m=0}^\infty \frac 1{(m+c-1)w^m} $$ If $|w|>2$, say, the series converges pretty fast, so your real trouble is not $z=0.1$ but $z\approx \frac 12$, where the series converges not so fast. However, let's say that you have $15$ decimal digit float point precision. Then your error with direct computation will be, roughly speaking $10^{-15}c|w|^c$ and the number of terms in the series that you should take to make the error coming from the truncation about $10^{-15}$ is going to be $8c$ if $|w|^c>100$, say. So I suggest as a rule of thumb comparing $|w|^c$ to $100$ and if it is less than that, then do the direct computation but if it is above that to take $8c$ terms in the infinite series. The guaranteed (relative) error is then about $10^{-13}c$ which is $<10^{-8}$ (the handheld calculator precision) for all $c<10^5$. If that is not enough, implement higher precision arithmetic yourself or use some ready package and adjust the splitting into cases accordingly.