/**************************************************************************/
/*                                                                        */
/*                                 OCaml                                  */
/*                                                                        */
/*             Xavier Leroy, projet Cristal, INRIA Rocquencourt           */
/*                                                                        */
/*   Copyright 1996 Institut National de Recherche en Informatique et     */
/*     en Automatique.                                                    */
/*                                                                        */
/*   All rights reserved.  This file is distributed under the terms of    */
/*   the GNU Lesser General Public License version 2.1, with the          */
/*   special exception on linking described in the file LICENSE.          */
/*                                                                        */
/**************************************************************************/

#define CAML_INTERNALS

/* The interface of this file is in "caml/mlvalues.h" and "caml/alloc.h" */

/* Needed for uselocale */
#define _XOPEN_SOURCE 700

/* Needed for strtod_l */
#define _GNU_SOURCE

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include <limits.h>

#include "caml/alloc.h"
#include "caml/fail.h"
#include "caml/memory.h"
#include "caml/mlvalues.h"
#include "caml/misc.h"
#include "caml/reverse.h"
#include "caml/stacks.h"

#if defined(HAS_LOCALE) || defined(__MINGW32__)

#if defined(HAS_LOCALE_H) || defined(__MINGW32__)
#include <locale.h>
#endif

#if defined(HAS_XLOCALE_H)
#include <xlocale.h>
#endif

#if defined(_MSC_VER)
#ifndef locale_t
#define locale_t _locale_t
#endif
#ifndef freelocale
#define freelocale _free_locale
#endif
#ifndef strtod_l
#define strtod_l _strtod_l
#endif
#endif

#endif /* defined(HAS_LOCALE) */

#ifdef _MSC_VER
#include <float.h>
#ifndef isnan
#define isnan _isnan
#endif
#ifndef isfinite
#define isfinite _finite
#endif
#ifndef nextafter
#define nextafter _nextafter
#endif
#endif

#ifdef ARCH_ALIGN_DOUBLE

CAMLexport double caml_Double_val(value val)
{
  union { value v[2]; double d; } buffer;

  CAMLassert(sizeof(double) == 2 * sizeof(value));
  buffer.v[0] = Field(val, 0);
  buffer.v[1] = Field(val, 1);
  return buffer.d;
}

CAMLexport void caml_Store_double_val(value val, double dbl)
{
  union { value v[2]; double d; } buffer;

  CAMLassert(sizeof(double) == 2 * sizeof(value));
  buffer.d = dbl;
  Field(val, 0) = buffer.v[0];
  Field(val, 1) = buffer.v[1];
}

#endif

/*
 OCaml runtime itself doesn't call setlocale, i.e. it is using
 standard "C" locale by default, but it is possible that
 third-party code loaded into process does.
*/
#ifdef HAS_LOCALE
locale_t caml_locale = (locale_t)0;
#endif

#if defined(_MSC_VER) || defined(__MINGW32__)
/* there is no analogue to uselocale in MSVC so just set locale for thread */
#define USE_LOCALE setlocale(LC_NUMERIC,"C")
#define RESTORE_LOCALE do {} while(0)
#elif defined(HAS_LOCALE)
#define USE_LOCALE locale_t saved_locale = uselocale(caml_locale)
#define RESTORE_LOCALE uselocale(saved_locale)
#else
#define USE_LOCALE do {} while(0)
#define RESTORE_LOCALE do {} while(0)
#endif

void caml_init_locale(void)
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _configthreadlocale(_ENABLE_PER_THREAD_LOCALE);
#endif
#ifdef HAS_LOCALE
  if ((locale_t)0 == caml_locale)
  {
#if defined(_MSC_VER)
    caml_locale = _create_locale(LC_NUMERIC, "C");
#else
    caml_locale = newlocale(LC_NUMERIC_MASK,"C",(locale_t)0);
#endif
  }
#endif
}

void caml_free_locale(void)
{
#ifdef HAS_LOCALE
  if ((locale_t)0 != caml_locale) freelocale(caml_locale);
  caml_locale = (locale_t)0;
#endif
}

CAMLexport value caml_copy_double(double d)
{
  value res;

#define Setup_for_gc
#define Restore_after_gc
  Alloc_small(res, Double_wosize, Double_tag);
#undef Setup_for_gc
#undef Restore_after_gc
  Store_double_val(res, d);
  return res;
}

