#include "grib_int.h"

void
confp3(double pval, int *kexp, int *kmant, int kbits, int kround)
{
  /*

    Purpose:
    --------

    Convert floating point number from machine
    representation to GRIB representation.

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

       pval    - Floating point number to be converted.
       kbits   - Number of bits in computer word.
       kround  - Conversion type.
                 0 , Closest number in GRIB format less than
                     original number.
                 1 , Closest number in GRIB format to the
                     original number (equal to, greater than or
                     less than original number).

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

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

    Method:
    -------

    Floating point number represented as 8 bit signed
    exponent and 24 bit mantissa in integer values.

    Externals.
    ----------

    decfp2    - Decode from IBM floating point format.

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

    WMO Manual on Codes re GRIB representation.

    Comments:
    ---------

    Routine aborts if an invalid conversion type parameter
    is used or if a 24 bit mantissa is not produced.

    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 1.6 on NEC SX6
        - replace 1.0 / pow(16.0, (double)(iexp - 70)) by rpow16m70tab[iexp]
  */

  // extern int CGRIBEX_Debug;

  /* ----------------------------------------------------------------- */
  /*   Section 1 . Initialise                                          */
  /* ----------------------------------------------------------------- */

  // Check conversion type parameter.

  int iround = kround;
  if (iround != 0 && iround != 1)
    {
      Error("Invalid conversion type = %d", iround);

      // If not aborting, arbitrarily set rounding to 'up'.
      iround = 1;
    }

  /* ----------------------------------------------------------------- */
  /*   Section 2 . Convert value of zero.                              */
  /* ----------------------------------------------------------------- */

  if (fabs(pval) <= 0)
    {
      *kexp = 0;
      *kmant = 0;
      goto LABEL900;
    }

  /* ----------------------------------------------------------------- */
  /*   Section 3 . Convert other values.                               */
  /* ----------------------------------------------------------------- */
  {
    double zeps = (kbits != 32) ? 1.0e-12 : 1.0e-8;
    double zref = pval;

    // Sign of value.
    int isign = (zref >= 0.0) ? 0 : 128;
    zref = fabs(zref);

    // Exponent.
    int iexp = (int) (log(zref) / log(16.0) + 65.0 + zeps);

    // only ANSI C99 has log2
    // iexp = (int) (log2(zref) * 0.25 + 65.0 + zeps);

    if (iexp < 0) iexp = 0;
    if (iexp > 127) iexp = 127;

    // double rpowref = zref / pow(16.0, (double)(iexp - 70));
    double rpowref = ldexp(zref, 4 * -(iexp - 70));

    // Mantissa.
    if (iround == 0)
      {
        /*  Closest number in GRIB format less than original number. */
        /*  Truncate for positive numbers. */
        /*  Round up for negative numbers. */
        *kmant = (isign == 0) ? (int) rpowref : (int) lround(rpowref + 0.5);
      }
    else
      {
        /*  Closest number in GRIB format to the original number   */
        /*  (equal to, greater than or less than original number). */
        *kmant = (int) lround(rpowref);
      }

    /*  Check that mantissa value does not exceed 24 bits. */
    /*  If it does, adjust the exponent upwards and recalculate the mantissa. */
    /*  16777215 = 2**24 - 1 */
    if (*kmant > 16777215)
      {

      LABEL350:

        ++iexp;

        // Check for exponent overflow during adjustment
        if (iexp > 127)
          {
            Message("Exponent overflow");
            Message("Original number = %30.20f", pval);
            Message("Sign = %3d, Exponent = %3d, Mantissa = %12d", isign, iexp, *kmant);

            Error("Exponent overflow");

            // If not aborting, arbitrarily set value to zero
            Message("Value arbitrarily set to zero.");
            *kexp = 0;
            *kmant = 0;
            goto LABEL900;
          }

        rpowref = ldexp(zref, 4 * -(iexp - 70));

        if (iround == 0)
          {
            /*  Closest number in GRIB format less than original number. */
            /*  Truncate for positive numbers. */
            /*  Round up for negative numbers. */
            *kmant = (isign == 0) ? (int) rpowref : (int) lround(rpowref + 0.5);
          }
        else
          {
            /*  Closest number in GRIB format to the original number */
            /*  (equal to, greater or less than original number). */
            *kmant = (int) lround(rpowref);
          }

        // Repeat calculation (with modified exponent) if still have mantissa overflow.
        if (*kmant > 16777215) goto LABEL350;
      }

    // Add sign bit to exponent.
    *kexp = iexp + isign;
  }

  /* ----------------------------------------------------------------- */
  /*   Section 9. Return                                               */
  /* ----------------------------------------------------------------- */

LABEL900:
  /*
  if ( CGRIBEX_Debug )
    {
      double zval;

      Message("Conversion type parameter = %4d", kround);
      Message("Original number = %30.20f", pval);

      zval = decfp2(*kexp, *kmant);

      Message("Converted to      %30.20f", zval);
      Message("Sign = %3d, Exponent = %3d, Mantissa = %12d", isign, iexp, *kmant);
    }
  */
  return;
} /* confp3 */
