On Mon, Sep 13, 2010 at 9:36 AM, Dan Drake <dr...@kaist.edu> wrote:
> I've found something strange with mpmath's hypergeometric functions:
>
> sage: f = lambda n: sum(binomial(n,k)^2 * 2^k for k in range(n+1))
> sage: g = lambda n: (-f(n) + 8*f(n-1) - 3*f(n-2))/8
>
> Now define:
>
> sage: h = lambda n: catalan_number(n)*mpmath.hyp3f2(-n, -n+1, 2, -2*n, 1, -1)
>
> Now, I have a proof that g(n+1) equals h(n) when n >= 1. But I find this:
>
> sage: g(25)
> 60529607948876437
> sage: h(24)
> mpf('60529607948876432.0')
> sage: g(25) - h(24)
> mpf('8.0')
>
> When I print the two numbers, it looks like they differ by 5 -- but when
> I subtract them, they differ by 8. Why is that?
This happens because g(25) is a Rational, so mpmath rounds it to a
floating-point value (with 53-bit precision) upon coercion. Note:
sage: type(g(25))
<type 'sage.rings.rational.Rational'>
sage: Integer(Rational(60529607948876437).n(53))
60529607948876440
sage: Integer(Rational(60529607948876437).n(53)) - 60529607948876432
8
sage: Integer(g(25)) - h(24)
mpf('5.0')
I could (indeed, should) fix mpmath to special-case Sage Rationals
with power-of-two denominators to convert losslessly, like Integers
already do.
Fredrik
--
To post to this group, send an email to sage-devel@googlegroups.com
To unsubscribe from this group, send an email to
sage-devel+unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URL: http://www.sagemath.org
