diff --git a/M2/Macaulay2/d/gmp1.d b/M2/Macaulay2/d/gmp1.d index fa73c99cc2..1452ec1a3f 100644 --- a/M2/Macaulay2/d/gmp1.d +++ b/M2/Macaulay2/d/gmp1.d @@ -1,6 +1,6 @@ -- this file is small, so the exported definitions can be imported into the engine -header ""; +header "#include \"gmp_aux.h\""; --This file contains gmp related functions. --Functions in this file may make calls to stdio. @@ -9,6 +9,7 @@ use gmp; use stdio; use err; + pow(x:ZZmutable, y:ZZ, n:ulong) ::= Ccode( void, "mpz_pow_ui(", x, ",", y, ",", n, ")" ); export (x:ZZ) ^ (n:int) : ZZ := ( @@ -62,6 +63,13 @@ getstr(returnexponent:long, base:int, sigdigs:int, x:RRi) ::= ( Ccode(void, "mpfr_free_str(", strptr, ")"); ret); +getstr_dragon4(returnexponent:long, x:RR) ::= ( + strptr := Ccode(charstarOrNull, + "mpfr_dragon4(&", returnexponent, ", ", x, ")"); + ret := tostring(strptr); + Ccode(void, "mpfr_free_str(", strptr, ")"); + ret); + export format( s:int, -- number of significant digits (0 means all) ac:int, -- accuracy, how far to right of point to go (-1 means all) @@ -83,7 +91,10 @@ export format( ); if x === 0 then return array(string)(if ng then "-" else "","0"); ex := long(0); - mantissa := getstr(ex, base, s, x); + mantissa := ( + if s == meaningful + then getstr_dragon4(ex, x) + else getstr(ex, base, s, x)); nt := 0; for i from length(mantissa)-1 to 0 by -1 do ( if mantissa.i != '0' then break; diff --git a/M2/Macaulay2/d/gmp_aux.c b/M2/Macaulay2/d/gmp_aux.c index 5aebc5c67b..2c8b478c55 100644 --- a/M2/Macaulay2/d/gmp_aux.c +++ b/M2/Macaulay2/d/gmp_aux.c @@ -42,6 +42,112 @@ void mp_free_str(char *str){ } +char *mpfr_dragon4(mpfr_exp_t *expptr, mpfr_srcptr x) +{ + int i, condition; + mpfr_exp_t e, k; + mpfr_prec_t p; + mpz_t f, r, s, m_plus, m_minus, one, tmp_zz; + mpfr_t tmp_rr; + char *str; + void *(*alloc_function)(size_t); + + mp_get_memory_functions(&alloc_function, NULL, NULL); + str = alloc_function(mpfr_get_str_ndigits(10, mpfr_get_prec(x)) + 1); + + mpz_inits(f, r, s, m_plus, m_minus, one, tmp_zz, NULL); + mpfr_init(tmp_rr); + + p = mpfr_get_prec(x); + e = mpfr_get_z_2exp(f, x); + mpz_set_ui(one, 1); + mpz_mul_2exp(tmp_zz, one, p - 1); + + /* burger/dybvig table 1 */ + if (e >= 0) { + if (mpz_cmp(f, tmp_zz) == 0) { + mpz_mul_2exp(r, f, e + 2); + mpz_set_ui(s, 4); + mpz_mul_2exp(m_plus, one, e + 1); + mpz_mul_2exp(m_minus, one, e); + } else { + mpz_mul_2exp(r, f, e + 1); + mpz_set_ui(s, 2); + mpz_mul_2exp(m_plus, one, e); + mpz_set(m_minus, m_plus); + } + } else { + if ((e > mpfr_get_emin()) && (mpz_cmp(f, tmp_zz) == 0)) { + mpz_mul_2exp(r, f, 2); + mpz_mul_2exp(s, one, 2 - e); + mpz_set_ui(m_plus, 2); + mpz_set(m_minus, one); + } else { + mpz_mul_2exp(r, f, 1); + mpz_mul_2exp(s, one, 1 - e); + mpz_set(m_minus, one); + mpz_set(m_plus, one); + } + } + + /* TODO: simplify this? see burger/dybvig section 3.2 */ + mpfr_set_z(tmp_rr, r, MPFR_RNDN); + mpfr_add_z(tmp_rr, tmp_rr, m_plus, MPFR_RNDN); + mpfr_div_z(tmp_rr, tmp_rr, s, MPFR_RNDN); + mpfr_log(tmp_rr, tmp_rr, MPFR_RNDN); + mpfr_div_d(tmp_rr, tmp_rr, 2.3025850929940459 /* log 10 */, MPFR_RNDN); + mpfr_ceil(tmp_rr, tmp_rr); + k = mpfr_get_si(tmp_rr, MPFR_RNDN); + *expptr = k; + + if (k >= 0) { + mpz_ui_pow_ui(tmp_zz, 10, k); + mpz_mul(s, s, tmp_zz); + } else { + mpz_ui_pow_ui(tmp_zz, 10, -k); + mpz_mul(r, r, tmp_zz); + mpz_mul(m_plus, m_plus, tmp_zz); + mpz_mul(m_minus, m_minus, tmp_zz); + } + + for(i = 0; ; i++) { + mpz_mul_ui(r, r, 10); + mpz_fdiv_qr(tmp_zz, r, r, s); + str[i] = '0' + mpz_get_ui(tmp_zz); + mpz_mul_ui(m_plus, m_plus, 10); + mpz_mul_ui(m_minus, m_minus, 10); + + condition = 0; + + /* condition (1) */ + if (mpz_cmp(r, m_minus) < 0) + condition |= 1; + + /* condition (2) */ + mpz_add(tmp_zz, r, m_plus); + if (mpz_cmp(tmp_zz, s) > 0) /* typo in burger/dybvig -- they use "<" */ + condition |= 2; + + if (condition == 1) { /* (1) true, (2) false */ + str[i + 1] = '\0'; + break; + } else if (condition == 2) { /* (1) false, (2) true */ + str[i]++; + str[i + 1] = '\0'; + break; + } else if (condition == 3) { /* (1) & (2) both true */ + /* TODO */ + str[i + 1] = '\0'; + break; + } + } + + mpz_clears(f, r, s, m_plus, m_minus, tmp_zz, NULL); + mpfr_clear(tmp_rr); + + return str; +} + /* Local Variables: compile-command: "echo \"make: Entering directory \\`$M2BUILDDIR/Macaulay2/d'\" && make -C $M2BUILDDIR/Macaulay2/d " diff --git a/M2/Macaulay2/d/gmp_aux.h b/M2/Macaulay2/d/gmp_aux.h index 471c4e6b9a..c05a33b721 100644 --- a/M2/Macaulay2/d/gmp_aux.h +++ b/M2/Macaulay2/d/gmp_aux.h @@ -4,3 +4,4 @@ extern int mpz_hash(mpz_srcptr x); extern int mpfr_hash(mpfr_srcptr x); extern int mpfi_hash(mpfi_srcptr x); extern void mp_free_str(char* str); +extern char *mpfr_dragon4(mpfr_exp_t *expptr, mpfr_srcptr x);