Bug in libc6

[Lei Yin <leiyin at 263 dot net>]

The attached file has the different results in libc5 and libc6.
The results in libc5 is correct. We checked the program step by
step and found that the bessel function (y0) in libc6 carries
out the different (in addition, the wrong) results compared 
with it derived by libc5.

Please help us to check it. Thanks.

Yours

---
Lei Yin                          __________________________________________
                                 )  M r. L e i   Y i n =================== )_
   ,-~~-.___.                   / State Key Laboratory of Millimeter Wave / /
  / |  '     \            ()   / Radio Eng. Dept. , Southeast University / /
 (  )         0         O     /........................................./ /
  _/-, ,----'        o      / E-mail: Lei Yin <leiyin@263.net>        / /
     ====           //      /     Tel: +86-25-3793276     (Office)    / /
    /  \-'~;    /~~~(O)    /           +86-550-3025280    (Home)     / /
 =(  _____| (_________ |  (_________________________________________(_/

#include <stdio.h>
#include <math.h>

#define PI 3.1415926535898 
#define N 10

int main(void)
{
	int n;
	double Esr,Esi;
	double a,phi,R;
	double bjn,bynr,byni;
	double byn1r,byn1i,byn2r,byn2i;
	double rcsdlmt;
	double temp,sign;
	
	R=10.0; 
	//a=1/2.0/PI;
	a=1.0/2.0/PI;

	printf("#N=%d,R=%f,a=%f\n",N,R,a);
	for (phi=0; phi<=360.0; phi+=1.0)
	{
		byn1r=j0(2.0*PI*R);
		byn1i=-y0(2.0*PI*R);
		byn2r=j0(2.0*PI*a);
		byn2i=-y0(2.0*PI*a);
		temp=byn2r*byn2r+byn2i*byn2i;
		bynr=(byn1r*byn2r+byn1i*byn2i)/temp;
		byni=(byn2r*byn1i-byn1r*byn2i)/temp;
		Esr=-j0(2.0*PI*a)*bynr;
		Esi=-j0(2.0*PI*a)*byni;
		sign=1.0;
		for (n=1; n<N; n+=2)
		{
			byn1r=jn(n,2.0*PI*R);
			byn1i=-yn(n,2.0*PI*R);
			byn2r=jn(n,2.0*PI*a);
			byn2i=-yn(n,2.0*PI*a);
		        temp=byn2r*byn2r+byn2i*byn2i;
			bynr=(byn1r*byn2r+byn1i*byn2i)/temp;
			byni=(byn2r*byn1i-byn1r*byn2i)/temp;
			temp=sign*2.0*jn(n,2.0*PI*a)*cos(n*phi*PI/180.0);
			Esr+=-byni*temp;
			Esi+=bynr*temp;
			sign*=-1.0;
		}
		sign=1.0;
		for (n=2; n<N; n+=2)
		{
			byn1r=jn(n,2.0*PI*R);
			byn1i=-yn(n,2.0*PI*R);
			byn2r=jn(n,2.0*PI*a);
			byn2i=-yn(n,2.0*PI*a);
		        temp=byn2r*byn2r+byn2i*byn2i;
			bynr=(byn1r*byn2r+byn1i*byn2i)/temp;
			byni=(byn2r*byn1i-byn1r*byn2i)/temp;
			temp=sign*2.0*jn(n,2.0*PI*a)*cos(n*phi*PI/180.0);
			Esr+=bynr*temp;
			Esi+=byni*temp;
			sign*=-1.0;
		}
		rcsdlmt=sqrt(2.0*PI*R*(Esr*Esr+Esi*Esi));
		printf("%f  %f \n",phi,rcsdlmt);
	}

	return 0;

}