I haven't had time to investigate your program and
findings, but I have no reason to suspect that your
results are flawed. I am somewhat puzzled by them.
After all the work you did, I hesitate to suggest
that I would rather have gfortran produce correct
results slowly (ie, n2 - n1 + 1 calls to mpfr_{jn,yn})
than to get possibly wrong results fast via
recurrence algorithm. If you want to change back
to your original algorithm that's fine with me;
otherwise, the current patch can be committed and
I'll see if I can find time to poke at it.