/* mpfr_get_str -- output a floating-point number to a string

Copyright (C) 1999, 2001 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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
#include "mpfr.h"
#include "mpfr-impl.h"

/*
  Convert op to a string in base 'base' with 'n' digits and writes the 
  mantissa in 'str', the exponent in 'expptr'.
  The result is rounded wrt 'rnd_mode'.

  For op = 3.1416 we get str = "31416" and expptr=1.
 */
char*
mpfr_get_str (char *str, mp_exp_t *expptr, int base, size_t n,
		  mpfr_srcptr op, mp_rnd_t rnd_mode)
{
  double d;
  long e, q, div, p, err, prec, sh;
  mpfr_t a, b;
  mpz_t bz;
  char *str0=NULL;
  mp_rnd_t rnd1;
  int f, pow2, ok=0, neg, str_is_null=(str==NULL);

  if (n == 0) /* determine n from precision of op */
    {
      n = (int) (__mp_bases[base].chars_per_bit_exactly * MPFR_PREC(op));
      if (n < MPFR_PREC_MIN)
        n = MPFR_PREC_MIN;
    }

  MPFR_ASSERTN(n >= MPFR_PREC_MIN);

  if (base < 2 || 36 < base)
    {
      fprintf (stderr,
               "Error: too small or too large base in mpfr_get_str: %d\n",
               base);
      exit (1);
    }
  
  neg = MPFR_SIGN(op) < 0;

  if (MPFR_IS_INF(op))
    { 
      if (str == NULL)
        str = (*__gmp_allocate_func)(neg + 4);
      str0 = str;
      if (neg)
        *str++ = '-';
      *str++ = 'I';
      *str++ = 'n';
      *str++ = 'f';
      *str = '\0';
      return str0; /* strlen(str0) = neg + 3 */
    }

  if (MPFR_IS_ZERO(op))
    {
      if (str == NULL)
        str = (*__gmp_allocate_func)(neg + n + 1);
      str0 = str;
      if (neg)
        *str++ = '-';
      for (f=0;f<n;f++)
        *str++ = '0';
      *str++ = '\0';
      *expptr = 1;
      return str0; /* strlen(str0) = neg + n */
    }

  count_leading_zeros (pow2, (mp_limb_t) base);
  pow2 = BITS_PER_MP_LIMB - pow2 - 1;
  if (base != (1<<pow2))
    pow2=0; 
  /* if pow2 <> 0, then base = 2^pow2 */

  /* first determines the exponent */
  e = MPFR_EXP(op);
  d = ABS(mpfr_get_d2(op, 0));
  /* the absolute value of op is between 1/2*2^e and 2^e */
  /* the output exponent f is such that base^(f-1) <= |op| < base^f
     i.e. f = 1 + floor(log(|op|)/log(base))
     = 1 + floor((log(|m|)+e*log(2))/log(base)) */
  /* f = 1 + (int) floor((log(d)/LOG2+(double)e)*LOG2/log((double)base)); */
  /* when base = 2^pow2, then |op| < 2^(pow2*f)
     i.e. e <= pow2*f and f = ceil(e/pow2) */
  if (pow2)
      f = ((e < 0) ? e : (e + pow2 - 1)) / pow2;
  else
    {
      d = ((double) e + (double) _mpfr_floor_log2 (d))
                      * __mp_bases[base].chars_per_bit_exactly;
      /* warning: (int) d rounds towards 0 */
      f = (int) d; /* f equals floor(d) if d >= 0 and ceil(d) if d < 0 */
      if (f <= d)
        f++;
    }

  mpfr_init (a);
  mpfr_init (b);
  mpz_init (bz);

  str0 = str;

 start_again:
  /* now the first n digits of the mantissa are obtained from
     rnd(op*base^(n-f)) */
  if (pow2)
    prec = n*pow2;
  else
    prec = 1 + (long) ((double) n / __mp_bases[base].chars_per_bit_exactly);
  err = 5;
  q = prec + err;
  /* one has to use at least q bits */
  q = (((q-1)/BITS_PER_MP_LIMB)+1) * BITS_PER_MP_LIMB;
  mpfr_set_prec (a, q);
  mpfr_set_prec (b, q);

  do {
    p = n - f;
    if ((div = (p<0)))
      p = -p;
    rnd1 = rnd_mode;
    if (div)
      {
        /* if div we divide by base^p so we have to invert the rounding mode */
        switch (rnd1)
          {
          case GMP_RNDN: rnd1=GMP_RNDN; break;
          case GMP_RNDZ: rnd1=GMP_RNDU; break;
          case GMP_RNDU: rnd1=GMP_RNDZ; break;
          case GMP_RNDD: rnd1=GMP_RNDZ; break;
          }
      }

    if (pow2)
      {
	if (div)
	  mpfr_div_2exp (b, op, pow2*p, rnd_mode);
	else
	  mpfr_mul_2exp (b, op, pow2*p, rnd_mode);
      }
    else
      {
        /* compute base^p with q bits and rounding towards zero */
        mpfr_set_prec (b, q);
        if (p == 0)
          {
            mpfr_set (b, op, rnd_mode);
            mpfr_set_ui (a, 1, rnd_mode);
          }
        else
          {
            mpfr_set_prec (a, q);
            mpfr_ui_pow_ui (a, base, p, rnd1);
            if (div)
              {
                mpfr_set_ui (b, 1, rnd_mode);
                mpfr_div (a, b, a, rnd_mode);
              }
            /* now a is an approximation by default of 1/base^(f-n) */
            mpfr_mul (b, op, a, rnd_mode);
          }
      }

    if (neg)
      MPFR_CHANGE_SIGN(b); /* put b positive */

    if (q > 2 * prec + BITS_PER_MP_LIMB)
      {
        /* if the intermediate precision exceeds twice that of the input,
           a worst-case for the division cannot occur */
        ok = 1;
        rnd_mode = GMP_RNDN;
      }
    else ok = pow2 || mpfr_can_round (b, q-err, rnd_mode, rnd_mode, prec);

  } while (ok==0 && (q+=BITS_PER_MP_LIMB) );

  if (ok)
    {
      mp_rnd_t rnd = rnd_mode;

      if (neg)
        switch (rnd_mode)
          {
          case GMP_RNDU: rnd = GMP_RNDZ; break;
          case GMP_RNDD: rnd = GMP_RNDU; break;
          }

      mpfr_round (b, rnd, MPFR_EXP(b));
    }

  prec = MPFR_EXP(b); /* may have changed due to rounding */

  /* now the mantissa is the integer part of b */
  q = 1 + (prec - 1) / BITS_PER_MP_LIMB;
  _mpz_realloc (bz, q);
  sh = prec % BITS_PER_MP_LIMB;
  e = 1 + (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB - q;
  if (sh)
    mpn_rshift (PTR(bz), MPFR_MANT(b) + e, q, BITS_PER_MP_LIMB - sh);
  else
    MPN_COPY(PTR(bz), MPFR_MANT(b) + e, q);
  bz->_mp_size = q;

  /* computes the number of characters needed */
  q = neg + n + 2; /* n+1 may not be enough for 100000... */
  if (str0 == NULL)
    str0 = (*__gmp_allocate_func) (q);

  str = str0; /* restore initial value in case we had to restart */

  if (neg)
    *str++ = '-';

  mpz_get_str (str, base, bz); /* n digits of mantissa */

  if (strlen(str) > n)
    {
      f++; /* possible due to rounding */
      goto start_again;
    }
  else if (strlen(str)==n-1) {
    f--;
    str[n-1]='0';
    str[n]='\0';
  }
  *expptr = f;
  mpfr_clear (a);
  mpfr_clear (b);
  mpz_clear (bz);

  /* if the given string was null, ensure we return a block which is exactly
     strlen(str)+1 bytes long (useful for __gmp_free_func and the C++ wrapper)
  */
  if (str_is_null && ((strlen(str0) + 1) != q))
    str0 = (char*) (*__gmp_reallocate_func) (str0, q, strlen(str0) + 1);

  return str0;
}