00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include <basix/double.hpp>
00014 #include <basix/port.hpp>
00015 #include <numerix/mmx_mpfr.hpp>
00016 #include <numerix/integer.hpp>
00017 #include <numerix/string_scnot.hpp>
00018 namespace mmx {
00019
00020 xnat mmx_significant_digits = 0;
00021 xnat mmx_bit_precision = 64;
00022 mpfr_rnd_t mmx_rounding_mode = GMP_RNDN;
00023 xnat mmx_discrepancy = (BITS_PER_LIMB >> 1);
00024 bool mmx_pretty_exponents = false;
00025
00026 double
00027 mpfr_get_magnitude (const mpfr_t arg) {
00028 if (mpfr_inf_p (arg) || !mpfr_number_p (arg))
00029 return (double) mpfr_get_emax ();
00030 if (mpfr_sgn (arg) == 0)
00031 return (double) mpfr_get_emin ();
00032 mp_exp_t expo= mpfr_get_exp (arg);
00033 mpfr_t aux;
00034 mpfr_init2 (aux, mpfr_get_prec (arg));
00035 mpfr_div_2si (aux, arg, expo, GMP_RNDN);
00036 mpfr_abs (aux, aux, GMP_RNDN);
00037 return ((double) expo) + log (mpfr_get_d (aux, GMP_RNDU)) / log (2.0);
00038
00039 }
00040
00041 string
00042 mpfr_to_string (const mpz_t x) {
00043 char* s= mpz_get_str (NULL, 10, x);
00044 string r= s;
00045 mpfr_free_str (s);
00046 return r;
00047 }
00048
00049 string
00050 mpfr_to_string (const mpfr_t x, xnat bits) {
00051 string r;
00052 if (mpfr_nan_p (x) != 0) return "NaN";
00053 if (mpfr_inf_p (x) != 0) {
00054 if (mpfr_sgn (x) < 0) return "-Infty";
00055 return "Infty";
00056 }
00057 if (mpfr_zero_p (x) != 0) return "0";
00058 xnat digs= bits==0? 0: max (((int) (0.301029995 * bits)) - 1, 2);
00059 mp_exp_t exponent;
00060 char* mantissa_s= mpfr_get_str (NULL, &exponent, 10, digs, x, GMP_RNDN);
00061 int first= 0;
00062 if (mantissa_s[0] == '-') { r << "-"; first++; }
00063 r << mantissa_s[first];
00064 r << ".";
00065 r << (mantissa_s + (first + 1));
00066 r << "e" << as_string (exponent-1);
00067 mpfr_free_str (mantissa_s);
00068 return trunc_digits (r, 0);
00069 }
00070
00071 string
00072 mpfr_double_as_string (double x) {
00073 mpfr_t aux;
00074 mpfr_init2 (aux, 64);
00075 mpfr_set_d (aux, x, GMP_RNDN);
00076 string r= mpfr_to_string (aux, 50);
00077 mpfr_clear (aux);
00078 return trunc_digits (r, 14);
00079 }
00080
00081 string
00082 zero_to_string (const mpfr_t err) {
00083 string r;
00084 mp_exp_t exponent;
00085 char* mantissa_s= mpfr_get_str (NULL, &exponent, 10, 2, err, GMP_RNDN);
00086 r << "0.0e" << as_string (exponent-1);
00087 mpfr_free_str (mantissa_s);
00088 return r;
00089 }
00090
00091 void
00092 mpfr_binary_write (const port& p, const mpfr_t arg) {
00093 binary_write<mpfr_prec_t> (p, arg->_mpfr_prec);
00094 binary_write<mpfr_sign_t> (p, arg->_mpfr_sign);
00095 binary_write<mp_exp_t > (p, arg->_mpfr_exp );
00096 xnat n= (arg->_mpfr_prec + BITS_PER_LIMB - 1) / BITS_PER_LIMB;
00097 for (xnat i=0; i<n; i++)
00098 binary_write<mp_limb_t> (p, arg->_mpfr_d[i]);
00099 }
00100
00101 void
00102 mpfr_binary_read (const port& p, mpfr_t dest) {
00103 mpfr_prec_t pr= binary_read<mpfr_prec_t> (p);
00104 mpfr_sign_t si= binary_read<mpfr_sign_t> (p);
00105 mp_exp_t ex= binary_read<mp_exp_t > (p);
00106 mpfr_set_prec (dest, pr);
00107 xnat n= (pr + BITS_PER_LIMB - 1) / BITS_PER_LIMB;
00108 for (xnat i=0; i<n; i++)
00109 dest->_mpfr_d[i]= binary_read<mp_limb_t> (p);
00110 dest->_mpfr_prec= pr;
00111 dest->_mpfr_sign= si;
00112 dest->_mpfr_exp = ex;
00113 }
00114
00115 }