#ifndef FLAT_FLOAT_ARRAY
CAMLexport void caml_Store_double_array_field(value val, mlsize_t i, double dbl)
{
  CAMLparam1 (val);
  value d = caml_copy_double (dbl);

  CAMLassert (Tag_val (val) != Double_array_tag);
  caml_modify (&Field(val, i), d);
  CAMLreturn0;
}
#endif /* ! FLAT_FLOAT_ARRAY */

CAMLprim value caml_format_float(value fmt, value arg)
{
  value res;
  double d = Double_val(arg);

#ifdef HAS_BROKEN_PRINTF
  if (isfinite(d)) {
#endif
    USE_LOCALE;
    res = caml_alloc_sprintf(String_val(fmt), d);
    RESTORE_LOCALE;
#ifdef HAS_BROKEN_PRINTF
  } else {
    if (isnan(d)) {
      res = caml_copy_string("nan");
    } else {
      if (d > 0)
        res = caml_copy_string("inf");
      else
        res = caml_copy_string("-inf");
    }
  }
#endif
  return res;
}

CAMLprim value caml_hexstring_of_float(value arg, value vprec, value vstyle)
{
  union { uint64_t i; double d; } u;
  int sign, exp;
  uint64_t m;
  char buffer[64];
  char * buf, * p;
  intnat prec;
  int d;
  value res;

  /* Allocate output buffer */
  prec = Long_val(vprec);
                  /* 12 chars for sign, 0x, decimal point, exponent */
  buf = (prec + 12 <= 64 ? buffer : caml_stat_alloc(prec + 12));
  /* Extract sign, mantissa, and exponent */
  u.d = Double_val(arg);
  sign = u.i >> 63;
  exp = (u.i >> 52) & 0x7FF;
  m = u.i & (((uint64_t) 1 << 52) - 1);
  /* Put sign */
  p = buf;
  if (sign) {
    *p++ = '-';
  } else {
    switch (Int_val(vstyle)) {
    case '+': *p++ = '+'; break;
    case ' ': *p++ = ' '; break;
    }
  }
  /* Treat special cases */
  if (exp == 0x7FF) {
    char * txt;
    if (m == 0) txt = "infinity"; else txt = "nan";
    memcpy(p, txt, strlen(txt));
    p[strlen(txt)] = 0;
    res = caml_copy_string(buf);
  } else {
    /* Output "0x" prefix */
    *p++ = '0'; *p++ = 'x';
    /* Normalize exponent and mantissa */
    if (exp == 0) {
      if (m != 0) exp = -1022;    /* denormal */
    } else {
      exp = exp - 1023;
      m = m | ((uint64_t) 1 << 52);
    }
    /* If a precision is given, and is small, round mantissa accordingly */
    prec = Long_val(vprec);
    if (prec >= 0 && prec < 13) {
      int i = 52 - prec * 4;
      uint64_t unit = (uint64_t) 1 << i;
      uint64_t half = unit >> 1;
      uint64_t mask = unit - 1;
      uint64_t frac = m & mask;
      m = m & ~mask;
      /* Round to nearest, ties to even */
      if (frac > half || (frac == half && (m & unit) != 0)) {
        m += unit;
      }
    }
    /* Leading digit */
    d = m >> 52;
    *p++ = (d < 10 ? d + '0' : d - 10 + 'a');
    m = (m << 4) & (((uint64_t) 1 << 56) - 1);
    /* Fractional digits.  If a precision is given, print that number of
       digits.  Otherwise, print as many digits as needed to represent
       the mantissa exactly. */
    if (prec >= 0 ? prec > 0 : m != 0) {
      *p++ = '.';
      while (prec >= 0 ? prec > 0 : m != 0) {
        d = m >> 52;
        *p++ = (d < 10 ? d + '0' : d - 10 + 'a');
        m = (m << 4) & (((uint64_t) 1 << 56) - 1);
        prec--;
      }
    }
    *p = 0;
    /* Add exponent */
    res = caml_alloc_sprintf("%sp%+d", buf, exp);
  }
  if (buf != buffer) caml_stat_free(buf);
  return res;
}

