Bug in besselk and similar functions

[Michael Stoll <Michael_Stoll@math.uni-bonn.de>]

In a 2.0.x PARI version, try (with the default realprecision of 28)
> besselk(0,15.999)
and then
> besselk(0,16.0)
> besselk(0,32.999)
> nesselk(0,33.0)
The first and fourth work OK, but the second and third give (at least
on my machine) a stack overflow, regardless of how big the stack is.
It did work in PARI 1.39, so I did a comparison and found that probably
the only essential difference is (line 48 in the new version)
--- old ---
  lbin=10-((l-2)<<TWOPOTBITS_IN_LONG);
*** new ***
  lbin = 5 - bit_accuracy(l);
-----------
In the body of kbessel(), there is a construct of the following type:
    for(;;)
    {
      [...]
      for (k=1; ; k++)
      {
	[...]
	if (gexpo(p1)-gexpo(f) <= lbin) break;
	avma=av1;
      }
      [...]
      if (expo(subrr(q,r)) <= lbin) break;
    }
If x = 15.999, the inner loop is executed 25,26,27,28,29,31,30 times,
then the outer break statement is triggered. If x = 16.0, the sequence
is 25,26,27,28,29,31,30,1,1,1,... , and the required accuracy is never
reached. Increasing the 5 above to 7 seems to work, at least in the cases
I've tested. What was the reason for changing 10 into 5?
The same holds for hyperu(). If you use besselk(0,x,1), you also get
an endless loop for x = 16.0, but no stack overflow (which is probably
due to a better memory organisation in hyperu()).
I'd suggest to use the old 10, unless somebody gives a good reason to
take a smaller number (and a proof that it works!). Since there are no
comments (as usual), it is difficult to know what lbin is intended to
do.

Cheers,
        Michael