/*
 *      l o g f 
 */

#define __LIBRARY__

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

float logf _LIB_F1_(float,x)

/* This is truncated version of series from Abramowitz & Stegun 
* 4.1.27 as used in  double log(x); error should be < E-9 */

{
	static float coef[] =  {2.0,				/* 2/1 */
							6.666666667E-1,		/* 2/3 */
							4.000000000E-1,		/* 2/5 */
							2.857142857E-1,		/* 2/7 */
							2.222222222E-1,		/* etc. */
							1.818181818E-1,
							1.538461538E-1,
							1.333333333E-1,
							1.176470588E-1};
							
	int ncf, exponent;
	float ans, fract, *ptr, ratio, r2;

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

  	fract = (float)frexp((double)x, &exponent);
	ratio = (fract - 1.0) / (fract + 1.0);
	r2 = ratio * ratio;
	
	if (fract < 0.6)			/* No. of terms required for worst case */
		ncf = 9;				/* fract = 0.5; fewer terms required for */
	else if (fract < 0.8)		/* same accuracy with larger values of fract */
		ncf = 7;
	else
		ncf = 5;	
	
	ptr = &coef[0];
	ptr += ncf;
	ans = *--ptr; 	
	while (--ncf > 0){
		ans = *--ptr +(r2 * ans);  
	}
	return (exponent * M_LN2 + ratio * ans);
}