/*
 * x : OUT : MPFR number extracted from the binary buffer, should have the same
 *           precision than the number in the binary format
 * buffer : IN : limb of the MPFR number x in a binary format,
 * buffer_size : IN/OUT : size of the buffer => size used in the buffer
 * return 0 if successful
 */
static int
mpfr_fpif_read_limbs (mpfr_t x, unsigned char *buffer, size_t *buffer_size)
{
  mpfr_prec_t precision;
  size_t nb_byte;
  size_t mp_bytes_per_limb;
  size_t nb_partial_byte;
  size_t i, j;

  precision = mpfr_get_prec (x);
  nb_byte = (precision + 7) >> 3;
  mp_bytes_per_limb = mp_bits_per_limb >> 3;
  nb_partial_byte = nb_byte % mp_bytes_per_limb;

  if ((buffer == NULL) || (*buffer_size < nb_byte))
    {
      *buffer_size = 0;
      return 1;
    }
  *buffer_size = nb_byte;

  if (nb_partial_byte > 0)
    {
      memset (MPFR_MANT(x), 0, sizeof(mp_limb_t));
      getBigEndianData ((unsigned char*) MPFR_MANT(x), buffer,
                        sizeof(mp_limb_t), nb_partial_byte);
    }
  for (i = nb_partial_byte, j = (nb_partial_byte == 0) ? 0 : 1; i < nb_byte;
       i += mp_bytes_per_limb, j++)
    getLittleEndianData ((unsigned char*) (MPFR_MANT(x) + j), buffer + i,
                         sizeof(mp_limb_t), sizeof(mp_limb_t));

  return 0;
}

void
mpfr_swap (mpfr_ptr u, mpfr_ptr v)
{
  mpfr_prec_t p1, p2;
  mpfr_sign_t s1, s2;
  mpfr_exp_t e1, e2;
  mp_limb_t *m1, *m2;

  p1 = MPFR_PREC(u);
  p2 = MPFR_PREC(v);
  MPFR_PREC(v) = p1;
  MPFR_PREC(u) = p2;

  s1 = MPFR_SIGN(u);
  s2 = MPFR_SIGN(v);
  MPFR_SIGN(v) = s1;
  MPFR_SIGN(u) = s2;

  e1 = MPFR_EXP(u);
  e2 = MPFR_EXP(v);
  MPFR_EXP(v) = e1;
  MPFR_EXP(u) = e2;

  m1 = MPFR_MANT(u);
  m2 = MPFR_MANT(v);
  MPFR_MANT(v) = m1;
  MPFR_MANT(u) = m2;
}

/*
 * buffer : OUT : store the limb of the MPFR number x in a binary format,
 *                can be null (may be reallocated if too small)
 * buffer_size : IN/OUT : size of the buffer => size used in the buffer
 * x : IN : MPFR number
 * return pointer to a buffer storing the limb of the MPFR number x in a binary
 *        format
 */
static unsigned char*
mpfr_fpif_store_limbs (unsigned char *buffer, size_t *buffer_size, mpfr_t x)
{
  unsigned char *result;
  mpfr_prec_t precision;
  size_t nb_byte;
  size_t nb_limb, mp_bytes_per_limb;
  size_t nb_partial_byte;
  size_t i, j;

  precision = mpfr_get_prec (x);
  nb_byte = (precision + 7) >> 3;
  mp_bytes_per_limb = mp_bits_per_limb >> 3;
  nb_partial_byte = nb_byte % mp_bytes_per_limb;
  nb_limb = (nb_byte + mp_bytes_per_limb - 1) / mp_bytes_per_limb;

  result = buffer;
  ALLOC_RESULT(result, buffer_size, nb_byte);

  putBigEndianData (result, (unsigned char*) MPFR_MANT(x),
                    sizeof(mp_limb_t), nb_partial_byte);
  for (i = nb_partial_byte, j = (nb_partial_byte == 0) ? 0 : 1; j < nb_limb;
       i += mp_bytes_per_limb, j++)
    putLittleEndianData (result + i, (unsigned char*) (MPFR_MANT(x) + j),
                         sizeof(mp_limb_t), sizeof(mp_limb_t));

  return result;
}

