EIC code displayed by LXR master/​include/​vdt/​log.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-02-21 10:16:42

 /*
  * 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_ */

This page was automatically generated by the 2.3.7 LXR engine. The LXR team