math: bessel cleanup (jn.c and jnf.c)
authorSzabolcs Nagy <nsz@port70.net>
Tue, 1 Jan 2013 21:20:45 +0000 (22:20 +0100)
committerSzabolcs Nagy <nsz@port70.net>
Tue, 1 Jan 2013 21:20:45 +0000 (22:20 +0100)
commit5652d70054daf3c2c9b6d475fdf9d24a940e51aa
treebe9e8823682335c7fda681166948d20e61031b23
parent5bb6b24952e3f95ede42b60ac64a33ac34b8e272
math: bessel cleanup (jn.c and jnf.c)

both jn and yn functions had integer overflow issues for large
and small n

to handle these issues nm1 (== |n|-1) is used instead of n and -n
in the code and some loops are changed to make sure the iteration
counter does not overflow

(another solution could be to use larger integer type or even double
but that has more size and runtime cost, on x87 loading int64_t or
even uint32_t into an fpu register is more than two times slower than
loading int32_t, and using double for n slows down iteration logic)

yn(-1,0) now returns inf

posix2008 specifies that on overflow and at +-0 all y0,y1,yn functions
return -inf, this is not consistent with math when n<0 odd integer in yn
(eg. when x->0, yn(-1,x)->inf, but historically yn(-1,0) seems to be
special cased and returned -inf)

some threshold values in jnf and ynf were fixed that seems to be
incorrectly copy-pasted from the double version
src/math/jn.c
src/math/jnf.c