/*
 * Set f to z, choosing the smallest precision for f
 * so that z = f*(2^BPML)*zs*2^(RetVal)
 */
static int
set_z (mpfr_ptr f, mpz_srcptr z, mp_size_t *zs)
{
  mp_limb_t *p;
  mp_size_t s;
  int c;
  mp_prec_t pf;

  MPFR_ASSERTD (mpz_sgn (z) != 0);

  /* Remove useless ending 0 */
  for (p = PTR (z), s = *zs = ABS (SIZ (z)) ; *p == 0; p++, s--)
    MPFR_ASSERTD (s >= 0);

  /* Get working precision */
  count_leading_zeros (c, p[s-1]);
  pf = s * BITS_PER_MP_LIMB - c;
  if (pf < MPFR_PREC_MIN)
    pf = MPFR_PREC_MIN;
  mpfr_init2 (f, pf);

  /* Copy Mantissa */
  if (MPFR_LIKELY (c))
    mpn_lshift (MPFR_MANT (f), p, s, c);
  else
    MPN_COPY (MPFR_MANT (f), p, s);

  MPFR_SET_SIGN (f, mpz_sgn (z));
  MPFR_SET_EXP (f, 0);

  return -c;
}

static void *
start_routine (void *arg)
{
  mpfr_prec_t p;
  mpfr_t x;
  mpfr_prec_t inc = *(int *) arg;
  mp_limb_t *m;

  for (p = 100; p < 20000; p += 64 + 100 * (inc % 10))
    {
      mpfr_init2 (x, p);
      m = MPFR_MANT (x);
      mpfr_const_pi (x, MPFR_RNDD);
      mpfr_prec_round (x, 53, MPFR_RNDD);
      if (mpfr_cmp_str1 (x, "3.141592653589793116"))
        {
          printf ("mpfr_const_pi failed with threading\n");
          mpfr_out_str (stdout, 10, 0, x, MPFR_RNDN); putchar('\n');
          exit (1);
        }
      /* Check that no reallocation has been performed */
      MPFR_ASSERTN (m == MPFR_MANT (x));
      mpfr_clear (x);
    }

  pthread_exit (NULL);
}

int
mpfr_cmpabs (mpfr_srcptr b, mpfr_srcptr c)
{
  mp_exp_t be, ce;
  mp_size_t bn, cn;
  mp_limb_t *bp, *cp;

  if (MPFR_ARE_SINGULAR (b, c))
    {
      if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c))
        {
          MPFR_SET_ERANGE ();
          return 0;
        }
      else if (MPFR_IS_INF (b))
        return ! MPFR_IS_INF (c);
      else if (MPFR_IS_INF (c))
        return -1;
      else if (MPFR_IS_ZERO (c))
        return ! MPFR_IS_ZERO (b);
      else /* b == 0 */
        return -1;
    }

  be = MPFR_GET_EXP (b);
  ce = MPFR_GET_EXP (c);
  if (be > ce)
    return 1;
  if (be < ce)
    return -1;

  /* exponents are equal */

  bn = MPFR_LIMB_SIZE(b)-1;
  cn = MPFR_LIMB_SIZE(c)-1;

  bp = MPFR_MANT(b);
  cp = MPFR_MANT(c);

  for ( ; bn >= 0 && cn >= 0; bn--, cn--)
    {
      if (bp[bn] > cp[cn])
        return 1;
      if (bp[bn] < cp[cn])
        return -1;
    }

  for ( ; bn >= 0; bn--)
    if (bp[bn])
      return 1;

  for ( ; cn >= 0; cn--)
    if (cp[cn])
      return -1;

   return 0;
}

