Why does Sage offer unusable SciPy?

hypergeometric

asked 2014-09-17 13:32:32 +0200

Peter Luschny
548 ●13 ●35 ●56

In the docs of special functions [1] there is a warning:

"SciPy’s versions are poorly documented and seem less accurate than the Maxima and PARI versions."

Less accurate? Where does Sage draw the border to unusable code?

I certainly do not expect such incorrect results produced with Sage:

A000262 = lambda n: (-1)^(n-1)*hypergeometric_U(1-n, 2, -1, "scipy")
[round(A000262(n).n(100)) for n in (0..19)]

[1, 1, 3, 13, 73, 501, 4051, 37633, 394353, 4596553, 58941091,
824073141, 12470162233, 202976401213, 3535017524400, 65573803186900,
1290434218670000, 26846616451200000, 588633468315000064,
13564373693600000000]

Peter

[1] http://www.sagemath.org/doc/reference...

[2] https://oeis.org/A000262

edit retag flag offensive close merge delete

add a comment

answered 2014-09-18 10:54:56 +0200

Dima
716 ●8 ●22 https://github.com/dim...

updated 2014-09-20 19:21:20 +0200

AFAIK, scipy calls a Fortran numerical routine to compute this. Naturally, it doesn't do arbitrary precision. The relative precision is probably still quite OK, I didn't check. If you try adding optional parameter prec, with prec not equal to 53, you get an error.

Not sure whether this is an acceptable solution, but perhaps Sage should return the result as RDF, not as an arbitrary precision real. (EDIT: this is now http://trac.sagemath.org/ticket/17011)

One might also just call hypergeometric instead. This will use Maxima (Maxima doesn't know this particular function, but it can do qFp, and this is precisely what is needed here). E.g. round(hypergeometric((-18,-19),(),1).n(100)) gives the correct 19th term of the sequence in question.

edit flag offensive delete link

Comments

Dima, yes of course it is simple to avoid the errors provided you are aware of them. As you noted U(1-n, 2, -1) = (-1)^(n-1)*hypergeometric([-n,-n+1],[],1) does the job. But this was not my question. I asked why Sage *offers* such dangerous functions. "Sage should return the result as RDF, not as an arbitrary precision real." Yes, at least this.

Peter Luschny ( 2014-09-20 13:13:40 +0200 )edit

Sage is also a system for doing research, that includes comparing different implementations. A lot of work went into picking suitable defaults. Why shouldn't we expose the other implementations, too? Maybe SciPy improves in some future version...

Volker Braun ( 2014-09-21 19:13:01 +0200 )edit

"Sage is also a system for doing research, that includes comparing different implementations." To compare comparable implementations is of course OK. Comparing Maxima and PARI versions is OK. However, for instance, creating the impression that a 32-bit integer implementation will give the same results as an arbitrary precision implementation is not OK. Similarly here: To sell SciPy's very limited validity in hypergeometric_U() as 'less accurate' is not sufficient as long the impression persists that it is comparable to PARI or Maxima. (contd.)

Peter Luschny ( 2014-09-23 14:24:46 +0200 )edit

"A lot of work went into picking suitable defaults." Did you notice that the default in our case, i.e. not using the option "scipy", does not work? "Why shouldn't we expose the other implementations, too?" If you offer something which will lead a normal user with a certain trust in Sage inevitably into making errors than this is a disservice. Better save your time and at the same time preserve your customers from mistakes.

Peter Luschny ( 2014-09-23 14:25:11 +0200 )edit

add a comment

Why does Sage offer unusable SciPy?

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Why does Sage offer unusable SciPy? edit

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Why does Sage offer unusable SciPy?