diff --git a/M2/Macaulay2/d/Makefile.mpfrtoa b/M2/Macaulay2/d/Makefile.mpfrtoa new file mode 100644 index 0000000000..6bffd20bdb --- /dev/null +++ b/M2/Macaulay2/d/Makefile.mpfrtoa @@ -0,0 +1,8 @@ +LDLIBS = -lgmp -lmpfr -lm + +all: mpfrtoa-test + +mpfrtoa-test: mpfrtoa.o + +clean: + rm -f mpfrtoa-test diff --git a/M2/Macaulay2/d/mpfrtoa-test.c b/M2/Macaulay2/d/mpfrtoa-test.c new file mode 100644 index 0000000000..8b2c888f3b --- /dev/null +++ b/M2/Macaulay2/d/mpfrtoa-test.c @@ -0,0 +1,49 @@ +/* mpfrtoa-test.c -- MPFR -> ASCII test suite + +Copyright 2024 Doug Torrance + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +1. Redistributions of source code must retain the above copyright +notice, this list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright +notice, this list of conditions and the following disclaimer in the +documentation and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its +contributors may be used to endorse or promote products derived from +this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +“AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ + +#include +#include +#include "mpfrtoa.h" + +int main(void) +{ + mpfr_t x; + mpfr_exp_t exp; + char *str; + + mpfr_init(x); + mpfr_set_d(x, 123, MPFR_RNDN); + str = mpfrtoa(NULL, &exp, 10, 0, x); + printf("str = %s, exp = %ld\n", str, exp); + mpfr_clear(x); + + return 0; +} diff --git a/M2/Macaulay2/d/mpfrtoa.c b/M2/Macaulay2/d/mpfrtoa.c new file mode 100644 index 0000000000..0419dfea00 --- /dev/null +++ b/M2/Macaulay2/d/mpfrtoa.c @@ -0,0 +1,174 @@ +/* mpfrtoa.c -- MPFR -> ASCII + +Copyright 2024 Doug Torrance + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +1. Redistributions of source code must retain the above copyright +notice, this list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright +notice, this list of conditions and the following disclaimer in the +documentation and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its +contributors may be used to endorse or promote products derived from +this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +“AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ + +/* Based on Burger & Dybvig's improvements to Steel & White's Dragon4 +algorithm. + +Burger, Robert G., and R. Kent Dybvig. "Printing floating-point + numbers quickly and accurately." ACM SIGPLAN Notices 31, no. 5 (1996): + 108-116. + +Steele Jr, Guy L., and Jon L. White. "How to print floating-point + numbers accurately." In Proceedings of the ACM SIGPLAN 1990 conference + on Programming language design and implementation, pp. 112-126. 1990. +*/ + +#include +#include +#include + +#include /* REMOVE ME */ + +static const char num_to_text36[] = "0123456789abcdefghijklmnopqrstuvwxyz"; +static const char num_to_text62[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" + "abcdefghijklmnopqrstuvwxyz"; + +char *mpfrtoa(char *str, mpfr_exp_t *expptr, int base, size_t num_digits, + mpfr_t op) +{ + const char *num_to_text; + int is_negative, i, condition; + mpfr_exp_t e, k; + mpfr_prec_t p; + mpz_t f, r, s, m_plus, m_minus, tmp; + mpfr_t v; + + /* if NaN, infinity, or zero, fall back to mpfr_get_str */ + if (!mpfr_number_p(op) || mpfr_zero_p(op)) + return mpfr_get_str(str, expptr, base, num_digits, op, MPFR_RNDN); + + is_negative = mpfr_sgn(op) < 0; + p = mpfr_get_prec(op); + + num_to_text = (2 <= base && base <= 36) ? num_to_text36 : num_to_text62; + base = (base > 0) ? base : -base; + + if (num_digits == 0) + num_digits = mpfr_get_str_ndigits(base, p); + + if (str == NULL) { + void *(*alloc_function)(size_t); + + mp_get_memory_functions(&alloc_function, NULL, NULL); + str = alloc_function(is_negative + num_digits + 1); + } + + mpz_inits(f, r, s, m_plus, m_minus, tmp, NULL); + + /* deal with negatives */ + i = 0; + if (is_negative) + str[i++] = '-'; + mpfr_init2(v, p); + mpfr_abs(v, op, MPFR_RNDN); /* rounding mode ignored since same precision */ + + e = mpfr_get_z_2exp(f, v); + mpz_ui_pow_ui(tmp, 2, p - 1); + + /* burger/dybvig table 1 */ + if (e >= 0) { + if (mpz_cmp(f, tmp) == 0) { + mpz_mul_2exp(r, f, e + 2); + mpz_set_ui(s, 4); + mpz_ui_pow_ui(m_plus, 2, e + 1); + mpz_ui_pow_ui(m_minus, 2, e); + } else { + mpz_mul_2exp(r, f, e + 1); + mpz_set_ui(s, 2); + mpz_ui_pow_ui(m_plus, 2, e); + mpz_set(m_minus, m_plus); + } + } else { + if ((e > mpfr_get_emin()) && (mpz_cmp(f, tmp) == 0)) { + mpz_mul_2exp(r, f, 2); + mpz_ui_pow_ui(s, 2, 2 - e); + mpz_set_ui(m_plus, 2); + mpz_set_ui(m_minus, 1); + } else { + mpz_mul_2exp(r, f, 1); + mpz_ui_pow_ui(s, 2, 1 - e); + mpz_set_ui(m_minus, 1); + mpz_set_ui(m_plus, 1); + } + } + + /* TODO: simplify this? see burger/dybvig section 3.2 */ + /* only works if e + p - 1< 10^14 */ + k = ceil((e + p - 1)*(1/log2(base))); + *expptr = k; + printf("k = %ld\n", k); + + if (k >= 0) { + mpz_ui_pow_ui(tmp, 10, k); + mpz_mul(s, s, tmp); + } else { + mpz_ui_pow_ui(tmp, 10, -k); + mpz_mul(r, r, tmp); + mpz_mul(m_plus, m_plus, tmp); + mpz_mul(m_minus, m_minus, tmp); + } + + for(;; i++) { + mpz_mul_ui(r, r, 10); + mpz_fdiv_qr(tmp, r, r, s); + str[i] = num_to_text[mpz_get_ui(tmp)]; + 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, r, m_plus); + if (mpz_cmp(tmp, 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, NULL); + + return str; +} diff --git a/M2/Macaulay2/d/mpfrtoa.h b/M2/Macaulay2/d/mpfrtoa.h new file mode 100644 index 0000000000..ddb0baafff --- /dev/null +++ b/M2/Macaulay2/d/mpfrtoa.h @@ -0,0 +1,38 @@ +/* mpfrtoa.h -- MPFR -> ASCII + +Copyright 2024 Doug Torrance + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +1. Redistributions of source code must retain the above copyright +notice, this list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright +notice, this list of conditions and the following disclaimer in the +documentation and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its +contributors may be used to endorse or promote products derived from +this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +“AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ + +#ifndef __MPFRTOA_H +#define __MPFRTOA_H + +char *mpfrtoa(char *str, mpfr_exp_t *expptr, int base, size_t num_digits, + mpfr_t op); + +#endif