/* set bit n of x to b, where bit 0 is the most significant one */
static void
set_bit (mpfr_t x, unsigned int n, int b)
{
  unsigned l;
  mp_size_t xn;

  xn = (MPFR_PREC(x) - 1) / mp_bits_per_limb;
  l = n / mp_bits_per_limb;
  n %= mp_bits_per_limb;
  n = mp_bits_per_limb - 1 - n;
  if (b)
    MPFR_MANT(x)[xn - l] |= (mp_limb_t) 1 << n;
  else
    MPFR_MANT(x)[xn - l] &= ~((mp_limb_t) 1 << n);
}

void
mpfr_random2 (mpfr_ptr x, mp_size_t size, mp_exp_t exp)
{
  mp_size_t xn;
  unsigned long cnt;
  mp_ptr xp = MPFR_MANT(x), yp[1];
  mp_size_t prec = (MPFR_PREC(x) - 1)/BITS_PER_MP_LIMB; 

  MPFR_CLEAR_FLAGS(x);
  xn = ABS (size);
  if (xn != 0)
    {
      if (xn > prec + 1)
	xn = prec + 1;

      mpn_random2 (xp, xn);
    }

  if (exp != 0) {
    /* use mpn_random instead of random since that function is not
       available on all platforms (for example HPUX, DEC OSF, ...) */
    mpn_random ((mp_limb_t*) yp, 1);
    exp = (mp_exp_t) yp[0] % (2 * exp) - exp;
  }

  count_leading_zeros(cnt, xp[xn - 1]); 
  if (cnt) mpn_lshift(xp, xp, xn, cnt); 
  MPFR_EXP(x) = exp-cnt; 
  cnt = xn*BITS_PER_MP_LIMB - prec; 
  /* cnt is the number of non significant bits in the low limb */
  xp[0] &= ~((MP_LIMB_T_ONE << cnt) - MP_LIMB_T_ONE);
}

MPFR_HOT_FUNCTION_ATTR void
mpfr_clear (mpfr_ptr m)
{
  (*__gmp_free_func) (MPFR_GET_REAL_PTR (m),
                      MPFR_MALLOC_SIZE (MPFR_GET_ALLOC_SIZE (m)));
  MPFR_MANT (m) = (mp_limb_t *) 0;
}

/* sets x to x+sign(x)*ulp(x) */
int
mpfr_add_one_ulp (mpfr_ptr x, mp_rnd_t rnd_mode)
{
  mp_size_t xn;
  int sh;
  mp_limb_t *xp;

  if (MPFR_IS_NAN(x))
    MPFR_RET_NAN;

  if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x))
    return 0;

  xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
  sh = xn * BITS_PER_MP_LIMB - MPFR_PREC(x);
  xp = MPFR_MANT(x);
  if (mpn_add_1 (xp, xp, xn, MP_LIMB_T_ONE << sh)) /* got 1.0000... */
    {
      mp_exp_t exp = MPFR_EXP(x);
      if (exp == __mpfr_emax)
        return mpfr_set_overflow(x, rnd_mode, MPFR_SIGN(x));
      else
        {
          MPFR_EXP(x)++;
          xp[xn-1] = GMP_LIMB_HIGHBIT;
        }
    }
  return 0;
}

int
mpfr_set_overflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign)
{
  int inex;

  MPFR_CLEAR_FLAGS(x);
  if ((rnd_mode == GMP_RNDU && sign < 0)
   || (rnd_mode == GMP_RNDD && sign > 0))
    {
      mp_size_t xn, i;
      int sh;
      mp_limb_t *xp;

      MPFR_EXP(x) = __mpfr_emax;
      xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB;
      sh = xn * BITS_PER_MP_LIMB - MPFR_PREC(x);
      xp = MPFR_MANT(x);
      xp[0] = MP_LIMB_T_MAX << sh;
      for (i = 1; i < xn; i++)
        xp[i] = MP_LIMB_T_MAX;
      inex = -1;
    }
  else
    {
      MPFR_SET_INF(x);
      inex = 1;
    }
  if (MPFR_SIGN(x) != sign)
    MPFR_CHANGE_SIGN(x);
  __mpfr_flags |= MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
  return sign > 0 ? inex : -inex;
}

