Evaluating $\binom{100}{i}a^i(1-a)^{(100-i)}$ in GMP-GNU

Question

I want to calculate

$\binom{100}{i}a^i(1-a)^{(100-i)}$ for different $i$ with $a=0.001$ using GMP-GNU.

How can this be done?

This was asked (and closed) in math.SE... – J. M. Dec 10 '11 at 05:53 — J. M., Dec 10 '11 at 05:53
Are R solutions allowed? :) – Roman Luštrik Dec 10 '11 at 06:40 — Roman Luštrik, Dec 10 '11 at 06:40

score 3 · Accepted Answer · answered Dec 11 '11 at 07:10

I hope this is what you were looking for. Here is a program to print out '1':

/* gcc -lgmp foo.c -o foo */
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

int main(int argc, char **argv) {
  int t, i, n=100;
  mpz_t fac_num, fac_den;
  mpq_t a, q, tq, sum;

  mpz_init(fac_num); mpz_init(fac_den); 
  mpq_init(a); mpq_init(q); mpq_init(tq); 
  mpq_init(sum);

  mpq_set_ui(sum, 0, 1);

  for (i=0; i<=n; i++) {

    mpq_set_ui(a, 1, 1000);
    mpq_set_ui(q, 999, 1000);
    mpq_set_ui(tq, 1, 1);

    mpz_set_ui(fac_num, 1);
    mpz_set_ui(fac_den, 1);
    for (t=0; t<i; t++) {
      mpz_mul_ui(fac_num, fac_num, n-t);
      mpz_mul_ui(fac_den, fac_den, t+1);
      mpq_mul(tq, tq, a);
    }

    for (; t<n; t++) 
      mpq_mul(tq, tq, q);

    mpq_set_ui(q, 1, 1);
    mpq_set_num(q, fac_num);
    mpq_set_den(q, fac_den);
    mpq_mul(q, q, tq);

    mpq_add(sum, sum, q);

  }

  mpq_canonicalize(sum);

  gmp_printf("%Qd\n", sum);

  mpz_clear(fac_num); mpz_clear(fac_den); 
  mpq_clear(a); mpq_clear(q); mpq_clear(tq);
  mpq_clear(sum);

}

There might be more efficient ways to do the factorial and you can do exponentiation by squaring, but I leave those for the reader.

Thank you for your help. As factorial is already for mpz_t integers, is it possible to use that function? — user12290, Dec 11 '11 at 18:01
You can use mpz_fac_ui and mpz_pow_ui for the factorial and (efficient) squared exponentiation. My way is slightly more efficient than calling mpz_fac_ui three times to calculate the choose function (unless there is an optimization gmp uses that I don't know about) but using the built in functions makes it cleaner. mpz_pow_ui is much more efficient than the above. — user834, Dec 11 '11 at 20:30

score 2 · Answer 2 · answered Dec 22 '11 at 20:09

2

I would also go for the evaluation of the logarithm of your expression log(binom(100,i)) + i*log(a)+(100-i)*log(1-a) and then take the exponential of the result.

Could you tell us why you need multiple precision arithmetic? What is the accuracy that you need on the final result? I think you can get quite far with good series approximations, certainly when a is close to zero...

answered Dec 22 '11 at 20:09

GertVdE

6,179
1
21
36

I need to multiply this with an 1000 bit integer. Hence I need accuracy as much as possible. – user12290 Dec 27 '11 at 05:48

score 1 · Answer 3 · answered Dec 22 '11 at 23:36

Depending on the accuracy you need, think about just approximating the result.

At least for the part {100 over i} you should be able to get a good approximation using techniques like Sterling's formula (http://en.wikipedia.org/wiki/Sterling's_approximation). The remainder can be computed relatively easily.

Evaluating $\binom{100}{i}a^i(1-a)^{(100-i)}$ in GMP-GNU

3 Answers3