source: trunk/HeuristicLab.ExtLibs/HeuristicLab.NativeInterpreter/0.1/NativeInterpreter-0.1/lib/vdt/log.h @ 16451

/*
 * 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;
}


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

} //vdt namespace

#endif /* LOG_H_ */