# Numerical problems floating-point arithmetic

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.

• One possible solution may be using the module mpmath in python.See mpmath.org or code.activestate.com/recipes/… – Dieter Kadelka Feb 12 at 12:27
• First observation is that you can remove the denominator outside the brace by decreasing expoents in the brace by $2$. Second observstion is that is better to take the difference of teo sums than summing over differences; summing what is to the right of the minus sign amounts to multiplying with the number of summands – Manfred Weis Feb 12 at 12:43
• If you further do the exponentiation separately for numerators an denominator, then you can reduce the powers of the deminator by $1$ if you also remove the $z$ factor of the numerator outside the summation. As a general advice: simplify as much as possible before looking for more elaborate techniques. – Manfred Weis Feb 12 at 12:51
• en.wikipedia.org/wiki/Interval_arithmetic – Steve Huntsman Feb 12 at 13:45
• @ManfredWeis: Thanks for your suggestions! I tried implementing see suggestions for reducing the term outside the bracks. Unfortunately, it does not seem to help much. I am not sure I understand your second observation but just to clarify, I only substract the term with the log once at the end. – Julian Karch Feb 12 at 19:05

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.