static int caml_float_of_hex(const char * s, const char * end, double * res)
{
  int64_t m = 0;                /* the mantissa - top 60 bits at most */
  int n_bits = 0;               /* total number of bits read */
  int m_bits = 0;               /* number of bits in mantissa */
  int x_bits = 0;               /* number of bits after mantissa */
  int dec_point = -1;           /* bit count corresponding to decimal point */
                                /* -1 if no decimal point seen */
  int exp = 0;                  /* exponent */
  char * p;                     /* for converting the exponent */
  double f;

  while (s < end) {
    char c = *s++;
    switch (c) {
    case '.':
      if (dec_point >= 0) return -1; /* multiple decimal points */
      dec_point = n_bits;
      break;
    case 'p': case 'P': {
      long e;
      if (*s == 0) return -1;   /* nothing after exponent mark */
      e = strtol(s, &p, 10);
      if (p != end) return -1;  /* ill-formed exponent */
      /* Handle exponents larger than int by returning 0/infinity directly.
         Mind that INT_MIN/INT_MAX are included in the test so as to capture
         the overflow case of strtol on Win64 -- long and int have the same
         size there. */
      if (e <= INT_MIN) {
        *res = 0.;
        return 0;
      }
      else if (e >= INT_MAX) {
        *res = m == 0 ? 0. : HUGE_VAL;
        return 0;
      }
      /* regular exponent value */
      exp = e;
      s = p;                    /* stop at next loop iteration */
      break;
    }
    default: {                  /* Nonzero digit */
      int d;
      if (c >= '0' && c <= '9') d = c - '0';
      else if (c >= 'A' && c <= 'F') d = c - 'A' + 10;
      else if (c >= 'a' && c <= 'f') d = c - 'a' + 10;
      else return -1;           /* bad digit */
      n_bits += 4;
      if (d == 0 && m == 0) break; /* leading zeros are skipped */
      if (m_bits < 60) {
        /* There is still room in m.  Add this digit to the mantissa. */
        m = (m << 4) + d;
        m_bits += 4;
      } else {
        /* We've already collected 60 significant bits in m.
           Now all we care about is whether there is a nonzero bit
           after. In this case, round m to odd so that the later
           rounding of m to FP produces the correct result. */
        if (d != 0) m |= 1;        /* round to odd */
        x_bits += 4;
      }
      break;
    }
    }
  }
  if (n_bits == 0) return -1;
  /* Convert mantissa to FP.  We use a signed conversion because we can
     (m has 60 bits at most) and because it is faster
     on several architectures. */
  f = (double) (int64_t) m;
  /* Adjust exponent to take decimal point and extra digits into account */
  {
    int adj = x_bits;
    if (dec_point >= 0) adj = adj + (dec_point - n_bits);
    /* saturated addition exp + adj */
    if (adj > 0 && exp > INT_MAX - adj)
      exp = INT_MAX;
    else if (adj < 0 && exp < INT_MIN - adj)
      exp = INT_MIN;
    else
      exp = exp + adj;
  }
  /* Apply exponent if needed */
  if (exp != 0) f = ldexp(f, exp);
  /* Done! */
  *res = f;
  return 0;
}

CAMLprim value caml_float_of_string(value vs)
{
  char parse_buffer[64];
  char * buf, * dst, * end;
  const char *src;
  mlsize_t len;
  int sign;
  double d;

  /* Remove '_' characters before conversion */
  len = caml_string_length(vs);
  buf = len < sizeof(parse_buffer) ? parse_buffer : caml_stat_alloc(len + 1);
  src = String_val(vs);
  dst = buf;
  while (len--) {
    char c = *src++;
    if (c != '_') *dst++ = c;
  }
  *dst = 0;
  if (dst == buf) goto error;
  /* Check for hexadecimal FP constant */
  src = buf;
  sign = 1;
  if (*src == '-') { sign = -1; src++; }
  else if (*src == '+') { src++; };
  if (src[0] == '0' && (src[1] == 'x' || src[1] == 'X')) {
    /* Convert using our hexadecimal FP parser */
    if (caml_float_of_hex(src + 2, dst, &d) == -1) goto error;
    if (sign < 0) d = -d;
  } else {
    /* Convert using strtod */
#if defined(HAS_STRTOD_L) && defined(HAS_LOCALE)
    d = strtod_l((const char *) buf, &end, caml_locale);
#else
    USE_LOCALE;
    d = strtod((const char *) buf, &end);
    RESTORE_LOCALE;
#endif /* HAS_STRTOD_L */
    if (end != dst) goto error;
  }
  if (buf != parse_buffer) caml_stat_free(buf);
  return caml_copy_double(d);
 error:
  if (buf != parse_buffer) caml_stat_free(buf);
  caml_failwith("float_of_string");
  return Val_unit; /* not reached */
}

