log.h

/*
 * log.h
 * The basic idea is to exploit Pade polynomials.
 * A lot of ideas were inspired by the cephes math library (by Stephen L. Moshier
 * moshier@na-net.ornl.gov) as well as actual code. 
 * The Cephes library can be found here:  http://www.netlib.org/cephes/
 * 
 *  Created on: Jun 23, 2012
 *      Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente
 */

/* 
 * VDT is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef LOG_H_
#define LOG_H_

#include "vdtcore_common.h"
#include <limits>

namespace vdt{

// local namespace for the constants/functions which are necessary only here
namespace details{

const double LOG_UPPER_LIMIT = 1e307;
const double LOG_LOWER_LIMIT = 0;

const double SQRTH = 0.70710678118654752440;

inline double get_log_px(const double x){
	const double PX1log = 1.01875663804580931796E-4;
	const double PX2log = 4.97494994976747001425E-1;
	const double PX3log = 4.70579119878881725854E0;
	const double PX4log = 1.44989225341610930846E1;
	const double PX5log = 1.79368678507819816313E1;
	const double PX6log = 7.70838733755885391666E0;

	double px = PX1log;
	px *= x;
	px += PX2log;
	px *= x;
	px += PX3log;
	px *= x;
	px += PX4log;
	px *= x;
	px += PX5log;
	px *= x;
	px += PX6log;
	return px;

}

inline double get_log_qx(const double x){
	const double QX1log = 1.12873587189167450590E1;
	const double QX2log = 4.52279145837532221105E1;
	const double QX3log = 8.29875266912776603211E1;
	const double QX4log = 7.11544750618563894466E1;
	const double QX5log = 2.31251620126765340583E1;

	double qx = x;
	qx += QX1log;
	qx *=x;
	qx += QX2log;
	qx *=x;
	qx += QX3log;
	qx *=x;
	qx += QX4log;
	qx *=x;
	qx += QX5log;
	return qx;
}

}

// Log double precision --------------------------------------------------------
inline double fast_log(double x){

	const double original_x = x;

	/* separate mantissa from exponent */
	double fe;
	x = details::getMantExponent(x,fe);

	// blending
	x > details::SQRTH? fe+=1. : x+=x ;
	x -= 1.0;

	/* rational form */
	double px =  details::get_log_px(x);

	//for the final formula
	const double x2 = x*x;
	px *= x;
	px *= x2;

	const double qx = details::get_log_qx(x);

	double res = px / qx ;

	res -= fe * 2.121944400546905827679e-4;
	res -= 0.5 * x2  ;

	res = x + res;
	res += fe * 0.693359375;

	if (original_x > details::LOG_UPPER_LIMIT)
		res = std::numeric_limits<double>::infinity();
	if (original_x < details::LOG_LOWER_LIMIT) // THIS IS NAN!
		res =  - std::numeric_limits<double>::quiet_NaN();

	return res;

}

// Log single precision --------------------------------------------------------



namespace details{

const float LOGF_UPPER_LIMIT = MAXNUMF;
const float LOGF_LOWER_LIMIT = 0;

const float PX1logf = 7.0376836292E-2f;
const float PX2logf	= -1.1514610310E-1f;
const float PX3logf	= 1.1676998740E-1f;
const float PX4logf	= -1.2420140846E-1f;
const float PX5logf	= 1.4249322787E-1f;
const float PX6logf	= -1.6668057665E-1f;
const float PX7logf	= 2.0000714765E-1f;
const float PX8logf	= -2.4999993993E-1f;
const float PX9logf	= 3.3333331174E-1f;

inline float get_log_poly(const float x){
	float y = x*PX1logf;
	y += PX2logf;
	y *= x;
	y += PX3logf;
	y *= x;
	y += PX4logf;
	y *= x;
	y += PX5logf;
	y *= x;
	y += PX6logf;
	y *= x;
	y += PX7logf;
	y *= x;
	y += PX8logf;
	y *= x;
	y += PX9logf;
	return y;
}

const float SQRTHF = 0.707106781186547524f;

}

// Log single precision --------------------------------------------------------
inline float fast_logf( float x ) {

	const float original_x = x;

	float fe;
	x = details::getMantExponentf( x, fe);

	x > details::SQRTHF? fe+=1.f : x+=x ;
	x -= 1.0f;

	const float x2 = x*x;

	float res = details::get_log_poly(x);
	res *= x2*x;

	res += -2.12194440e-4f * fe;
	res +=  -0.5f * x2;

	res= x + res;

	res += 0.693359375f * fe;

	if (original_x > details::LOGF_UPPER_LIMIT)
		res = std::numeric_limits<float>::infinity();
	if (original_x < details::LOGF_LOWER_LIMIT)
		res = -std::numeric_limits<float>::quiet_NaN();

	return res;
}


//------------------------------------------------------------------------------

void logv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
void fast_logv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
void logfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
void fast_logfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);

} //vdt namespace

#endif /* LOG_H_ */