I would like to calculate the following expression with large $m$:

$$\sum^{m}_{q=1} \frac{(-1)^{q+1}}{q+1} {{m}\choose{q}} e^{-\frac{q}{q+1}\Gamma}.$$

But, due to the binomial, the computer gets stuck when $m$ grows large. In this problem we have that $\Gamma > 0$. I am trying to find a simplification or a way around, but I didn't find anything that could help me. Can anyone give me some hints?

  • $\begingroup$ For $m = 100$ and $\Gamma = 1/2$, Mathematica gives -0.963334247861897822795. I think almost all of those digits are correct. How large is $m$? $\endgroup$ – irchans Feb 15 at 19:56
  • $\begingroup$ For $m=4096$ and $\Gamma=1/2$, Mathematica gives the result -0.99792 after about 5 seconds of computation. $\endgroup$ – irchans Feb 15 at 20:00
  • $\begingroup$ In comments below, OP says they need $m=2^{13}$, so floating point computation is out of the picture. Maybe the only way forward is an efficient-as-possible arbitrary precision implementation but even then those combinations get very large. $\endgroup$ – parsiad Feb 15 at 20:05
  • $\begingroup$ Hello friends! My Mathematica gives me impossible values from $m=2^7$. This results comes from a probability which has to be between 0 and 1. $Gamma>0$, thus, it can assume values like 1, 10, 100... $\endgroup$ – Cláudio Ferreira Dias Feb 18 at 10:50

Computing the binomial coefficients efficiently

If by "gets stuck" you mean that the computation is slow, I would guess that you are computing the binomial term inefficiently. Indeed, you shouldn't recompute the binomial term for every summand, but instead use the fact that $$ \binom{m}{q}=\frac{m!}{q!\left(m-q\right)!}=\frac{m-\left(q-1\right)}{q}\frac{m!}{\left(q-1\right)!\left(m-\left(q-1\right)\right)!}=\frac{m-q+1}{q}\binom{m}{q-1}. $$ Defining $$ C_{q}=\frac{m-q+1}{q}C_{q-1}\text{ if }q\geq1\qquad\text{and}\qquad C_{0}=1, $$ it follows from the previous claim that $C_{q}=\binom{m}{q}$. Therefore, you can rewrite the sum you are interested as $$ S\equiv \sum_{q=1}^{m}\frac{\left(-1\right)^{q+1}}{q+1}C_{q}\exp\left(-\frac{q}{q+1}\Gamma\right). $$

Removing some terms by symmetry

We can use the fact that $C_{q}=C_{m-q}$ to reduce the number of terms. Note that $$ S-1=\sum_{q=0}^{m}\frac{\left(-1\right)^{q+1}}{q+1}\exp\left(-\frac{q}{q+1}\Gamma\right)C_{q}. $$ Assuming $m=2j+1$ is odd, we get $$ S-1=\sum_{q=0}^{j}\left(-1\right)^{q+1}\left(\frac{1}{q+1}\exp\left(-\frac{q}{q+1}\Gamma\right)-\frac{1}{m-q+1}\exp\left(-\frac{m-q}{m-q+1}\Gamma\right)\right)C_{q}. $$ Assuming $m=2j$ is even, we get \begin{multline*} S-1=\frac{\left(-1\right)^{j+1}}{j+1}\exp\left(-\frac{j}{j+1}\Gamma\right)C_{j}\\ +\sum_{q=0}^{j}\left(-1\right)^{q+1}\left(\frac{1}{q+1}\exp\left(-\frac{q}{q+1}\Gamma\right)+\frac{1}{m-q+1}\exp\left(-\frac{m-q}{m-q+1}\Gamma\right)\right)C_{q}. \end{multline*}

  • $\begingroup$ Hello Parsiad, Thanks for your suggestion! Your instructions were very clear and helpful! $\endgroup$ – Cláudio Ferreira Dias Feb 15 at 17:39
  • $\begingroup$ While the bit about computing the binomial term efficiently is fine, I'm going to remove what I wrote about approximating the sum and think about this problem a little more. The approach I used there is not practical. $\endgroup$ – parsiad Feb 15 at 17:49
  • $\begingroup$ Ok. I am looking forward to it. The second part seemed very interesting. I think I could choose an arbitrary k and verify the error. With luck, I could get small enough error on a few amounts of trials. $\endgroup$ – Cláudio Ferreira Dias Feb 15 at 18:10
  • $\begingroup$ @CláudioFerreiraDias: Just edited. I think this should work better. I will code it up to check. $\endgroup$ – parsiad Feb 15 at 18:23
  • $\begingroup$ @CláudioFerreiraDias: Added code. Even at $m = 50$ I'm only seeing a few digits of accuracy. At $m = 100$, the approximation blows up. What values of $m$ do you need this to work for? $\endgroup$ – parsiad Feb 15 at 19:45

This computation can be done numerically without a problem up to $m=2^{13}$. It took about 15 seconds of cpu time using Mathematica on my 2 year old Mac laptop for $m=2^{13}$.

Here is an example of the Mathematica code:

s[m_, g_] := Sum[ (-1)^(q+1)/(q + 1) Binomial[ m , q] Exp[ - q g/(q + 1)], {q, 1, m}];
Print[ Timing[ s[10000, N[1/4, 10000]]//N]; 

The output for the program above is {27.7445,0.999574} indicating that it took 27 seconds to compute the answer. Note that ${1000\choose 500}$ has about 3000 digits, so the program used 10000 digits of precision. The running time is order $m^3$.

The answer is usually close to 1 when $0<q<1$ and $m> 2^{10}$.

I wrote the code in Python and got the same result for $m=2^{13}$ and $q=1/4$.

from mpmath import mp

mp.dps =5000;

m = 2**13;

mp.pretty = True

rS = mp.mpf(0);

g = mp.mpf(1)/mp.mpf(4);

for q in range(1, m+1):
    rS = rS + mp.mpf((-1)**(q+1))* mp.binomial(m, q)/(q+1)*mp.exp(-q*g/(q+1));

mp.nprint(rS, 10)
  • $\begingroup$ +1. Nice! I don't have Mathematica but it would be interesting to see if it gets any faster if you compute the binomial using the recurrence and use the symmetry of the binomial to reduce the number of terms in the summation (see my recently edited answer). $\endgroup$ – parsiad Feb 15 at 21:51
  • $\begingroup$ @parsiad I have got to go somewhere, but when I get back, I can try your tricks to improve the speed. $\endgroup$ – irchans Feb 15 at 22:12
  • $\begingroup$ For most of the values of $m$ and $q$ that I tried, the sum was between 0.9 and 1. I don't know why. Any thoughts? $\endgroup$ – irchans Feb 15 at 22:18
  • 1
    $\begingroup$ If you use default Matlab without any special toolboxes, all calculations will be limited to 16 digits of precision which is inadequate. With Mathematica, you can specify the number of digits needed by not using any floating point values (decimals), and writing N[ <integer value>, <precisions>] for each integer or fraction. That is why I wrote N[ 1/4, 10000] in my code above. You should be able to just cut and paste the code above into Mathematica and it will work. For this application, Mathematica was about twice the speed of Python when I ran both on the same problem. $\endgroup$ – irchans Feb 18 at 18:32
  • 1
    $\begingroup$ If you want to run higher values of $m$ in Python, you can write mp.dps = m/2 $\endgroup$ – irchans Feb 18 at 18:35

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.