/* mpfr_set_ld -- convert a machine long double to
                  a multiple precision floating-point number

Copyright 2002, 2003, 2004, 2005 Free Software Foundation, Inc.

This file is part of the MPFR Library.

The MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.

The MPFR Library 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 General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the MPFR Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */

#include <float.h>

#define MPFR_NEED_LONGLONG_H
#include "mpfr-impl.h"

/* Various i386 systems have been seen with float.h LDBL constants equal to
   the DBL ones, whereas they ought to be bigger, reflecting the 10-byte
   IEEE extended format on that processor.  gcc 3.2.1 on FreeBSD and Solaris
   has been seen with the problem, and gcc 2.95.4 on FreeBSD 4.7.  */

#if HAVE_LDOUBLE_IEEE_EXT_LITTLE
static const struct {
  char         bytes[10];
  long double  dummy;  /* for memory alignment */
} ldbl_max_struct = {
  { '\377','\377','\377','\377',
    '\377','\377','\377','\377',
    '\376','\177' }, 0.0
};
#define MPFR_LDBL_MAX   (* (const long double *) ldbl_max_struct.bytes)
#else
#define MPFR_LDBL_MAX   LDBL_MAX
#endif

#ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE

/* Generic code */
int
mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode)
{
  mpfr_t t, u;
  int inexact, shift_exp = 0, inexact2;
  MPFR_SAVE_EXPO_DECL (expo);

  /* Check for NAN */
  LONGDOUBLE_NAN_ACTION (d, goto nan);
  
  /* Check for INF */
  if (d > MPFR_LDBL_MAX)
    {
      mpfr_set_inf (r, 1);
      return 0;
    }
  else if (d < -MPFR_LDBL_MAX)
    {
      mpfr_set_inf (r, -1);
      return 0;
    }
  /* Check for ZERO */
  else if (d == 0.0)
    return mpfr_set_d (r, (double) d, rnd_mode);

  mpfr_init2 (t, MPFR_LDBL_MANT_DIG);
  mpfr_init2 (u, IEEE_DBL_MANT_DIG);

  MPFR_SET_ZERO (t);
  MPFR_SAVE_EXPO_MARK (expo);

  while (d != (long double) 0.0)
    {
      /* Check overflow of Double */
      if (d > (long double) DBL_MAX || (-d) > (long double) DBL_MAX)
	{
          long double div9, div10, div11, div12, div13;

#define TWO_64 18446744073709551616.0 /* 2^64 */
#define TWO_128 (TWO_64 * TWO_64)
#define TWO_256 (TWO_128 * TWO_128)
          div9 = (long double) (double) (TWO_256 * TWO_256); /* 2^(2^9) */
          div10 = div9 * div9;
          div11 = div10 * div10; /* 2^(2^11) */
          div12 = div11 * div11; /* 2^(2^12) */
          div13 = div12 * div12; /* 2^(2^13) */
          if (ABS (d) >= div13)
            {
              d = d / div13; /* exact */
              shift_exp += 8192;
            }
          if (ABS (d) >= div12)
            {
              d = d / div12; /* exact */
              shift_exp += 4096;
            }
          if (ABS (d) >= div11)
            {
              d = d / div11; /* exact */
              shift_exp += 2048;
            }
          if (ABS (d) >= div10)
            {
              d = d / div10; /* exact */
              shift_exp += 1024;
            }
          /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < d < 2^1024,
             therefore we have one extra exponent reduction step */
          if (ABS (d) >= div9)
            {
              d = d / div9; /* exact */
              shift_exp += 512;
            }
        }
      /* Check underflow of double */
      else if (d < (long double)DBL_MIN && (-d) < (long double)DBL_MIN)
        {
	  long double div10, div11, div12, div13;
	  div10 = (long double) (double) 5.5626846462680034577255e-309; /* 2^(-2^10) */
	  div11 = div10 * div10; /* 2^(-2^11) */
	  div12 = div11 * div11; /* 2^(-2^12) */
	  div13 = div12 * div12; /* 2^(-2^13) */
	  if (ABS (d) <= div13)
	    {
	      d = d / div13; /* exact */
	      shift_exp -= 8192;
	    }
	  if (ABS (d) <= div12)
	    {
	      d = d / div12; /* exact */
	      shift_exp -= 4096;
	    }
	  if (ABS (d) <= div11)
	    {
	      d = d / div11; /* exact */
	      shift_exp -= 2048;
	    }
	  if (ABS (d) <= div10)
	    {
	      d = d / div10; /* exact */
	      shift_exp -= 1024;
	    }
	}
      else
	{
	  /* since DBL_MIN < ABS(d) < DBL_MAX, the cast to double should not
	     overflow here */
	  MPFR_ASSERTD ((long double) DBL_MIN <= ABS (d));
	  MPFR_ASSERTD ((long double) DBL_MAX >= ABS (d));
	  inexact = mpfr_set_d (u, (double) d, GMP_RNDN);
	  MPFR_ASSERTD (inexact == 0);
	  MPFR_ASSERTD (!MPFR_IS_ZERO (u));
	  inexact = mpfr_add (t, t, u, GMP_RNDN); /* exact */
	  MPFR_ASSERTD (inexact == 0);
	  d = d - (long double) mpfr_get_d1 (u);
	}
    }
  /* now t is exactly the input value d */
  inexact = mpfr_set (r, t, rnd_mode);
  inexact2 = mpfr_mul_2si (r, r, shift_exp, rnd_mode);
  if (inexact2) /* overflow or underflow */
    inexact = inexact2;
  mpfr_clear (t);
  mpfr_clear (u);

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (r, inexact, rnd_mode);

 nan:
  MPFR_SET_NAN(r);
  MPFR_RET_NAN;
}

