/*
 *      l o g
 */

#define __LIBRARY__

#include <math.h>
#include <stdlib.h>
#include <errno.h>

/* These coefficients come from Abramowitz & Stegun 4.1.27 which is  
*  the series used for the original algorithm. Error < E-19 
*  S E Peterson 3/92 */


static double logcoeff[] =
       {
		2.0,						/* 2/1 */
		6.6666666666666666667E-1,	/* 2/3 */
		4.0000000000000000000E-1,   /* 2/5 */
		2.8571428571428571429E-1,	/* 2/7 */
		2.2222222222222222222E-1,	/* etc. */
		1.8181818181818181818E-1,
		1.5384615384615384615E-1,
		1.3333333333333333333E-1,
		1.1764705882352941176E-1,
		1.0526315789473684210E-1,
		9.5238095238095238095E-2,
		8.6956521739130434782E-2,
		8.0000000000000000000E-2,
		7.4074074074074074074E-2,
		6.8965517241379310345E-2,
		6.4516129032258064516E-2,
		6.0606060606060606061E-2,
		5.7142857142857142857E-2,
		5.4054054054054054054E-2
       };

double log  _LIB_F1_(double,x)
{
 double fract;
 double ratio;
 int exponent,ncf;

 if(x <= 0.0)
   {
    errno = EDOM;
    return -HUGE_VAL;
   }

 fract = frexp(x, &exponent);
 ratio = (fract - 1.0) / (fract + 1.0);

if (fract < 0.58)			/* convergence is very slow when fract ==0.5 */
	ncf = 19;				/* so reduce order of polynomial as fract    */ 
else if (fract <0.70)		/* goes toward 1.0							 */
	ncf = 16;
else if (fract <0.81)
	ncf = 13;
else if (fract <0.9)
	ncf = 10;
else
	ncf = 8;	
				
return(exponent * M_LN2 + ratio * _poly(ratio * ratio, logcoeff, ncf));
}