int
mpfr_set_underflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign)
{
  int inex;

  MPFR_CLEAR_FLAGS(x);
  if ((rnd_mode == GMP_RNDU && sign > 0)
   || (rnd_mode == GMP_RNDD && sign < 0))
    {
      mp_size_t xn;
      mp_limb_t *xp;

      MPFR_EXP(x) = __mpfr_emin;
      xn = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB;
      xp = MPFR_MANT(x);
      xp[xn] = GMP_LIMB_HIGHBIT;
      MPN_ZERO(xp, xn);
      inex = 1;
    }
  else
    {
      MPFR_SET_ZERO(x);
      inex = -1;
    }
  if (MPFR_SIGN(x) != sign)
    MPFR_CHANGE_SIGN(x);
  __mpfr_flags |= MPFR_FLAGS_INEXACT | MPFR_FLAGS_UNDERFLOW;
  return sign > 0 ? inex : -inex;
}

int
mpfr_integer_p (mpfr_srcptr x)
{
  mpfr_exp_t expo;
  mpfr_prec_t prec;
  mp_size_t xn;
  mp_limb_t *xp;

  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x)))
    return (MPFR_IS_ZERO(x));

  expo = MPFR_GET_EXP (x);
  if (expo <= 0)
    return 0;

  prec = MPFR_PREC(x);
  if ((mpfr_uexp_t) expo >= (mpfr_uexp_t) prec)
    return 1;

  /* 0 < expo < prec */

  xn =  (mp_size_t) ((prec - 1) / GMP_NUMB_BITS);  /* index of last limb */
  xn -= (mp_size_t) (expo / GMP_NUMB_BITS);
  /* now the index of the last limb containing bits of the fractional part */

  xp = MPFR_MANT(x);
  MPFR_ASSERTN(xn >= 0);
  if (xp[xn] << (expo % GMP_NUMB_BITS) != 0)
    return 0;
  while (--xn >= 0)
    if (xp[xn] != 0)
      return 0;
  return 1;
}

void
mpfr_custom_init_set (mpfr_ptr x, int kind, mp_exp_t exp,
                     mp_prec_t prec, void *mantissa)
{
  mpfr_kind_t t;
  int s;
  mp_exp_t e;

  if (kind >= 0)
    {
      t = (mpfr_kind_t) kind;
      s = MPFR_SIGN_POS;
    }
  else
    {
      t = (mpfr_kind_t) -kind;
      s = MPFR_SIGN_NEG;
    }
  MPFR_ASSERTD (t <= MPFR_REGULAR_KIND);
  e = MPFR_LIKELY (t == MPFR_REGULAR_KIND) ? exp :
    MPFR_UNLIKELY (t == MPFR_NAN_KIND) ? MPFR_EXP_NAN :
    MPFR_UNLIKELY (t == MPFR_INF_KIND) ? MPFR_EXP_INF : MPFR_EXP_ZERO;

  MPFR_PREC (x) = prec;
  MPFR_SET_SIGN (x, s);
  MPFR_EXP (x) = e;
  MPFR_MANT (x) = (mp_limb_t*) mantissa;
  return;
}

// remove trailing zeros (only by limbs)
void mpfr_remove_trailing_zeros(mpfr_t x) {
  unsigned int xn = MPFR_LIMB_SIZE(x); 
  mp_limb_t* xp = MPFR_MANT(x);
  unsigned int i = 0;
  while (i < xn && xp[i] == 0) i++;
  if (i > 0 && i < xn) mpfr_round_prec(x, GMP_RNDN, (xn-i)*BITS_PER_MP_LIMB); 
}

