Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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?

share|improve this question

closed as off topic by Did, AD., J. M., Grigory M, Asaf Karagila Dec 8 '11 at 13:56

Questions on Mathematics Stack Exchange are expected to relate to math within the scope defined by the community. Consider editing the question or leaving comments for improvement if you believe the question can be reworded to fit within the scope. Read more about reopening questions here.If this question can be reworded to fit the rules in the help center, please edit the question.

    
Since the OP specifically says "using GMP-GNU", I'm not sure whether it's appropriate to remove that from the title. I understand the broader appeal, though, and that it may otherwise not be on-topic for math.SE (if it is, even now?). –  r.e.s. Dec 8 '11 at 7:44
    
What have you tried? –  lhf Dec 8 '11 at 11:08
    
mpz_bin_uiui can be used for the binomial calculation. But how can it be possible to calculate mpf_t number in GMP? –  user12290 Dec 8 '11 at 12:51
    
This might fit on Computational Science. –  David Z Dec 9 '11 at 0:42
add comment

3 Answers 3

up vote 2 down vote accepted

What you have is called a Bernstein polynomial. You can use the recursion formula for binomial coefficients,

$$\binom{n}{k+1}=\frac{n-k}{k+1}\binom{n}{k}$$

along with repeated squaring to evaluate your polynomial.

share|improve this answer
    
Mathematica has the function built-in: BernsteinBasis[100, i, a] –  J. M. Dec 8 '11 at 6:46
add comment

The OP specifically wants to use GMP-GNU. For example, with Sage/gmpy, the following prints all the terms and verifies that they sum to 1:

import gmpy
terms = [long(gmpy.bincoef(100L, i))*0.001^i*0.999^(100-i) for i in [0..100]]
print terms, sum(terms)

(The gmpy library may need to be installed first, using the command ./sage -i gmpy-1.0.1, the current folder being the one in which Sage is installed.)

share|improve this answer
add comment

So what's the problem?

In Maple:

seq(binomial(100,i)0.001^i(1-0.001)^(100-i),i=0..100);

.9047921471, .9056978450e-1, .4487692025e-2, .1467446842e-3, .3562120711e-5, .6846117884e-7, .1085053719e-8, .1458530667e-10, .1697239140e-12, .1736692257e-14, .1581971926e-16, .1295636303e-18, .9618921502e-21, .6517787728e-23, .4054393911e-25, .2326846021e-27, .1237374323e-29, .6120205097e-32, .2824919491e-34, .1220396176e-36, .4947552068e-39, .1886668408e-41, .6781636376e-44, .2302161455e-46, .7393494828e-49, .2249872300e-51, .6496512762e-54, .1782308028e-56, .4651383028e-59, .1155982113e-61, .2738562897e-64, .6190041746e-67, .1336063815e-69, .2755857052e-72, .5436095580e-75, .1026118428e-77, .1854568396e-80, .3211113204e-83, .5329016698e-86, .8480250385e-89, .1294532717e-91, .1896334456e-94, .2666564967e-97, .3600362412e-100, .4668774626e-103, .5815846492e-106, .6960690192e-109, .8005394125e-112, .8848137484e-115, .9399259447e-118, .9596841473e-121, .9418086201e-124, .8883618693e-127, .8053595052e-130, .7016627135e-133, .5874326112e-136, .4725165792e-139, .3651147549e-142, .2709594846e-145, .1930794922e-148, .1320697228e-151, .8668978669e-155, .5458525753e-158, .3295739844e-161, .1907256854e-164, .1057384257e-167, .5612953713e-171, .2851215787e-174, .1385063312e-177, .6429911943e-181, .2850382814e-184, .1205592697e-187, .4860720196e-191, .1866252081e-194, .6816114246e-198, .2365284890e-201, .7788330733e-205, .2429961619e-208, .7172443883e-212, .1999388825e-215, .5253649314e-219, .1298495857e-222, .3011721635e-226, .6537982492e-230, .1324487611e-233, .2495648799e-237, .4357233045e-241, .7018658041e-245, .1037885096e-248, .1400796431e-252, .1713798326e-256, .1885180042e-260, .1846043912e-264, .1589584347e-268, .1184917943e-272, .7491183452e-277, .3905563612e-281, .1612153849e-285, .4940104950e-290, .99900e-295, .10e-299

share|improve this answer
add comment

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