Fixed Bessel I implementation for Scipy

Aims to address
http://scipy.org/scipy/scipy/ticket/849
http://scipy.org/scipy/scipy/ticket/503

It's a reimplementation of Iv, based on code converted to C from Boost. Licenses should be OK.

I haven't checked if the netlib.org/specfun (not the Specfun in Scipy!) codes would be better than this. Do they have license issues?

[njs, 2009-02-26 03:48:50]

Read through, didn't see any issues.

[charlesr.harris, 2009-02-26 04:13:12]

Do you want to leave the fixed up function in cephes? If we are going to have a
verified set of functions it might be better to move it to a separate directory.

http://codereview.appspot.com/20078/diff/1/5
File scipy/special/cephes/scipy_iv.c (right):

http://codereview.appspot.com/20078/diff/1/5#newcode96
Line 96: v = -v;		/* symmetry */
Trailing comments tend to be fragile. I'm not sure this one is needed, but
prefer inline comments where possible.

http://codereview.appspot.com/20078/diff/1/5#newcode107
Line 107: if (v != 2.0 * floor(v / 2.0))
This is a style question. I like all the statements in if statements to be
bracketed as it makes them stand out more. Folks differ about this, but we
should settle on a style.

http://codereview.appspot.com/20078/diff/1/5#newcode118
Line 118: } else
When one branch of an if statement is bracketed, the other should be also, even
when it is only a single statement.

http://codereview.appspot.com/20078/diff/1/5#newcode124
Line 124: /* Uniform asymptotic expansion for large orders.
Multiline comments like this:

/*
 * blah
 * blah
 */

I'm more or less following linux here, with modifications for python.

http://codereview.appspot.com/20078/diff/1/5#newcode130
Line 130: } else {
Another one way or the other thingie. I prefer

}
else {

But in any case we should settle on a style.

http://codereview.appspot.com/20078/diff/1/5#newcode159
Line 159: factor = (mu - (2 * k - 1) * (2 * k - 1)) / (8 * x) / k;
Do we want to put spaces around * and / ?

http://codereview.appspot.com/20078/diff/1/5#newcode175
Line 175: * Computed with:
I wonder if we shouldn't consider sympy or something with arbitrary precision
for these sort of things so that they can be extended beyond double precision.
Not complaining, just suggesting.

http://codereview.appspot.com/20078/diff/1/5#newcode398
Line 398: /* |x| <= |v|, CF1_ik converges rapidly
See above on multiline comments.

http://codereview.appspot.com/20078/diff/1/5#newcode445
Line 445: BOOST_ASSERT(fabs(x) > 1);
BOOST_ASSERT is nice, but boost specific. Is it defined somewhere?

http://codereview.appspot.com/20078/diff/1/5#newcode490
Line 490: need_i = 1,
Ummm.. Ok. Why not hex as these are used as flags? Maybe a comment to that
effect also. Static const int need_i = 0x1 might be an alternative.

[charlesr.harris, 2009-02-26 04:16:27]

http://codereview.appspot.com/20078/diff/1/2
File scipy/special/cephes/const.c (right):

http://codereview.appspot.com/20078/diff/1/2#newcode64
Line 64: double EULER  = 0.57721566490153286061;        /* Euler constant */
More precision might be helpful if we later want to implement extended/quad
precision variants. Thirty two digits might be a good precision if you can find
it.

It might also be useful to move the constants up to a higher level so they can
be used in more that cephes.

[charlesr.harris, 2009-02-26 04:19:35]

http://codereview.appspot.com/20078/diff/1/3
File scipy/special/cephes/hyperg.c (left):

http://codereview.appspot.com/20078/diff/1/3#oldcode87
Line 87: double hyperg( a, b, x)
The names hyperg and hy1f1p are a bit cryptic. I don't know what we want to do
about that, though.

[stefanv, 2009-02-26 22:15:46]

Hey all,

Charles, a question for you.  Pauli, wondering about how to check for nan's. 
Good patch!

Cheers
Stéfan

http://codereview.appspot.com/20078/diff/1/3
File scipy/special/cephes/hyperg.c (right):

http://codereview.appspot.com/20078/diff/1/3#newcode207
Line 207: *err = 1.0;
Looks like the whole file needs some re-indentation.  According to the Python 3k
standard that we follow, 4 spaces per indent.

http://codereview.appspot.com/20078/diff/1/3#newcode288
Line 288: if (acanc != acanc)
I saw that the GNU C manual says:

x == x is false if the value of x is NaN. You can use this to test whether a
value is NaN or not, but the recommended way to test for NaN is with the isnan
function.

http://codereview.appspot.com/20078/diff/1/5
File scipy/special/cephes/scipy_iv.c (right):

http://codereview.appspot.com/20078/diff/1/5#newcode96
Line 96: v = -v;		/* symmetry */
On 2009/02/26 04:13:12, charlesr.harris wrote:
> Trailing comments tend to be fragile. I'm not sure this one is needed, but
> prefer inline comments where possible.

Hey Charles, out of curiosity, what do you mean by "fragile"?

http://codereview.appspot.com/20078/diff/1/6
File scipy/special/tests/test_basic.py (right):

http://codereview.appspot.com/20078/diff/1/6#newcode1648
Line 1648: for v, z in self._cephes_vs_amos_points():
I like this refactoring!

[pav, 2009-03-01 19:42:59]

Charles wrote:
> Do you want to leave the fixed up function in cephes? If 
> we are going to have a verified set of functions it might
> be better to move it to a separate directory.

It uses the error handling, constants, and gamma function from Cephes, so we'd
have to copy those to the core library, too. I changed the file name to indicate
that it's a significantly revised version, though.

I think the 'core' library should be introduced separately later, as there are
build issues etc. to think about.

http://codereview.appspot.com/20078/diff/1/3
File scipy/special/cephes/hyperg.c (left):

http://codereview.appspot.com/20078/diff/1/3#oldcode87
Line 87: double hyperg( a, b, x)
On 2009/02/26 04:19:35, charlesr.harris wrote:
> The names hyperg and hy1f1p are a bit cryptic. I don't know what we want to do
> about that, though.

These are the Cephes names. I'd rather not go changing them.

http://codereview.appspot.com/20078/diff/1/3
File scipy/special/cephes/hyperg.c (right):

http://codereview.appspot.com/20078/diff/1/3#newcode207
Line 207: *err = 1.0;
On 2009/02/26 22:15:46, Stefan wrote:
> Looks like the whole file needs some re-indentation.  According to the Python
3k
> standard that we follow, 4 spaces per indent.

Yes, the Cephes indentation is horrible. But I wanted to make any
re-indentations as a separate SVN commit, so that it's more clear what has been
changed and what hasn't.

http://codereview.appspot.com/20078/diff/1/3#newcode288
Line 288: if (acanc != acanc)
> but the recommended way to test for NaN is with the isnan

isnan is a C99 extension, so no other alternatives than testing for (x == x).

http://codereview.appspot.com/20078/diff/1/5
File scipy/special/cephes/scipy_iv.c (right):

http://codereview.appspot.com/20078/diff/1/5#newcode96
Line 96: v = -v;		/* symmetry */
On 2009/02/26 22:15:46, Stefan wrote:
> On 2009/02/26 04:13:12, charlesr.harris wrote:
> > Trailing comments tend to be fragile. I'm not sure this one is needed, but
> > prefer inline comments where possible.
> 
> Hey Charles, out of curiosity, what do you mean by "fragile"?

I'd rather not spend time on changing these to inline comments, unless there's a
good reason to do it.

http://codereview.appspot.com/20078/diff/1/5#newcode107
Line 107: if (v != 2.0 * floor(v / 2.0))
On 2009/02/26 04:13:12, charlesr.harris wrote:
> This is a style question. I like all the statements in if statements to be
> bracketed as it makes them stand out more. Folks differ about this, but we
> should settle on a style. 

Done.

http://codereview.appspot.com/20078/diff/1/5#newcode118
Line 118: } else
On 2009/02/26 04:13:12, charlesr.harris wrote:
> When one branch of an if statement is bracketed, the other should be also,
even
> when it is only a single statement.