CAMLprim value caml_int_of_float(value f)
{
  return Val_long((intnat) Double_val(f));
}

CAMLprim value caml_float_of_int(value n)
{
  return caml_copy_double((double) Long_val(n));
}

CAMLprim value caml_neg_float(value f)
{
  return caml_copy_double(- Double_val(f));
}

CAMLprim value caml_abs_float(value f)
{
  return caml_copy_double(fabs(Double_val(f)));
}

CAMLprim value caml_add_float(value f, value g)
{
  return caml_copy_double(Double_val(f) + Double_val(g));
}

CAMLprim value caml_sub_float(value f, value g)
{
  return caml_copy_double(Double_val(f) - Double_val(g));
}

CAMLprim value caml_mul_float(value f, value g)
{
  return caml_copy_double(Double_val(f) * Double_val(g));
}

CAMLprim value caml_div_float(value f, value g)
{
  return caml_copy_double(Double_val(f) / Double_val(g));
}

CAMLprim value caml_exp_float(value f)
{
  return caml_copy_double(exp(Double_val(f)));
}

CAMLexport double caml_trunc(double x)
{
#ifdef HAS_C99_FLOAT_OPS
  return trunc(x);
#else
  return (x >= 0.0)? floor(x) : ceil(x);
#endif
}

CAMLprim value caml_trunc_float(value f)
{
  return caml_copy_double(caml_trunc(Double_val(f)));
}

CAMLexport double caml_round(double f)
{
#ifdef HAS_C99_FLOAT_OPS
  return round(f);
#else
  union { uint64_t i; double d; } u, pred_one_half; /* predecessor of 0.5 */
  int e;  /* exponent */
  u.d = f;
  e = (u.i >> 52) & 0x7ff; /* - 0x3ff for the actual exponent */
  pred_one_half.i = 0x3FDFFFFFFFFFFFFF; /* 0x1.FFFFFFFFFFFFFp-2 */

  if (isfinite(f) && f != 0.) {
    if (e >= 52 + 0x3ff) return f; /* f is an integer already */
    if (f > 0.0)
      /* If we added 0.5 instead of its predecessor, then the
         predecessor of 0.5 would be rounded to 1. instead of 0. */
      return floor(f + pred_one_half.d);
    else
      return ceil(f - pred_one_half.d);
  }
  else
    return f;
#endif
}

CAMLprim value caml_round_float(value f)
{
  return caml_copy_double(caml_round(Double_val(f)));
}

CAMLprim value caml_floor_float(value f)
{
  return caml_copy_double(floor(Double_val(f)));
}

CAMLexport double caml_nextafter(double x, double y)
{
  return nextafter(x, y);
}

CAMLprim value caml_nextafter_float(value x, value y)
{
  return caml_copy_double(caml_nextafter(Double_val(x), Double_val(y)));
}

#ifndef HAS_WORKING_FMA
union double_as_int64 { double d; uint64_t i; };
#define IEEE754_DOUBLE_BIAS 0x3ff
#define IEEE_EXPONENT(N) (((N) >> 52) & 0x7ff)
#define IEEE_NEGATIVE(N) ((N) >> 63)
//C99 hexa float literals cannot be used, use pow() instead.
#define FL53    (pow(2,53)) //0x1p53
#define FLM53   (pow(2,-53)) //0x1p-53
#define FL54    (pow(2,54)) //0x1p54
#define FLM54   (pow(2,-54)) //0x1p-54
#define FL108   (pow(2,108)) //0x1p108
#define FLM108  (pow(2,-108)) //0x1p-108
#define FLM1074 (pow(2,-1074)) //0x1p-1074
#endif

