# Simple precision issue [closed]

Considering binomial expansion

$1 = (p+(1-p))^n = \sum_{k=0}^n p^n (1-p)^{n-k},$

which is true for any real value of $p$, I'm getting diverging results when I numerically evaluate the above even with very high precision. The code I use is a simple one, say for $n=50$ and $p=10$:

SetPrecision[Sum[Binomial[n, k] p^k (1. - p)^(n - k) /.
{n -> 50., p -> 10.}, {k,0., 50.}], 1000] // N


For the above computation, each $p^k$ requires at most 50 digits of precision, and same with $(1-p)^k$ and $\binom{n}{k}$ which is at most $2^n$, so in total a couple of hundred digits should be plenty but even with 1000 digits I'm getting $-1.48779 \times 10^{49}$ instead of 1. What am I doing wrong, and how can I change my code to get the correct result 1?

• The computation as written is done at machine precision, before SetPrecision is applied. Compare with bsum[n_, p_] := Sum[Binomial[n, k] p^k (1 - p)^(n - k), {k, 0, n}] and evaluating bsum[50, 10], bsum[50, 10.] and bsum[50, SetPrecision[10, 80]]. – ilian Jun 23 '17 at 20:04
• Thanks, but your example evaluates the sum symbolically to 1 and then applies the parameters (doing nothing). I gave this example to make a point for a much more complex one (that cannot be done symbolically but suffers from the same precision issue as the one I gave). Could you write an example that indeed forces Mathematica to use numerical computation but also gives the correct answer 1.0000? Thanks! – MCH Jun 23 '17 at 21:05
• I think the result of bsum[50, 10.] shows there is no such thing as first evaluating the sum symbolically to 1... – ilian Jun 23 '17 at 21:11
• Yes, seems so. I'll try this trick on the more complex calculation and see if it still gives nonsense. Thanks again. – MCH Jun 23 '17 at 21:18
• Maybe you want something like this: Sum[Binomial[n, k] p^k (1.100 - p)^(n - k) /. {n -> 50.100, p -> 10.100}, {k, 0, 50.100}] -- Or Hold[Sum[Binomial[n, k] p^k (1. - p)^(n - k) /. {n -> 50., p -> 10.}, {k, 0., 50.}]] /. {x_Real /; x != 0 :> SetPrecision[x, 100], z_ /; z == 0 -> 0} // ReleaseHold – Michael E2 Jun 24 '17 at 0:34

## 1 Answer

You can only get high procession results when every number in the computation is a high-precison inexact number or is exact. If even one machine precision number, such as 1., gets included, the whole computation gets reduced to machine precision. So try this

With[{pr = 50},
N[Sum[Binomial[n, k] p^k (1 - p)^(n - k) /. {n -> 50, p -> 10}, {k, 0, 50}], pr]]


or better yet

With[{pr = 50}, N[Sum[Binomial[50, k] 10^k (1 - 10)^(50 - k), {k, 0, 50}], pr]]


Both give

1.0000000000000000000000000000000000000000000000000

• Thanks. To make sure that Mathematica doesn't "cheat" by symbolically evaluating the sum to the exact value 1 and then numerically estimate just 1, I tried to change 1-p to something like SetPrecision[1, 500] - p. Seems to work. – MCH Jun 24 '17 at 18:44