Anyway, I got held up and I think it's good to have new folks looking
into this.

It looks good except that you need to uglify k.

I looked at the GSL implementation, based on same reference, and their
loop is cleaner. What about porting that implementation here? Possible?

My implementation is also using one more term for P than for Q, which
is discouraged in GSL, according to their comments.

Ed, do you have any comment about this point?

Regarding the 2nd line, after my observations, usually term stops
contributing to P, Q before k > nu/2, so actually, an offset by one is
most likely without consequences.

GSL implementation is nevertheless more elegant.

Also, keep in mind that these series are asymptotic - they'll
eventually blow up.

You should watch the magnitude of sequential terms and bail out of
the loop if either term starts growing.

Ed

Yes, you are right. As nu increases, the multiplying '__term' gets
larger, and the result will lose accuracy.

This cannot be used for large nu.

I have not been able to figure out how to detect when it starts to
lose accuracy though, other than empirically.
GSL has no checks on terms and uses this method only under certain
conditions, on nu and x.