CAMLexport double caml_fma(double x, double y, double z)
{
#ifdef HAS_WORKING_FMA
  return fma(x, y, z);
#else // Emulation of FMA, from S. Boldo and G. Melquiond, "Emulation
      // of a FMA and Correctly Rounded Sums: Proved Algorithms Using
      // Rounding to Odd," in IEEE Transactions on Computers, vol. 57,
      // no. 4, pp. 462-471, April 2008. Special cases implementation
      // comes from glibc's IEEE754 FMA emulation.
      // Only valid for double precision and round-to-nearest mode.

  union double_as_int64 u, v, w;
  union double_as_int64 ora;
  double mh, ml, xh, xl, yh, yl, t;
  double ah, al;
  double orah, oral;
  double t1, t2;
  double tiny;
  int neg, adjust = 0;
  u.d = x;
  v.d = y;
  w.d = z;

  if ( IEEE_EXPONENT(u.i) + IEEE_EXPONENT(v.i) >= 0x7FF +
       IEEE754_DOUBLE_BIAS - DBL_MANT_DIG
       || IEEE_EXPONENT(u.i) >= 0x7ff - DBL_MANT_DIG
       || IEEE_EXPONENT(v.i) >= 0x7ff - DBL_MANT_DIG
       || IEEE_EXPONENT(w.i) >= 0x7ff - DBL_MANT_DIG
       || IEEE_EXPONENT(u.i) + IEEE_EXPONENT(v.i) <=
         IEEE754_DOUBLE_BIAS + DBL_MANT_DIG )
    {
      /* If z is Inf, but x and y are finite, the result should be z
       * rather than NaN. */
      if (IEEE_EXPONENT(w.i) == 0x7ff &&
          IEEE_EXPONENT(u.i) != 0x7ff &&
          IEEE_EXPONENT(v.i) != 0x7ff)
              return (z + x) + y;
      /* If z is zero and x and y are nonzero, compute the result as
         x * y to avoid the wrong sign of a zero result if x * y
         underflows to 0. */
      if (z == 0 && x != 0 && y != 0)
        return x * y;
      /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
         x * y + z.  */
      if (IEEE_EXPONENT(u.i) == 0x7ff
          || IEEE_EXPONENT(v.i) == 0x7ff
          || IEEE_EXPONENT(w.i) == 0x7ff
          || x == 0
          || y == 0)
        return x * y + z;
      /* If fma will certainly overflow, compute as x * y. */
      if ((IEEE_EXPONENT(u.i) + IEEE_EXPONENT(v.i))
          > 0x7ff + IEEE754_DOUBLE_BIAS)
        return x * y;
      /* If x * y is less than 1/4 of DBL_TRUE_MIN, neither the result
         nor whether there is underflow depends on its exact value,
         only on its sign. */
      if (IEEE_EXPONENT(u.i) + IEEE_EXPONENT(v.i)
          < IEEE754_DOUBLE_BIAS - DBL_MANT_DIG - 2)
        {
          neg = IEEE_NEGATIVE(u.i) ^ IEEE_NEGATIVE(v.i) ;
          tiny = neg ? -FLM1074 : FLM1074;
          if (IEEE_EXPONENT(w.i) >= 3)
            return tiny + z;
          /* Scaling up, adding TINY and scaling down produces the
             correct result, because in round-to-nearest mode adding
             TINY has no effect and in other modes double rounding is
             harmless. But it may not produce required underflow
             exceptions. */
          v.d = z * FL54 + tiny;
          return v.d * FLM54;
        }
      if (IEEE_EXPONENT(u.i) + IEEE_EXPONENT(v.i)
          >= 0x7ff + IEEE754_DOUBLE_BIAS - DBL_MANT_DIG)
        {
          /* Compute 1p-53 times smaller result and multiply at the
             end.  */
          if (IEEE_EXPONENT(u.i) > IEEE_EXPONENT(v.i))
            x *= FLM53;
          else
            y *= FLM53;
          /* If x + y exponent is very large and z exponent is very small,
             it doesn't matter if we don't adjust it.  */
          if (IEEE_EXPONENT(w.i) > DBL_MANT_DIG)
            z *= FLM53;
          adjust = 1;
        }
      else if (IEEE_EXPONENT(w.i) >= 0x7ff - DBL_MANT_DIG)
        {
          /* Similarly. If z exponent is very large and x and y
             exponents are very small, adjust them up to avoid
             spurious underflows, rather than down.  */
          if (IEEE_EXPONENT(u.i) + IEEE_EXPONENT(v.i)
              <= IEEE754_DOUBLE_BIAS + 2 * DBL_MANT_DIG)
            {
              if (IEEE_EXPONENT(u.i) > IEEE_EXPONENT(v.i))
                x *= FL108;
              else
                y *= FL108;
            }
          else if (IEEE_EXPONENT(u.i) > IEEE_EXPONENT(v.i))
            {
              if (IEEE_EXPONENT(u.i) > DBL_MANT_DIG)
                x *= FLM53;
            }
          else if (IEEE_EXPONENT(v.i) > DBL_MANT_DIG)
            y *= FLM53;
          z *= FLM53;
          adjust = 1;
        }
      else if (IEEE_EXPONENT(u.i) >= 0x7ff - DBL_MANT_DIG)
        {
          x *= FLM53;
          y *= FL53;
        }
      else if (IEEE_EXPONENT(v.i) >= 0x7ff - DBL_MANT_DIG)
        {
          y *= FLM53;
          x *= FL53;
        }
      else /* if (IEEE_EXPONENT(u.i) + IEEE_EXPONENT(v.i) <=
              IEEE754_DOUBLE_BIAS + DBL_MANT_DIG) */
        {
          if (IEEE_EXPONENT(u.i) > IEEE_EXPONENT(v.i))
            x *= FL108;
          else
            y *= FL108;
          if (IEEE_EXPONENT(w.i) <= 4 * DBL_MANT_DIG + 6)
            {
              z *= FL108;
              adjust = -1;
            }
        }
    }

  /* Ensure correct sign of exact 0 + 0.  */
  if ((x == 0 || y == 0) && z == 0)
    return x * y + z;

  // Error-free multiplication: mh + ml = x * y
  mh = x * y;
  t = x * 134217729.0;
  xh = t - (t - x);
  xl = x - xh;
  t = y * 134217729.0;
  yh = t - (t - y);
  yl = y - yh;
  ml = xl * yl - (((mh - xh * yh) - xl * yh) - xh * yl);
  // Error-free addition: ah + al = z + mh
  ah = z + mh;
  t = ah - z;
  al = (z - (ah - t)) + (mh - t);

  /* If the result is an exact zero, ensure it has the correct sign. */
  if (ah == 0 && ml == 0)
    return z + mh;

  // Normalize ah, al, ml.
  t1 = al + ml;
  t = t1 - al;
  t2 = (al - (t1 - t)) + (ml - t);
  al = t1;
  ml = t2;
  t1 = ah + al;
  t = t1 - ah;
  t2 = (ah - (t1 - t)) + (al - t);
  ah = t1;
  al = t2;

  // Odd-rounded addition: ora = al + ml.
  orah = al + ml;
  oral = (al - orah) + ml;

  if ( oral != 0.0 )
    {
      ora.d = orah;
      if ( !(ora.i & 1) )
        {
          if ( (oral > 0.0) ^ (orah < 0.0) )
            ora.i++;
          else
            ora.i--;
          orah = ora.d;
        }
    }

  // Rounded addition: ra = ah + orah.
  if ( adjust > 0 )
    return (ah + orah) * FL53;
  else if ( adjust < 0 )
    return (ah + orah) * FLM108;
  else
    return ah + orah;
#endif
}

