hx = __HI(x); lx = __LO(x);
u.d = x; hx = __HI(u); lx = __LO(u);
	double __ieee754_lgamma_r(double x, int *signgamp)#else	double __ieee754_lgamma_r(x,signgamp)	double x; int *signgamp;#endif{	double t,y,z,nadj,p,p1,p2,p3,q,r,w;	int i,hx,lx,ix;	hx = __HI(x);	lx = __LO(x);    /* purge off +-inf, NaN, +-0, and negative arguments */	*signgamp = 1;	ix = hx&0x7fffffff;	if(ix>=0x7ff00000) return x*x;	if((ix|lx)==0) return one/zero;	if(ix<0x3b900000) {	/* |x|<2**-70, return -log(|x|) */	    if(hx<0) {	        *signgamp = -1;	        return -__ieee754_log(-x);	    } else return -__ieee754_log(x);	}	if(hx<0) {	    if(ix>=0x43300000) 	/* |x|>=2**52, must be -integer */		return one/zero;	    t = sin_pi(x);	    if(t==zero) return one/zero; /* -integer */	    nadj = __ieee754_log(pi/fd_fabs(t*x));	    if(t<zero) *signgamp = -1;	    x = -x;	}    /* purge off 1 and 2 */	if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;    /* for x < 2.0 */	else if(ix<0x40000000) {	    if(ix<=0x3feccccc) { 	/* lgamma(x) = lgamma(x+1)-log(x) */		r = -__ieee754_log(x);		if(ix>=0x3FE76944) {y = one-x; i= 0;}		else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}	  	else {y = x; i=2;}	    } else {	  	r = zero;	        if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */	        else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */		else {y=x-one;i=2;}	    }	    switch(i) {	      case 0:		z = y*y;		p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));		p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));		p  = y*p1+p2;		r  += (p-0.5*y); break;	      case 1:		z = y*y;		w = z*y;		p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12)));	/* parallel comp */		p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));		p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));		p  = z*p1-(tt-w*(p2+y*p3));		r += (tf + p); break;	      case 2:			p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));		p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));		r += (-0.5*y + p1/p2);	    }	}	else if(ix<0x40200000) { 			/* x < 8.0 */	    i = (int)x;	    t = zero;	    y = x-(double)i;	    p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));	    q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));	    r = half*y+p/q;	    z = one;	/* lgamma(1+s) = log(s) + lgamma(s) */	    switch(i) {	    case 7: z *= (y+6.0);	/* FALLTHRU */	    case 6: z *= (y+5.0);	/* FALLTHRU */	    case 5: z *= (y+4.0);	/* FALLTHRU */	    case 4: z *= (y+3.0);	/* FALLTHRU */	    case 3: z *= (y+2.0);	/* FALLTHRU */		    r += __ieee754_log(z); break;	    }    /* 8.0 <= x < 2**58 */	} else if (ix < 0x43900000) {	    t = __ieee754_log(x);	    z = one/x;	    y = z*z;	    w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));	    r = (x-half)*(t-one)+w;	} else     /* 2**58 <= x <= inf */	    r =  x*(__ieee754_log(x)-one);	if(hx<0) r = nadj - r;	return r;}