#else /* IEEE Extended Little Endian Code */

int
mpfr_set_ld (mpfr_ptr r, long double d, mp_rnd_t rnd_mode)
{
  int inexact, i, k, cnt;
  mpfr_t tmp;
  mp_limb_t tmpmant[MPFR_LIMBS_PER_LONG_DOUBLE];
  mpfr_long_double_t x;
  mp_exp_t exp;
  int signd;
  MPFR_SAVE_EXPO_DECL (expo);

  /* Check for NAN */
  if (MPFR_UNLIKELY (d != d))
    { 
      MPFR_SET_NAN (r);
      MPFR_RET_NAN;
    }  
  /* Check for INF */
  else if (MPFR_UNLIKELY (d > MPFR_LDBL_MAX))
    {
      MPFR_SET_INF (r);
      MPFR_SET_POS (r);
      return 0;
    }
  else if (MPFR_UNLIKELY (d < -MPFR_LDBL_MAX))
    {
      MPFR_SET_INF (r);
      MPFR_SET_NEG (r);
      return 0;
    }
  /* Check for ZERO */
  else if (MPFR_UNLIKELY (d == 0.0))
    {
      x.ld = d;
      MPFR_SET_ZERO (r);
      if (x.s.sign == 1)
	MPFR_SET_NEG(r);
      else
	MPFR_SET_POS(r);
      return 0;
    }
  
  /* now d is neither 0, nor NaN nor Inf */
  MPFR_SAVE_EXPO_MARK (expo);
  
  MPFR_MANT (tmp) = tmpmant;
  MPFR_PREC (tmp) = 64;

  /* Extract sign */
  x.ld = d;
  signd = MPFR_SIGN_POS;
  if (x.ld < 0.0) {
    signd = MPFR_SIGN_NEG;
    x.ld = -x.ld;
  }

  /* Extract mantissa */
#if BITS_PER_MP_LIMB >= 64
  tmpmant[0] = ((mp_limb_t) x.s.manh << 32) | ((mp_limb_t) x.s.manl);
#else
  tmpmant[0] = (mp_limb_t) x.s.manl;
  tmpmant[1] = (mp_limb_t) x.s.manh;
#endif
  
  /* Normalize mantissa */
  i = MPFR_LIMBS_PER_LONG_DOUBLE;
  MPN_NORMALIZE_NOT_ZERO (tmpmant, i);
  k = MPFR_LIMBS_PER_LONG_DOUBLE - i;
  count_leading_zeros (cnt, tmpmant[i - 1]);
  if (MPFR_LIKELY (cnt != 0))
    mpn_lshift (tmpmant + k, tmpmant, i, cnt);
  else if (k != 0)
    MPN_COPY (tmpmant + k, tmpmant, i);  
  if (MPFR_UNLIKELY (k != 0))
    MPN_ZERO (tmpmant, k);

  /* Set exponent */
  if (x.s.exph == 0 && x.s.expl == 0)
    exp = -0x3FFD;
  else
    exp = (x.s.exph << 8) + x.s.expl - 0x3FFE;

  MPFR_SET_EXP (tmp, exp - cnt - k * BITS_PER_MP_LIMB);

  /* tmp is exact */
  inexact = mpfr_set4 (r, tmp, rnd_mode, signd);

  MPFR_SAVE_EXPO_FREE (expo);
  return mpfr_check_range (r, inexact, rnd_mode);
}

#endif