CAMLprim value caml_fma_float(value f1, value f2, value f3)
{
  return caml_copy_double(caml_fma(Double_val(f1),
                                   Double_val(f2), Double_val(f3)));
}

CAMLprim value caml_fmod_float(value f1, value f2)
{
  return caml_copy_double(fmod(Double_val(f1), Double_val(f2)));
}

CAMLprim value caml_frexp_float(value f)
{
  CAMLparam1 (f);
  CAMLlocal2 (res, mantissa);
  int exponent;

  mantissa = caml_copy_double(frexp (Double_val(f), &exponent));
  res = caml_alloc_tuple(2);
  Field(res, 0) = mantissa;
  Field(res, 1) = Val_int(exponent);
  CAMLreturn (res);
}

// Seems dumb but intnat could not correspond to int type.
double caml_ldexp_float_unboxed(double f, intnat i)
{
  return ldexp(f, (int) i);
}


CAMLprim value caml_ldexp_float(value f, value i)
{
  return caml_copy_double(ldexp(Double_val(f), Int_val(i)));
}

CAMLprim value caml_log_float(value f)
{
  return caml_copy_double(log(Double_val(f)));
}

CAMLprim value caml_log10_float(value f)
{
  return caml_copy_double(log10(Double_val(f)));
}

CAMLprim value caml_modf_float(value f)
{
  double frem;

  CAMLparam1 (f);
  CAMLlocal3 (res, quo, rem);

  quo = caml_copy_double(modf (Double_val(f), &frem));
  rem = caml_copy_double(frem);
  res = caml_alloc_tuple(2);
  Field(res, 0) = quo;
  Field(res, 1) = rem;
  CAMLreturn (res);
}

