# Why does Sage offer unusable SciPy?

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

edit retag close merge delete

Sort by » oldest newest most voted

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.

more

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.

( 2014-09-20 06:13:40 -0600 )edit
2

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...

( 2014-09-21 12:13:01 -0600 )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.)

( 2014-09-23 07:24:46 -0600 )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.

( 2014-09-23 07:25:11 -0600 )edit