0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ One possible solution may be using the module mpmath in python.See mpmath.org or code.activestate.com/recipes/… $\endgroup$ – Dieter Kadelka Feb 12 at 12:27
  • $\begingroup$ 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 $\endgroup$ – Manfred Weis Feb 12 at 12:43
  • $\begingroup$ 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. $\endgroup$ – Manfred Weis Feb 12 at 12:51
  • $\begingroup$ en.wikipedia.org/wiki/Interval_arithmetic $\endgroup$ – Steve Huntsman Feb 12 at 13:45
  • $\begingroup$ @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. $\endgroup$ – Julian Karch Feb 12 at 19:05
5
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Great, this solves my problem thanks! I have no clue about the first equal sign, which translates the finite into an infinite sum. What is happening there? Also, I do not understand how you arrive at the errors. I guess this is related to my lack of knowledge about the Taylor series. I tried reading up on that but was unsuccessful. Could you point me to resources that will allow me to understand your answer? $\endgroup$ – Julian Karch Feb 13 at 18:27
  • $\begingroup$ Oh, and sorry for posting on MathOverflow. I mixed up the two. I actually wanted to post on MSE. I am aware that this is not research level. $\endgroup$ – Julian Karch Feb 13 at 18:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.