17  *   1. Argument Reduction: find k and f such that 
18  *			x = 2^k * (1+f), 
19  *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
21  *   2. Approximation of log(1+f).
22  *	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
23  *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
24  *	     	 = 2s + s*R
27  *	of this polynomial approximation is bounded by 2**-58.45. In
29  *		        2      4      6      8      10      12      14
33  *	    |      2          14          |     -58.45
34  *	    | Lg1*s +...+Lg7*s    -  R(z) | <= 2 
36  *	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
39  *		log(1+f) = f - s*(f - R)	(if f is not too large)
40  *		log(1+f) = f - (hfsq - s*(hfsq+R)).	(better accuracy)
42  *	3. Finally,  log(x) = k*ln2 + log(1+f).  
43  *			    = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
87 	double hfsq,f,s,z,R,w,t1,t2,dk;  in log()  local
94 	if (hx < 0x00100000) {			/* x < 2**-1022  */  in log()
105 	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */  in log()
107 	f = x-1.0;  in log()
108 	if((0x000fffff&(2+hx))<3) {	/* -2**-20 <= f < 2**-20 */  in log()
109 	    if(f==zero) {  in log()
117 	    R = f*f*(0.5-0.33333333333333333*f);  in log()
118 	    if(k==0) return f-R; else {dk=(double)k;  in log()
119 	    	     return dk*ln2_hi-((R-dk*ln2_lo)-f);}  in log()
121  	s = f/(2.0+f);   in log()
132 	    hfsq=0.5*f*f;  in log()
133 	    if(k==0) return f-(hfsq-s*(hfsq+R)); else  in log()
134 		     return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);  in log()
136 	    if(k==0) return f-s*(f-R); else  in log()
137 		     return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);  in log()