# Issues with real precision

83 views

### Fernando Gouvea

Feb 23, 2024, 5:00:50 PMFeb 23
to sage-support

In an introductory probability class, one computes the probability of getting all of n possible coupons in r individual purchases. The naive approach with inclusion-exclusion leads to the awful formula

f(n,r) = \sum_{k=0}^n \binom{n,k} \frac{(n-k)^r}{n^r}

Computing this in GP is straightforward:

(10:32) gp > g(n,k,r)=(-1)^(k)*binomial(n,k)*(n-k)^r/n^r
%3 = (n,k,r)->(-1)^(k)*binomial(n,k)*(n-k)^r/n^r
(16:43) gp > f(n,r)=sum(k=0,n,1.0*g(n,k,r))
%4 = (n,r)->sum(k=0,n,1.0*g(n,k,r))
(16:43) gp > f(6,12)
%5 = 0.43781568062117902081322291656082236787
(16:43) gp > f(50,100)
%6 = 0.00016616318861823418331310572392871711604
(16:45) gp > f(50,200)
%7 = 0.39828703196689437232172533821416865689
(16:47) gp > f(365,2236)

Doing it in SageMath seems to be much more delicate. This, for example, fails:

sage: g(n,k,r)=(-1)^(k)*binomial(n,k)*(n-k)^r/n^r
sage: f(n,r)=sum(1.0*g(n,k,r),k,0,n)
sage: f(365,2000).unhold()
0.0

It works better if I rewrite the function with ((n-k)/n)^r, of course. I presume the issue is that the default precision for real numbers is not enough to handle the computation. How do I tell SageMath to use a higher precision?

Fernando

--
=============================================================
Fernando Q. Gouvea         http://www.colby.edu/~fqgouvea
Carter Professor of Mathematics
Dept. of Mathematics
Colby College
5836 Mayflower Hill
Waterville, ME 04901

Reached by phone, a Westbury woman who identified herself as Onorato's
granddaughter said he was a retired businessman. He lived alone and
never married or had children, she said.
-- New York Newsday, December 27, 2003



### dim...@gmail.com

Feb 23, 2024, 5:16:37 PMFeb 23
here is how to do it exactly, and convert the result to a float.

sage: sage: g(n,k,r)=(-1)^(k)*binomial(n,k)*(n-k)^r/n^r
....: sage: f(n,r)=sum(g(n,k,r),k,0,n)
....: sage: f(50,200).unhold().n()
0.398287031966894
....: sage: f(365,2000).unhold().n()
0.216119451633218

There should be less brute-force ways too,
bu this is the 1st that comes to mind.

HTH
Dima

>
> Fernando
>
> --
> =============================================================
> Fernando Q. Gouveahttp://www.colby.edu/~fqgouvea
> Carter Professor of Mathematics
> Dept. of Mathematics
> Colby College
> 5836 Mayflower Hill
> Waterville, ME 04901
>
> Reached by phone, a Westbury woman who identified herself as Onorato's
> granddaughter said he was a retired businessman. He lived alone and
> never married or had children, she said.
> -- New York Newsday, December 27, 2003
>
> --
> You received this message because you are subscribed to the Google Groups "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/fd3261d8-3ee5-48a8-ba6e-40530790b92f%40colby.edu.
signature.asc

### Dima Pasechnik

Feb 23, 2024, 5:23:20 PMFeb 23
to sage-support
the normal Python way, without any symbolic sum, would be like this:

sage: sage: g(n,k,r)=(-1)^(k)*binomial(n,k)*(n-k)^r/n^r
....: sage: def f(n,r): return math.fsum([1.0*g(n,k,r) for k in range(n+1)])
....: sage: f(365,2000)
0.21611945163321847

HTH
Dima

### Emmanuel Charpentier

Feb 24, 2024, 8:54:51 AMFeb 24
to sage-support

Le vendredi 23 février 2024 à 23:23:20 UTC+1, Dima Pasechnik a écrit :

[ Snip…]

the normal Python way, without any symbolic sum, would be like this:

sage: sage: g(n,k,r)=(-1)^(k)*binomial(n,k)*(n-k)^r/n^r
....: sage: def f(n,r): return math.fsum([1.0*g(n,k,r) for k in range(n+1)])
....: sage: f(365,2000)
0.21611945163321847

This works in Sagemath. It wouldn’t work in Python : the range of the magnitudes of the terms of the alternating sum are way too large for the precision of Python’s floats, necessary if you want to use math.comb. Programming this in Python would need some serious analytical work, or using a multiple-precision integer library, which Sage does for you…

[ Re-snip... ]

HTH,

### Gareth Ma

Feb 24, 2024, 12:11:25 PMFeb 24
to sage-support
Note that you can wrap it in Decimal or Fraction, which are both builtin Python libraries.

### Emmanuel Charpentier

Feb 24, 2024, 5:36:34 PMFeb 24
to sage-support

Le samedi 24 février 2024 à 18:11:25 UTC+1, Gareth Ma a écrit :

Note that you can wrap it in Decimal or Fraction, which are both builtin Python libraries.

I stand corrected. Python ints are bignums, not 32-bits integers. and fractions.Fractions are rationals. Ordinary integer arithmetic wroks with these rationals. We just need a binomial coefficient function returning an int (math.comb returns a float) :

>>> from math import factorial >>> from fractions import Fraction >>> def ibin(n, k): return Fraction(factorial(n), factorial(k)*factorial(n-k)) if n>=k else 0 ... >>> def g(n, k, r): return (-1)**k*ibin(n,k)*Fraction((n-k)**r, n**r) ... >>> def f(n, r): return float(sum([g(n,k,r) for k in range(n+1)])) ... >>> f(365, 2000) 0.21611945163321844

A tad more muscular than using Sage…

Thanks a lot ! I learned something useful today…