[SciPy-dev] Bessel functions from Boost
Sun Feb 8 22:46:34 CST 2009
Michael Abshoff wrote:
> Since I assume some people around here are interested in arbitrary
> precisions and after looking some more at the documentation it seems
> that that library only supports this via using an NTL type which in turn
> uses GMP. NTL itself is GPLed, GMP is LGPL, so either one does not fit
> the licensing requirements of Scipy.
> David mentioned to write a library from scratch and I also assume that
> you want arbitrary precisions.
Actually, I did not think about arbitrary precision at all. Not that it
would not be nice, but I simply know nothing about the topic :)
I am far from being entirely convinced by my own suggestion of writing
from scratch: it is a huge amount of work, and the possibilities to get
it wrong are numerous. But I am wondering whether we have a choice: the
only existing code which is licence compatible with scipy is the one we
use now (cephes, specfun, etc...), and the code is not great (cephes is
not even ANSI C, for example), except for toms maybe.
My suggestion is based on how R is doing things; I tend to consider R as
a reference, or at least a pretty good baseline when precision and
accuracy are concerned. R does not use Cephes nor specfun. We of course
can't use R code, but we could at least take inspiration of their
references (e.g. citations for implementation), and use R for testing.
I also think one of the problem of scipy.special is its size - there are
so many functions, with little to no testing. So maybe we could make up
a list of a small subset of functions which are 'essential', and focus
on them (we would of course keep the current code). I don't know whether
such a small subset exists. For information, nmath (the core maths
routines in R) is ~7500 LOC according to sloccount, and they have ~ 100
functions, maybe that would be a good subset ?
I also have no idea how to test those functions: when someones says
their function is precise up to 1e-6, is it against some theoretical
values which are computable (like gamma(0.5) = euler constant), from
theoretical consideration on the implementation ? Pauli talked about
point tests, but it sounds hard to get the right points for testing ?
More information about the Scipy-dev