int
mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
{
  mp_ptr rp;
  mp_prec_t nbits;
  mp_size_t nlimbs;
  mp_size_t k; /* number of high zero limbs */
  mp_exp_t exp;
  int cnt;

  MPFR_CLEAR_FLAGS (rop);

  rp = MPFR_MANT (rop);
  nbits = MPFR_PREC (rop);
  nlimbs = MPFR_LIMB_SIZE (rop);
  MPFR_SET_POS (rop);

  /* Uniform non-normalized significand */
  _gmp_rand (rp, rstate, nlimbs * BITS_PER_MP_LIMB);

  /* If nbits isn't a multiple of BITS_PER_MP_LIMB, mask the low bits */
  cnt = nlimbs * BITS_PER_MP_LIMB - nbits;
  if (MPFR_LIKELY (cnt != 0))
    rp[0] &= ~MPFR_LIMB_MASK (cnt);

  /* Count the null significant limbs and remaining limbs */
  exp = 0;
  k = 0;
  while (nlimbs != 0 && rp[nlimbs - 1] == 0)
    {
      k ++;
      nlimbs --;
      exp -= BITS_PER_MP_LIMB;
    }

  if (MPFR_LIKELY (nlimbs != 0)) /* otherwise value is zero */
    {
      count_leading_zeros (cnt, rp[nlimbs - 1]);
      /* Normalization */
      if (mpfr_set_exp (rop, exp - cnt))
        {
          /* If the exponent is not in the current exponent range, we
             choose to return a NaN as this is probably a user error.
             Indeed this can happen only if the exponent range has been
             reduced to a very small interval and/or the precision is
             huge (very unlikely). */
          MPFR_SET_NAN (rop);
          __gmpfr_flags |= MPFR_FLAGS_NAN; /* Can't use MPFR_RET_NAN */
          return 1;
        }
      if (cnt != 0)
        mpn_lshift (rp + k, rp, nlimbs, cnt);
      if (k != 0)
        MPN_ZERO (rp, k);
    }
  else
    MPFR_SET_ZERO (rop);

  return 0;
}

int
mpfr_urandomb (mpfr_ptr rop, gmp_randstate_t rstate)
{
    mpfr_limb_ptr rp;
    mpfr_prec_t nbits;
    mp_size_t nlimbs;
    mp_size_t k; /* number of high zero limbs */
    mpfr_exp_t exp;
    int cnt;

    rp = MPFR_MANT (rop);
    nbits = MPFR_PREC (rop);
    nlimbs = MPFR_LIMB_SIZE (rop);
    MPFR_SET_POS (rop);
    cnt = nlimbs * GMP_NUMB_BITS - nbits;

    /* Uniform non-normalized significand */
    /* generate exactly nbits so that the random generator stays in the same
       state, independent of the machine word size GMP_NUMB_BITS */
    mpfr_rand_raw (rp, rstate, nbits);
    if (MPFR_LIKELY (cnt != 0)) /* this will put the low bits to zero */
        mpn_lshift (rp, rp, nlimbs, cnt);

    /* Count the null significant limbs and remaining limbs */
    exp = 0;
    k = 0;
    while (nlimbs != 0 && rp[nlimbs - 1] == 0)
    {
        k ++;
        nlimbs --;
        exp -= GMP_NUMB_BITS;
    }

    if (MPFR_LIKELY (nlimbs != 0)) /* otherwise value is zero */
    {
        count_leading_zeros (cnt, rp[nlimbs - 1]);
        /* Normalization */
        if (mpfr_set_exp (rop, exp - cnt))
        {
            /* If the exponent is not in the current exponent range, we
               choose to return a NaN as this is probably a user error.
               Indeed this can happen only if the exponent range has been
               reduced to a very small interval and/or the precision is
               huge (very unlikely). */
            MPFR_SET_NAN (rop);
            __gmpfr_flags |= MPFR_FLAGS_NAN; /* Can't use MPFR_RET_NAN */
            return 1;
        }
        if (cnt != 0)
            mpn_lshift (rp + k, rp, nlimbs, cnt);
        if (k != 0)
            MPN_ZERO (rp, k);
    }
    else
        MPFR_SET_ZERO (rop);

    return 0;
}

mpfr_prec_t
mpfr_min_prec (mpfr_srcptr x)
{
  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
    return 0;

  /* from a suggestion by Andreas Enge (2010-11-18) */
  return MPFR_LIMB_SIZE (x) * GMP_NUMB_BITS - mpn_scan1 (MPFR_MANT (x), 0);
}

