Ask Your Question
2

Why does Sage offer unusable SciPy?

asked 2014-09-17 06:32:32 -0500

Peter Luschny gravatar image

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

1 answer

Sort by » oldest newest most voted
1

answered 2014-09-18 03:54:56 -0500

updated 2014-09-20 12:21:20 -0500

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 more

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 gravatar imagePeter Luschny ( 2014-09-20 06:13:40 -0500 )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...

Volker Braun gravatar imageVolker Braun ( 2014-09-21 12:13:01 -0500 )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 gravatar imagePeter Luschny ( 2014-09-23 07:24:46 -0500 )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 gravatar imagePeter Luschny ( 2014-09-23 07:25:11 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2014-09-17 06:32:32 -0500

Seen: 141 times

Last updated: Sep 20 '14