Subversion Repositories fp_log

/*******************************************************************************************

 *

 * Double-Precision ICSILog algorithm implementation

 * written by Nikolaos Alchiotis and Alexandros Stamatakis

 *

 * The Exelixis Lab

 * Bioinformatics Unit (I12)

 * Department of Computer Science

 * Technical University of Munich

 *

 * Emails: alachiot@cs.tum.edu, stamatak@cs.tum.edu

 * WWW:   http://wwwkramer.in.tum.de/exelixis/

 *

 * This code is made available under GNU GPL

 *

 *

*******************************************************************************************/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

/********************************************/

// Union to Access the bits

typedef union

{

  double value;

  struct

  {

    unsigned int rght_p;

    unsigned int lft_p;

  }

    parts;

}

  ieee_double_shape_type;

/********************************************/

//Constant Value

double con_val;

// Variables used for error checking

double value_minus_inf,

       value_nan,

       value_inf,

       value_inf_inf;

// Prec shows how many MSBs of the mantissa field will be used to index the mantissa LUT

unsigned int prec = 12;

// Lut_entries is the number of entries of the mantissa LUT

int lut_entries;//=pow(2,prec);

//ShiftBits is used by the main function to perform the shifting operation

int shiftBits;//= 12 + 20 - prec;

/********************************************/

// Function to initialize the mantissa LUT . It is the same one that the official ICSILog uses.

static void DP_ICSILog_INIT(const unsigned precision, double* const pTable)

{

  /*

     step along table elements and x-axis positions

     (start with extra half increment, so the steps intersect at their midpoints.)

  */

  con_val = log(2);

  value_minus_inf = -1.0 / 0.0;

  value_nan = 0.0 / 0.0;

  value_inf = pow(10, 308);

  value_inf_inf = value_inf + value_inf;

  double oneToTwo = 1.0f + (1.0f / (double)( 1 << (precision + 1)));

  int i;

  for(i = 0;  i < (1 << precision);  ++i )

    {

      // make y-axis value for table element

      pTable[i] = logf(oneToTwo);

      oneToTwo += 1.0f / (double)(1 << precision);

    }

}

/********************************************/

// The DP_myICSILog main function.

static inline double DP_myICSILog(const double input , const int num_of_bits, register double * const t)

{

  register double result;

  if(input < 0.0)

    result = value_nan;

  else

    {

      if(input == 0.0)

        result = value_minus_inf;

      else

        {

          if(input > value_inf)

            result = value_inf_inf;

          else

            {

              ieee_double_shape_type model_input;

              model_input.value = input;

              register unsigned int left = model_input.parts.lft_p;

              register double m1 = ((left << 1) >> 21) - 1023.0;

              register int index = (left << 12) >> num_of_bits;

              result = (con_val * m1  + t[index]);

            }

        }

    }

  return result;

}

int main()

{

  double

    value = 123456789,

    *DP_myICSILog_LUT,

    result;

  /*********************************************************************/

  // How to Initialize the function

  prec = 12;

  lut_entries = pow(2,prec);

  shiftBits = 12 + 20 - prec;

  DP_myICSILog_LUT = (double*)malloc(sizeof(double) * lut_entries);

  DP_ICSILog_INIT(prec, DP_myICSILog_LUT);

  /********************************************************************/

  /********************************************************************/

  // How to Call the function

  result = DP_myICSILog(value, shiftBits, DP_myICSILog_LUT);

  printf("\n\nDP_myICSILog(%f) = %f \n",value, result);

  printf("\nLog(%f) = %f \n\n\n", value, log(value));

  return 0;

}