#include <math.h>
#include "grib_int.h"

double
decfp2(int kexp, int kmant)
{
  /*

    Purpose:
    --------

    Convert GRIB representation of a floating point
    number to machine representation.

    Input Parameters:
    -----------------

    kexp    - 8 Bit signed exponent.
    kmant   - 24 Bit mantissa.

    Output Parameters:
    ------------------

    Return value   - Floating point number represented
                     by kexp and kmant.

    Method:
    -------

    Floating point number represented as 8 bit exponent
    and 24 bit mantissa in integer values converted to
    machine floating point format.

    Externals:
    ----------

    None.

    Reference:
    ----------

    WMO Manual on Codes re GRIB representation.

    Comments:
    ---------

    Rewritten from DECFP, to conform to programming standards.
    Sign bit on 0 value now ignored, if present.
    If using 32 bit reals, check power of 16 is not so small as to
    cause overflows (underflows!); this causes warning to be given
    on Fujitsus.

    Author:
    -------

    John Hennessy   ECMWF   18.06.91

    Modifications:
    --------------

    Uwe Schulzweida   MPIfM   01/04/2001

     - Convert to C from EMOS library version 130

    Uwe Schulzweida   MPIfM   02/08/2002

     - speed up by factor 2 on NEC SX6
        - replace pow(2.0, -24.0) by constant POW_2_M24
        - replace pow(16.0, (double)(iexp - 64)) by pow16m64tab[iexp]
  */

  /* ----------------------------------------------------------------- */
  /*   Section 1 . Convert value of 0.0. Ignore sign bit.              */
  /* ----------------------------------------------------------------- */

  if ((kexp == 128) || (kexp == 0) || (kexp == 255)) return 0.0;

  /* ----------------------------------------------------------------- */
  /*   Section 2 . Convert other values.                               */
  /* ----------------------------------------------------------------- */

  //  Sign of value.

  int iexp = kexp;
  int isign = (iexp < 128) * 2 - 1;

  iexp -= iexp < 128 ? 0 : 128;

  //  Decode value.

  // double pval = isign * pow(2.0, -24.0) * kmant * pow(16.0, (double)(iexp - 64));

  iexp -= 64;

  double pval = ldexp(1.0, 4 * iexp) * isign * POW_2_M24 * kmant;

  /* ----------------------------------------------------------------- */
  /*   Section 9. Return to calling routine.                           */
  /* ----------------------------------------------------------------- */

  return pval;
}