CAMLprim value caml_sqrt_float(value f)
{
  return caml_copy_double(sqrt(Double_val(f)));
}

CAMLprim value caml_power_float(value f, value g)
{
  return caml_copy_double(pow(Double_val(f), Double_val(g)));
}

CAMLprim value caml_sin_float(value f)
{
  return caml_copy_double(sin(Double_val(f)));
}

CAMLprim value caml_sinh_float(value f)
{
  return caml_copy_double(sinh(Double_val(f)));
}

CAMLprim value caml_cos_float(value f)
{
  return caml_copy_double(cos(Double_val(f)));
}

CAMLprim value caml_cosh_float(value f)
{
  return caml_copy_double(cosh(Double_val(f)));
}

CAMLprim value caml_tan_float(value f)
{
  return caml_copy_double(tan(Double_val(f)));
}

CAMLprim value caml_tanh_float(value f)
{
  return caml_copy_double(tanh(Double_val(f)));
}

CAMLprim value caml_asin_float(value f)
{
  return caml_copy_double(asin(Double_val(f)));
}

CAMLprim value caml_acos_float(value f)
{
  return caml_copy_double(acos(Double_val(f)));
}

CAMLprim value caml_atan_float(value f)
{
  return caml_copy_double(atan(Double_val(f)));
}

CAMLprim value caml_atan2_float(value f, value g)
{
  return caml_copy_double(atan2(Double_val(f), Double_val(g)));
}

CAMLprim value caml_ceil_float(value f)
{
  return caml_copy_double(ceil(Double_val(f)));
}

CAMLexport double caml_hypot(double x, double y)
{
#ifdef HAS_C99_FLOAT_OPS
  return hypot(x, y);
#else
  double tmp, ratio;
  x = fabs(x); y = fabs(y);
  if (x != x) /* x is NaN */
    return y > DBL_MAX ? y : x;  /* PR#6321 */
  if (y != y) /* y is NaN */
    return x > DBL_MAX ? x : y;  /* PR#6321 */
  if (x < y) { tmp = x; x = y; y = tmp; }
  if (x == 0.0) return 0.0;
  ratio = y / x;
  return x * sqrt(1.0 + ratio * ratio);
#endif
}

CAMLprim value caml_hypot_float(value f, value g)
{
  return caml_copy_double(caml_hypot(Double_val(f), Double_val(g)));
}

/* These emulations of expm1() and log1p() are due to William Kahan.
   See http://www.plunk.org/~hatch/rightway.php */
CAMLexport double caml_expm1(double x)
{
#ifdef HAS_C99_FLOAT_OPS
  return expm1(x);
#else
  double u = exp(x);
  if (u == 1.)
    return x;
  if (u - 1. == -1.)
    return -1.;
  return (u - 1.) * x / log(u);
#endif
}

CAMLexport double caml_log1p(double x)
{
#ifdef HAS_C99_FLOAT_OPS
  return log1p(x);
#else
  double u = 1. + x;
  if (u == 1.)
    return x;
  else
    return log(u) * x / (u - 1.);
#endif
}

CAMLprim value caml_expm1_float(value f)
{
  return caml_copy_double(caml_expm1(Double_val(f)));
}

