/Users/mourrain/Devel/mmx/numerix/src/mmx_mpfr.cpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : mmx_gmp.cpp
* DESCRIPTION: Additional code for interface to gmp and mpfr
* COPYRIGHT  : (C) 2003  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#include <basix/double.hpp>
#include <basix/port.hpp>
#include <numerix/mmx_mpfr.hpp>
#include <numerix/integer.hpp>
#include <numerix/string_scnot.hpp>
namespace mmx {

xnat       mmx_significant_digits = 0;
xnat       mmx_bit_precision      = 64;
mpfr_rnd_t mmx_rounding_mode      = GMP_RNDN;
xnat       mmx_discrepancy        = (BITS_PER_LIMB >> 1);
bool       mmx_pretty_exponents   = false;

double
mpfr_get_magnitude (const mpfr_t arg) {
  if (mpfr_inf_p (arg) || !mpfr_number_p (arg))
    return (double) mpfr_get_emax ();
  if (mpfr_sgn (arg) == 0)
    return (double) mpfr_get_emin ();
  mp_exp_t expo= mpfr_get_exp (arg);
  mpfr_t aux;
  mpfr_init2 (aux, mpfr_get_prec (arg));
  mpfr_div_2si (aux, arg, expo, GMP_RNDN);
  mpfr_abs (aux, aux, GMP_RNDN);
  return ((double) expo) + log (mpfr_get_d (aux, GMP_RNDU)) / log (2.0);
  //return ((double) expo) + log2 (mpfr_get_d (aux, GMP_RNDU));
}

string
mpfr_to_string (const mpz_t x) {
  char* s= mpz_get_str (NULL, 10, x);
  string r= s;
  mpfr_free_str (s);
  return r;
}

string
mpfr_to_string (const mpfr_t x, xnat bits) {
  string r;
  if (mpfr_nan_p (x) != 0) return "NaN";
  if (mpfr_inf_p (x) != 0) {
    if (mpfr_sgn (x) < 0) return "-Infty";
    return "Infty";
  }
  if (mpfr_zero_p (x) != 0) return "0";
  xnat digs= bits==0? 0: max (((int) (0.301029995 * bits)) - 1, 2);
  mp_exp_t exponent;
  char* mantissa_s= mpfr_get_str (NULL, &exponent, 10, digs, x, GMP_RNDN);
  int first= 0;
  if (mantissa_s[0] == '-') { r << "-"; first++; }
  r << mantissa_s[first];
  r << ".";
  r << (mantissa_s + (first + 1));
  r << "e" << as_string (exponent-1);
  mpfr_free_str (mantissa_s);
  return trunc_digits (r, 0);
}

string
mpfr_double_as_string (double x) {
  mpfr_t aux;
  mpfr_init2 (aux, 64);
  mpfr_set_d (aux, x, GMP_RNDN);
  string r= mpfr_to_string (aux, 50);
  mpfr_clear (aux);
  return trunc_digits (r, 14);
}

string
zero_to_string (const mpfr_t err) {
  string r;
  mp_exp_t exponent;
  char* mantissa_s= mpfr_get_str (NULL, &exponent, 10, 2, err, GMP_RNDN);
  r << "0.0e" << as_string (exponent-1);
  mpfr_free_str (mantissa_s);
  return r;
}

void
mpfr_binary_write (const port& p, const mpfr_t arg) {
  binary_write<mpfr_prec_t> (p, arg->_mpfr_prec);
  binary_write<mpfr_sign_t> (p, arg->_mpfr_sign);
  binary_write<mp_exp_t   > (p, arg->_mpfr_exp );
  xnat n= (arg->_mpfr_prec + BITS_PER_LIMB - 1) / BITS_PER_LIMB;
  for (xnat i=0; i<n; i++)
    binary_write<mp_limb_t> (p, arg->_mpfr_d[i]);
}

void
mpfr_binary_read (const port& p, mpfr_t dest) {
  mpfr_prec_t pr= binary_read<mpfr_prec_t> (p);
  mpfr_sign_t si= binary_read<mpfr_sign_t> (p);
  mp_exp_t    ex= binary_read<mp_exp_t   > (p);
  mpfr_set_prec (dest, pr);
  xnat n= (pr + BITS_PER_LIMB - 1) / BITS_PER_LIMB;
  for (xnat i=0; i<n; i++)
    dest->_mpfr_d[i]= binary_read<mp_limb_t> (p);
  dest->_mpfr_prec= pr;
  dest->_mpfr_sign= si;
  dest->_mpfr_exp = ex;
}

} // namespace mmx