Done.

http://codereview.appspot.com/20078/diff/1/5#newcode124
Line 124: /* Uniform asymptotic expansion for large orders.
On 2009/02/26 04:13:12, charlesr.harris wrote:
> Multiline comments like this:
> 
> /*
>  * blah
>  * blah
>  */
> 
> I'm more or less following linux here, with modifications for python.

Done.

http://codereview.appspot.com/20078/diff/1/5#newcode130
Line 130: } else {
On 2009/02/26 04:13:12, charlesr.harris wrote:
> Another one way or the other thingie. I prefer
> 
> }
> else {
> 
> But in any case we should settle on a style.

I think that's ugly, but changed nevertheless.

http://codereview.appspot.com/20078/diff/1/5#newcode159
Line 159: factor = (mu - (2 * k - 1) * (2 * k - 1)) / (8 * x) / k;
On 2009/02/26 04:13:12, charlesr.harris wrote:
> Do we want to put spaces around * and / ?

This is what 'indent -kr' gave. Certainly I'd put a space between (...) / (...),
but possibly not in 2 * k.

http://codereview.appspot.com/20078/diff/1/5#newcode175
Line 175: * Computed with:
On 2009/02/26 04:13:12, charlesr.harris wrote:
> I wonder if we shouldn't consider sympy or something with arbitrary precision
> for these sort of things so that they can be extended beyond double precision.
> Not complaining, just suggesting.

I think we could then just use Numpy's longdoubles. But it's probably not a big
trouble to rewrite this in Sympy/mpmath.

http://codereview.appspot.com/20078/diff/1/5#newcode398
Line 398: /* |x| <= |v|, CF1_ik converges rapidly
On 2009/02/26 04:13:12, charlesr.harris wrote:
> See above on multiline comments.

Done.

http://codereview.appspot.com/20078/diff/1/5#newcode445
Line 445: BOOST_ASSERT(fabs(x) > 1);
On 2009/02/26 04:13:12, charlesr.harris wrote:
> BOOST_ASSERT is nice, but boost specific. Is it defined somewhere?

Defined in this file, left in place since it's boost code.

http://codereview.appspot.com/20078/diff/1/5#newcode490
Line 490: need_i = 1,
On 2009/02/26 04:13:12, charlesr.harris wrote:
> Ummm.. Ok. Why not hex as these are used as flags? Maybe a comment to that
> effect also. Static const int need_i = 0x1 might be an alternative.

Done.