17  *	Given x, find r and integer k such that
19  *               x = k*ln2 + r,  |r| <= 0.5*ln2 ~ 0.34658
24  *   2. Approximating expm1(r) by a special rational function on
27  *	    r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
29  *	    r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
31  *	    R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
33  *		     = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
37  *	by 2**-61. In other words,
38  *	    R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
39  *	where 	Q1  =  -1.6666666666666567384E-2,
47  *	    | 1.0+Q1*z+...+Q5*z   -  R1(z) | <= 2
52  *			       2     3
53  *			      r     r    [ 3 - (R1 + R1*r/2)  ]
55  *		              2     2    [ 6 - r*(3 - R1*r/2) ]
63  *		        (      2                                    2 )
64  *		        ({  ( r    [ R1 -  (3 - R1*r/2) ]  )  }    r  )
66  *	                ({  ( 2    [ 6 - r*(3 - R1*r/2) ]  )  }    2  )
72  *	   expm1(x) = either 2^k*[expm1(r)+1] - 1
73  *		    = or     2^k*[expm1(r) + (1-2^-k)]
76  *	     to Qi*2^i, and replace z by (x^2)/2.
79  *	  (ii)  if k=0, return r-E
80  *	  (iii) if k=-1, return 0.5*(r-E)-0.5
81  *        (iv)	if k=1 if r < -0.25, return 2*((r+0.5)- E)
83  *	  (v)   if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
84  *	  (vi)  if k <= 20, return 2^k((1-2^-k)-(E-r)), else
85  *	  (vii) return 2^k(1-((E+2^-k)-r))
119 /* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */
132 	int32_t k,xsb;  in expm1()  local
161 		    {hi = x - ln2_hi; lo =  ln2_lo;  k =  1;}  in expm1()
163 		    {hi = x + ln2_hi; lo = -ln2_lo;  k = -1;}  in expm1()
165 		k  = invln2*x+((xsb==0)?0.5:-0.5);  in expm1()
166 		t  = k;  in expm1()
173 	else if(hx < 0x3c900000) {  	/* when |x|<2**-54, return x */  in expm1()
177 	else k = 0;  in expm1()
185 	if(k==0) return x - (x*e-hxs);		/* c is 0 */  in expm1()
187 	    INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20,0);	/* 2^k */  in expm1()
190 	    if(k== -1) return 0.5*(x-e)-0.5;  in expm1()
191 	    if(k==1) {  in expm1()
195 	    if (k <= -2 || k>56) {   /* suffice to return exp(x)-1 */  in expm1()
197 		if (k == 1024) y = y*2.0*0x1p1023;  in expm1()
202 	    if(k<20) {  in expm1()
203 	        SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k));  /* t=1-2^-k */  in expm1()
207 		SET_HIGH_WORD(t,((0x3ff-k)<<20));	/* 2^-k */  in expm1()