CAMLprim value caml_log1p_float(value f)
{
  return caml_copy_double(caml_log1p(Double_val(f)));
}

union double_as_two_int32 {
    double d;
#if defined(ARCH_BIG_ENDIAN) || (defined(__arm__) && !defined(__ARM_EABI__))
    struct { uint32_t h; uint32_t l; } i;
#else
    struct { uint32_t l; uint32_t h; } i;
#endif
};

CAMLexport double caml_copysign(double x, double y)
{
#ifdef HAS_C99_FLOAT_OPS
  return copysign(x, y);
#else
  union double_as_two_int32 ux, uy;
  ux.d = x;
  uy.d = y;
  ux.i.h &= 0x7FFFFFFFU;
  ux.i.h |= (uy.i.h & 0x80000000U);
  return ux.d;
#endif
}

CAMLprim value caml_copysign_float(value f, value g)
{
  return caml_copy_double(caml_copysign(Double_val(f), Double_val(g)));
}

CAMLprim value caml_signbit(double x)
{
#ifdef HAS_C99_FLOAT_OPS
  return Val_bool(signbit(x));
#else
  union double_as_two_int32 ux;
  ux.d = x;
  return Val_bool(ux.i.h >> 31);
#endif
}

CAMLprim value caml_signbit_float(value f)
{
  return caml_signbit(Double_val(f));
}

CAMLprim value caml_neq_float(value f, value g)
{
  return Val_bool(Double_val(f) != Double_val(g));
}

#define DEFINE_NAN_CMP(op) (value f, value g) \
{ \
  return Val_bool(Double_val(f) op Double_val(g)); \
}

intnat caml_float_compare_unboxed(double f, double g)
{
  /* If one or both of f and g is NaN, order according to the convention
     NaN = NaN and NaN < x for all other floats x. */
  /* This branchless implementation is from GPR#164.
     Note that [f == f] if and only if f is not NaN.
     We expand each subresult of the expression to
     avoid sign-extension on 64bit. GPR#2250.
     See also translation of Pcompare_floats in asmcomp/cmmgen.ml  */
  intnat res =
    (intnat)(f > g) - (intnat)(f < g) + (intnat)(f == f) - (intnat)(g == g);
  return res;
}

CAMLprim value caml_eq_float DEFINE_NAN_CMP(==)
CAMLprim value caml_le_float DEFINE_NAN_CMP(<=)
CAMLprim value caml_lt_float DEFINE_NAN_CMP(<)
CAMLprim value caml_ge_float DEFINE_NAN_CMP(>=)
CAMLprim value caml_gt_float DEFINE_NAN_CMP(>)

CAMLprim value caml_float_compare(value vf, value vg)
{
  return Val_int(caml_float_compare_unboxed(Double_val(vf),Double_val(vg)));
}

enum { FP_normal, FP_subnormal, FP_zero, FP_infinite, FP_nan };

value caml_classify_float_unboxed(double vd)
{
#ifdef ARCH_SIXTYFOUR
  union { double d; uint64_t i; } u;
  uint64_t n;
  uint32_t e;

  u.d = vd;
  n = u.i << 1;                 /* shift sign bit off */
  if (n == 0) return Val_int(FP_zero);
  e = n >> 53;                  /* extract exponent */
  if (e == 0) return Val_int(FP_subnormal);
  if (e == 0x7FF) {
    if (n << 11 == 0)           /* shift exponent off */
      return Val_int(FP_infinite);
    else
      return Val_int(FP_nan);
  }
  return Val_int(FP_normal);
#else
  union double_as_two_int32 u;
  uint32_t h, l;

  u.d = vd;
  h = u.i.h;  l = u.i.l;
  l = l | (h & 0xFFFFF);
  h = h & 0x7FF00000;
  if ((h | l) == 0)
    return Val_int(FP_zero);
  if (h == 0)
    return Val_int(FP_subnormal);
  if (h == 0x7FF00000) {
    if (l == 0)
      return Val_int(FP_infinite);
    else
      return Val_int(FP_nan);
  }
  return Val_int(FP_normal);
#endif
}

CAMLprim value caml_classify_float(value vd)
{
  return caml_classify_float_unboxed(Double_val(vd));
}