static void
print_mpfr (mpfr_srcptr x, const char *name)
{
  unsigned char temp[16]; /* buffer for the base-256 string */
  unsigned char *ptr;     /* pointer to its first non-zero byte */
  int size;               /* size of the string */
  int i;                  /* bits2use index */
  int j;                  /* output limb index */
  int k;                  /* byte index (in output limb) */
  int r;                  /* digit index, relative to ptr */
  char prefix[12];        /* "0x" or "UINT64_C(0x" */
  char suffix[2];         /* "" or ")" */

  if (printf ("#if 0\n") < 0)
    { fprintf (stderr, "Error in printf\n"); exit (1); }
  for (i = 0; i < size_of_bits2use; i++)
    {
      if (bits2use[i] == 64)
        {
          strcpy (prefix, "UINT64_C(0x");
          strcpy (suffix, ")");
        }
      else
        {
          strcpy (prefix, "0x");
          strcpy (suffix, "");
        }
      if (printf ("#elif GMP_NUMB_BITS == %d\n"
                  "const mp_limb_t %s__tab[] = { %s", bits2use[i], name,
                  prefix) < 0)
        { fprintf (stderr, "Error in printf\n"); exit (1); }
      size = mpn_get_str (temp, 256, MPFR_MANT (x), MPFR_LIMB_SIZE (x));
      MPFR_ASSERTN (size <= 16);
      ptr = temp;
      /* Skip leading zeros. */
      while (*ptr == 0)
        {
          ptr++;
          size--;
          MPFR_ASSERTN (size > 0);
        }
      MPFR_ASSERTN (*ptr >= 128);
      for (j = (MPFR_PREC (x) - 1) / bits2use[i]; j >= 0; j--)
        {
          r = j * (bits2use[i] / 8);
          for (k = 0; k < bits2use[i] / 8; k++)
            if (printf ("%02x", r < size ? ptr[r++] : 0) < 0)
              { fprintf (stderr, "Error in printf\n"); exit (1); }
          if (j == 0 && printf ("%s };\n", suffix) < 0)
            { fprintf (stderr, "Error in printf\n"); exit (1); }
          else if (j > 0 && printf ("%s, %s", suffix, prefix) < 0)
            { fprintf (stderr, "Error in printf\n"); exit (1); }
        }
    }
  if (printf ("#endif\n\n") < 0)
    { fprintf (stderr, "Error in printf\n"); exit (1); }
}

int
mpfr_set_si (mpfr_ptr x, long i, mp_rnd_t rnd_mode)
{
  int inex;
  mp_size_t xn;
  unsigned int cnt, nbits;
  mp_limb_t ai, *xp;

  MPFR_CLEAR_FLAGS(x);
  if (i == 0)
    {
      MPFR_SET_ZERO(x);
      MPFR_SET_POS(x);
      MPFR_RET(0);
    }

  xn = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB;
  ai = SAFE_ABS(long, i);
  count_leading_zeros(cnt, ai);

  xp = MPFR_MANT(x);
  xp[xn] = ai << cnt;
  /* don't forget to put zero in lower limbs */
  MPN_ZERO(xp, xn);
  /* set sign */
  if ((i < 0) ^ (MPFR_SIGN(x) < 0))
    MPFR_CHANGE_SIGN(x);

  MPFR_EXP(x) = nbits = BITS_PER_MP_LIMB - cnt;
  inex = mpfr_check_range(x, rnd_mode);
  if (inex)
    return inex; /* underflow or overflow */

  /* round if MPFR_PREC(x) smaller than length of i */
  if (MPFR_PREC(x) < nbits)
    {
      int carry;

      carry = mpfr_round_raw(xp+xn, xp+xn, nbits, (i < 0), MPFR_PREC(x),
                             rnd_mode, &inex);
      if (carry)
        {
          mp_exp_t exp = MPFR_EXP(x);

          if (exp == __mpfr_emax)
            return mpfr_set_overflow(x, rnd_mode, (i < 0 ? -1 : 1));

          MPFR_EXP(x)++;
          xp[xn] = GMP_LIMB_HIGHBIT;
        }
    }

  MPFR_RET(inex);
}