Copyright 1999-2020 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. The GNU 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 3 of the License, or (at your option) any later version. The GNU 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 GNU MPFR Library; see the file COPYING.LESSER. If not, see https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. Table of contents: 1. Documentation 2. Compiler/library detection 3. Changes in existing functions 4. New functions to implement 5. Efficiency 6. Miscellaneous 7. Portability ############################################################################## 1. Documentation ############################################################################## - add a description of the algorithms used and a proof of correctness ############################################################################## 2. Compiler/library detection ############################################################################## - update ICC detection. * Use only __INTEL_COMPILER instead of the obsolete macro __ICC? ############################################################################## 3. Changes in existing functions ############################################################################## - export mpfr_overflow and mpfr_underflow as public functions - many functions currently taking into account the precision of the *input* variable to set the initial working precision (acosh, asinh, cosh, ...). This is nonsense since the "average" working precision should only depend on the precision of the *output* variable (and maybe on the *value* of the input in case of cancellation). -> remove those dependencies from the input precision. - mpfr_can_round: change the meaning of the 2nd argument (err). Currently the error is at most 2^(MPFR_EXP(b)-err), i.e. err is the relative shift wrt the most significant bit of the approximation. I propose that the error is now at most 2^err ulps of the approximation, i.e. 2^(MPFR_EXP(b)-MPFR_PREC(b)+err). - mpfr_set_q first tries to convert the numerator and the denominator to mpfr_t. But this conversion may fail even if the correctly rounded result is representable. New way to implement: Function q = a/b. nq = PREC(q) na = PREC(a) nb = PREC(b) If na < nb a <- a*2^(nb-na) n <- na-nb+ (HIGH(a,nb) >= b) if (n >= nq) bb <- b*2^(n-nq) a = q*bb+r --> q has exactly n bits. else aa <- a*2^(nq-n) aa = q*b+r --> q has exactly n bits. If RNDN, takes nq+1 bits. (See also the new division function). - revisit the conversion functions between a MPFR number and a native floating-point value. * Consequences if some exception is trapped? * Specify under which conditions (current rounding direction and precision of the FPU, whether a format has been recognized...), correct rounding is guaranteed. Fix the code if need be. Do not forget subnormals. * Provide mpfr_buildopt_* functions to tell whether the format of a native type (float / double / long double) has been recognized and which format it is? * For functions that return a native floating-point value (mpfr_get_flt, mpfr_get_d, mpfr_get_ld, mpfr_get_decimal64), in case of underflow or overflow, follow the convention used for the functions in ? See §7.12.1 "Treatment of error conditions" of ISO C11, which provides two ways of handling error conditions, depending on math_errhandling: errno (to be set to ERANGE here) and floating-point exceptions. If floating-point exceptions need to be generated, do not use feraiseexcept(), as this function may require the math library (-lm); use a floating-point expression instead, such as DBL_MIN * DBL_MIN (underflow) or DBL_MAX * DBL_MAX (overflow), which are probably safe as used in the GNU libc implementation. * For testing the lack of subnormal support: see the -mfpu GCC option for ARM and https://en.wikipedia.org/wiki/Denormal_number#Disabling_denormal_floats_at_the_code_level ############################################################################## 4. New functions to implement ############################################################################## - a function to compute the hash of a floating-point number (suggested by Patrick Pelissier) - implement new functions from the C++17 standard: http://en.cppreference.com/w/cpp/numeric/special_math assoc_laguerre, assoc_legendre, comp_ellint_1, comp_ellint_2, comp_ellint_3, cyl_bessel_i, cyl_bessel_j, cyl_bessel_k, cyl_neumann, ellint_1, ellint_2, ellint_3, hermite, legendre, laguerre, sph_bessel, sph_legendre, sph_neumann. Already in mpfr4: beta and riemann_zeta. See also https://isocpp.org/files/papers/P0226R1.pdf and §29.9.5 in the C++17 draft: https://github.com/cplusplus/draft/blob/master/source/numerics.tex - implement mpfr_q_sub, mpfr_z_div, mpfr_q_div? - implement mpfr_pow_q and variants with two integers (native or mpz) instead of a rational? See IEEE P1788. - implement functions for random distributions, see for example https://sympa.inria.fr/sympa/arc/mpfr/2010-01/msg00034.html (suggested by Charles Karney , 18 Jan 2010): * a Bernoulli distribution with prob p/q (exact) * a general discrete distribution (i with prob w[i]/sum(w[i]) (Walker algorithm, but make it exact) * a uniform distribution in (a,b) * exponential distribution (mean lambda) (von Neumann's method?) * normal distribution (mean m, s.d. sigma) (ratio method?) - wanted for Magma [John Cannon , Tue, 19 Apr 2005]: HypergeometricU(a,b,s) = 1/gamma(a)*int(exp(-su)*u^(a-1)*(1+u)^(b-a-1), u=0..infinity) JacobiThetaNullK PolylogP, PolylogD, PolylogDold: see http://arxiv.org/abs/math.CA/0702243 and the references herein. JBessel(n, x) = BesselJ(n+1/2, x) KBessel, KBessel2 [2nd kind] JacobiTheta (see http://www.ams.org/journals/mcom/0000-000-00/S0025-5718-2017-03245-2/home.html) LogIntegral ExponentialIntegralEn (formula 5.1.4 of Abramowitz and Stegun) DawsonIntegral GammaD(x) = Gamma(x+1/2) - new functions of IEEE 754-2008, and more generally functions of the C binding draft TS 18661-4: http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1946.pdf Some propositions about rootn: mpfr_rootn_si, mpfr_rootn_sj, mpfr_rootn_z, and versions with an unsigned integer: mpfr_rootn_ui (now implemented, as similar to mpfr_root) and mpfr_rootn_uj. - functions defined in the LIA-2 standard + minimum and maximum (5.2.2): max, min, max_seq, min_seq, mmax_seq and mmin_seq (mpfr_min and mpfr_max correspond to mmin and mmax); + rounding_rest, floor_rest, ceiling_rest (5.2.4); + remr (5.2.5): x - round(x/y) y; + error functions from 5.2.7 (if useful in MPFR); + power1pm1 (5.3.6.7): (1 + x)^y - 1; + logbase (5.3.6.12): \log_x(y); + logbase1p1p (5.3.6.13): \log_{1+x}(1+y); + rad (5.3.9.1): x - round(x / (2 pi)) 2 pi = remr(x, 2 pi); + axis_rad (5.3.9.1) if useful in MPFR; + cycle (5.3.10.1): rad(2 pi x / u) u / (2 pi) = remr(x, u); + axis_cycle (5.3.10.1) if useful in MPFR; + sinu, cosu, tanu, cotu, secu, cscu, cossinu, arcsinu, arccosu, arctanu, arccotu, arcsecu, arccscu (5.3.10.{2..14}): sin(x 2 pi / u), etc.; [from which sinpi(x) = sin(Pi*x), ... are trivial to implement, with u=2.] + arcu (5.3.10.15): arctan2(y,x) u / (2 pi); + rad_to_cycle, cycle_to_rad, cycle_to_cycle (5.3.11.{1..3}). - From GSL, missing special functions (if useful in MPFR): (cf https://www.gnu.org/software/gsl/manual/gsl-ref.html#Special-Functions) + The Airy functions Ai(x) and Bi(x) defined by the integral representations: * Ai(x) = (1/\pi) \int_0^\infty \cos((1/3) t^3 + xt) dt * Bi(x) = (1/\pi) \int_0^\infty (e^(-(1/3) t^3) + \sin((1/3) t^3 + xt)) dt * Derivatives of Airy Functions + The Bessel functions for n integer and n fractional: * Regular Modified Cylindrical Bessel Functions I_n * Irregular Modified Cylindrical Bessel Functions K_n * Regular Spherical Bessel Functions j_n: j_0(x) = \sin(x)/x, j_1(x)= (\sin(x)/x-\cos(x))/x & j_2(x)= ((3/x^2-1)\sin(x)-3\cos(x)/x)/x Note: the "spherical" Bessel functions are solutions of x^2 y'' + 2 x y' + [x^2 - n (n+1)] y = 0 and satisfy j_n(x) = sqrt(Pi/(2x)) J_{n+1/2}(x). They should not be mixed with the classical Bessel Functions, also noted j0, j1, jn, y0, y1, yn in C99 and mpfr. Cf https://en.wikipedia.org/wiki/Bessel_function#Spherical_Bessel_functions *Irregular Spherical Bessel Functions y_n: y_0(x) = -\cos(x)/x, y_1(x)= -(\cos(x)/x+\sin(x))/x & y_2(x)= (-3/x^3+1/x)\cos(x)-(3/x^2)\sin(x) * Regular Modified Spherical Bessel Functions i_n: i_l(x) = \sqrt{\pi/(2x)} I_{l+1/2}(x) * Irregular Modified Spherical Bessel Functions: k_l(x) = \sqrt{\pi/(2x)} K_{l+1/2}(x). + Clausen Function: Cl_2(x) = - \int_0^x dt \log(2 \sin(t/2)) Cl_2(\theta) = \Im Li_2(\exp(i \theta)) (dilogarithm). + Dawson Function: \exp(-x^2) \int_0^x dt \exp(t^2). + Debye Functions: D_n(x) = n/x^n \int_0^x dt (t^n/(e^t - 1)) + Elliptic Integrals: * Definition of Legendre Forms: F(\phi,k) = \int_0^\phi dt 1/\sqrt((1 - k^2 \sin^2(t))) E(\phi,k) = \int_0^\phi dt \sqrt((1 - k^2 \sin^2(t))) P(\phi,k,n) = \int_0^\phi dt 1/((1 + n \sin^2(t))\sqrt(1 - k^2 \sin^2(t))) * Complete Legendre forms are denoted by K(k) = F(\pi/2, k) E(k) = E(\pi/2, k) * Definition of Carlson Forms RC(x,y) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1) RD(x,y,z) = 3/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-3/2) RF(x,y,z) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2) RJ(x,y,z,p) = 3/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2) (t+p)^(-1) + Elliptic Functions (Jacobi) + N-relative exponential: exprel_N(x) = N!/x^N (\exp(x) - \sum_{k=0}^{N-1} x^k/k!) + exponential integral: E_2(x) := \Re \int_1^\infty dt \exp(-xt)/t^2. Ei_3(x) = \int_0^x dt \exp(-t^3) for x >= 0. Ei(x) := - PV(\int_{-x}^\infty dt \exp(-t)/t) + Hyperbolic/Trigonometric Integrals Shi(x) = \int_0^x dt \sinh(t)/t Chi(x) := Re[ \gamma_E + \log(x) + \int_0^x dt (\cosh[t]-1)/t] Si(x) = \int_0^x dt \sin(t)/t Ci(x) = -\int_x^\infty dt \cos(t)/t for x > 0 AtanInt(x) = \int_0^x dt \arctan(t)/t [ \gamma_E is the Euler constant ] + Fermi-Dirac Function: F_j(x) := (1/r\Gamma(j+1)) \int_0^\infty dt (t^j / (\exp(t-x) + 1)) + Pochhammer symbol (a)_x := \Gamma(a + x)/\Gamma(a) : see [Smith01] in algorithms.bib logarithm of the Pochhammer symbol + Gegenbauer Functions + Laguerre Functions + Eta Function: \eta(s) = (1-2^{1-s}) \zeta(s) Hurwitz zeta function: \zeta(s,q) = \sum_0^\infty (k+q)^{-s}. + Lambert W Functions, W(x) are defined to be solutions of the equation: W(x) \exp(W(x)) = x. This function has multiple branches for x < 0 (2 funcs W0(x) and Wm1(x)) From Fredrik Johansson: See https://cs.uwaterloo.ca/research/tr/1993/03/W.pdf, in particular formulas 5.2 and 5.3 for the error bound: one first computes an approximation w, and then evaluates the residual w e^w - x. There is an expression for the error in terms of the residual and the derivative W'(t), where the derivative can be bounded by piecewise simple functions, something like min(1, 1/t) when t >= 0. See https://arxiv.org/abs/1705.03266 for rigorous error bounds. + Trigamma Function psi'(x). and Polygamma Function: psi^{(m)}(x) for m >= 0, x > 0. - functions from ISO/IEC 24747:2009 (Extensions to the C Library, to Support Mathematical Special Functions). Standard: http://www.iso.org/iso/catalogue_detail.htm?csnumber=38857 Draft: http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1292.pdf Rationale: http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1244.pdf See also: http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2010/n3060.pdf (similar, for C++). Also check whether the functions that are already implemented in MPFR match this standard. - from gnumeric (www.gnome.org/projects/gnumeric/doc/function-reference.html): - incomplete beta function, see message from Martin Maechler on 18 Jan 2016, and Section 6.6 in Abramowitz & Stegun - betaln - degrees - radians - sqrtpi - mpfr_inp_raw, mpfr_out_raw (cf mail "Serialization of mpfr_t" from Alexey and answer from Granlund on mpfr list, May 2007) - [maybe useful for SAGE] implement companion frac_* functions to the rint_* functions. For example mpfr_frac_floor(x) = x - floor(x). (The current mpfr_frac function corresponds to mpfr_rint_trunc.) - scaled erfc (https://sympa.inria.fr/sympa/arc/mpfr/2009-05/msg00054.html) - asec, acsc, acot, asech, acsch and acoth (mail from Björn Terelius on mpfr list, 18 June 2009) - function to reduce the precision of a variable, with a ternary value in input, i.e. taking care of double rounding. Two possible forms: like mpfr_set (i.e. with input and output) or like mpfr_prec_round (i.e. with a single variable). mpfr_subnormalize and mpfr_round_nearest_away_end could use it. - UBF functions for +, -, *, fmma, /, sqrt. Support UBF in mpfr_check_range or add mpfr_ubf_check_range? Make this available in the API, e.g. for MPC. - mpfr_cmp_uj and mpfr_cmp_sj. They would be useful to test MPFR with _MPFR_EXP_FORMAT=4. - base conversion with the round-trip property using a minimal precision, such as the to_chars functions from the C++ standard: The functions [...] ensure that the string representation consists of the smallest number of characters such that there is at least one digit before the radix point (if present) and parsing the representation using the corresponding from_chars function recovers value exactly. [Note: This guarantee applies only if to_chars and from_chars are executed on the same implementation. — end note] If there are several such representations, the representation with the smallest difference from the floating-point argument value is chosen, resolving any remaining ties using rounding according to round_to_nearest. Text from: https://www.zsh.org/mla/workers/2019/msg01138.html - Serialization / Deserialization. Suggested by Frédéric Pétrot: https://sympa.inria.fr/sympa/arc/mpfr/2020-02/msg00006.html like mpfr_fpif_{import,export}, but with memory instead of file. Idea of implementation to reuse most of the code and change very little: Instead of passing a FILE *fh, pass a struct ext_data *h, and instead of using fread and fwrite, use h->read (h, buffer, size) h->write (h, buffer, size) respectively. The struct ext_data structure could contain the following fields: * read: pointer to a wrapper function for the read method. * write: pointer to a wrapper function for the write method. * FILE *fh: to be used for operations with files. * unsigned char *arena: to be used for operations with memory. The wrapper functions for the read method could be: static int read_from_file (struct ext_data *h, unsigned char *buffer, size_t size) { return fread (buffer, size, 1, h->fh) != 1; } static int read_from_memory (struct ext_data *h, unsigned char *buffer, size_t size) { if (h->arena == NULL) return 1; memcpy (buffer, h->arena, size); h->arena += size; return 0; } So I expect very few changes in the existing code: * Write a few wrapper functions. * Rename mpfr_fpif_export to mpfr_fpif_export_aux and mpfr_fpif_import to mpfr_fpif_import_aux. * In the existing functions, replace FILE *fh, and fread/fwrite calls as mentioned above. * Add new mpfr_fpif_export, mpfr_fpif_import, mpfr_fpif_export_mem, mpfr_fpif_import_mem. ############################################################################## 5. Efficiency ############################################################################## - Fredrik Johansson reports that mpfr_ai is slow for large arguments: an asymptotic expansion should be used (once done, remove REDUCE_EMAX from tests/tai.c and update the description in mpfr.texi). - for exp(x), Fredrik Johansson reports a 20% speed improvement starting from 4000 bits, and up to a 75% memory improvement in his Arb implementation, by using recursive instead of iterative binary splitting: https://github.com/fredrik-johansson/arb/blob/master/elefun/exp_sum_bs_powtab.c - improve mpfr_grandom using the algorithm in http://arxiv.org/abs/1303.6257 - implement a mpfr_sqrthigh algorithm based on Mulders' algorithm, with a basecase variant - use mpn_div_q to speed up mpfr_div. However mpn_div_q, which is new in GMP 5, is not documented in the GMP manual, thus we are not sure it guarantees to return the same quotient as mpn_tdiv_qr. Also mpfr_div uses the remainder computed by mpn_divrem. A workaround would be to first try with mpn_div_q, and if we cannot (easily) compute the rounding, then use the current code with mpn_divrem. - improve atanh(x) for small x by using atanh(x) = log1p(2x/(1-x)), and log1p should also be improved for small arguments. - compute exp by using the series for cosh or sinh, which has half the terms (see Exercise 4.11 from Modern Computer Arithmetic, version 0.3) The same method can be used for log, using the series for atanh, i.e., atanh(x) = 1/2*log((1+x)/(1-x)). - improve mpfr_gamma (see https://code.google.com/p/fastfunlib/). A possible idea is to implement a fast algorithm for the argument reconstruction gamma(x+k): instead of performing k products by x+i, we could precompute x^2, ..., x^m for m ~ sqrt(k), and perform only sqrt(k) products. One could also use the series for 1/gamma(x), see for example http://dlmf.nist.gov/5/7/ or formula (36) from http://mathworld.wolfram.com/GammaFunction.html - improve the computation of Bernoulli numbers: instead of computing just one B[2n] at a time in mpfr_bernoulli_internal, we could compute several at a time, sharing the expensive computation of the 1/p^(2n) series. - fix regression with mpfr_mpz_root (from Keith Briggs, 5 July 2006), for example on 3Ghz P4 with gmp-4.2, x=12.345: prec=50000 k=2 k=3 k=10 k=100 mpz_root 0.036 0.072 0.476 7.628 mpfr_mpz_root 0.004 0.004 0.036 12.20 See also mail from Carl Witty on mpfr list, 09 Oct 2007. - for sparse input (say x=1 with 2 bits), mpfr_exp is not faster than for full precision when precision <= MPFR_EXP_THRESHOLD. The reason is that argument reduction kills sparsity. Maybe avoid argument reduction for sparse input? - speed up mpfr_atan for large arguments (to speed up mpc_log) see FR #6198 - improve mpfr_sin on values like ~pi (do not compute sin from cos, because of the cancellation). For instance, reduce the input modulo pi/2 in [-pi/4,pi/4], and define auxiliary functions for which the argument is assumed to be already reduced (so that the sin function can avoid unnecessary computations by calling the auxiliary cos function instead of the full cos function). This will require a native code for sin, for example using the reduction sin(3x)=3sin(x)-4sin(x)^3. See https://sympa.inria.fr/sympa/arc/mpfr/2007-08/msg00001.html and the following messages. - improve generic.c to work for number of terms <> 2^k - rewrite mpfr_greater_p... as native code. - mpf_t uses a scheme where the number of limbs actually present can be less than the selected precision, thereby allowing low precision values (for instance small integers) to be stored and manipulated in an mpf_t efficiently. Perhaps mpfr should get something similar, especially if looking to replace mpf with mpfr, though it'd be a major change. Alternately perhaps those mpfr routines like mpfr_mul where optimizations are possible through stripping low zero bits or limbs could check for that (this would be less efficient but easier). - try the idea of the paper "Reduced Cancellation in the Evaluation of Entire Functions and Applications to the Error Function" by W. Gawronski, J. Mueller and M. Reinhard, to be published in SIAM Journal on Numerical Analysis: to avoid cancellation in say erfc(x) for x large, they compute the Taylor expansion of erfc(x)*exp(x^2/2) instead (which has less cancellation), and then divide by exp(x^2/2) (which is simpler to compute). - replace the *_THRESHOLD macros by global (TLS) variables that can be changed at run time (via a function, like other variables)? One benefit is that users could use a single MPFR binary on several machines (e.g., a library provided by binary packages or shared via NFS) with different thresholds. On the default values, this would be a bit less efficient than the current code, but this isn't probably noticeable (this should be tested). Something like: long *mpfr_tune_get(void) to get the current values (the first value is the size of the array). int mpfr_tune_set(long *array) to set the tune values. int mpfr_tune_run(long level) to find the best values (the support for this feature is optional, this can also be done with an external function). - better distinguish different processors (for example Opteron and Core 2) and use corresponding default tuning parameters (as in GMP). This could be done in configure.ac to avoid hacking config.guess, for example define MPFR_HAVE_CORE2. Note (VL): the effect on cross-compilation (that can be a processor with the same architecture, e.g. compilation on a Core 2 for an Opteron) is not clear. The choice should be consistent with the build target (e.g. -march or -mtune value with gcc). Also choose better default values. For instance, the default value of MPFR_MUL_THRESHOLD is 40, while the best values that have been found are between 11 and 19 for 32 bits and between 4 and 10 for 64 bits! - during the Many Digits competition, we noticed that (our implantation of) Mulders short product was slower than a full product for large sizes. This should be precisely analyzed and fixed if needed. - for various functions, check the timings as a function of the magnitude of the input (and the input and/or output precisions?), and use better thresholds for asymptotic expansions. - improve the special case of mpfr_{add,sub} (x, x, y, ...) when |x| > |y| to do the addition in-place and have a complexity of O(prec(y)) in most cases. The mpfr_{add,sub}_{d,ui} functions should automatically benefit from this change. - in gmp_op.c, for functions with mpz_srcptr, check whether mpz_fits_slong_p is really useful in all cases (see TODO in this file). - optimize code that uses a test based on the fact that x >> s is undefined in C for s == width of x but the result is expected to be 0. ARM and PowerPC could benefit from such an optimization, but not x86. This needs support from the compiler. For PowerPC: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=79233 - deal with MPFR_RNDF in mpfr_round_near_x (replaced by MPFR_RNDZ). - instead of a fixed mparam.h, optionally use function multiversioning (FMV), currently only available with the GNU C++ front end: https://gcc.gnu.org/wiki/FunctionMultiVersioning According to https://lwn.net/Articles/691932/ the dispatch resolution is now done by the dynamic loader, so that this should be fast enough (the cost would be the reading of a static variable, initialized at load time, instead of a constant). In particular, binary package distributions would benefit from FMV as only one binary is generated for different processor families. ############################################################################## 6. Miscellaneous ############################################################################## - [suggested by Tobias Burnus and Asher Langton , Wed, 01 Aug 2007] support quiet and signaling NaNs in mpfr: * functions to set/test a quiet/signaling NaN: mpfr_set_snan, mpfr_snan_p, mpfr_set_qnan, mpfr_qnan_p * correctly convert to/from double (if encoding of s/qNaN is fixed in 754R) Note: Signaling NaNs are not specified by the ISO C standard and may not be supported by the implementation. GCC needs the -fsignaling-nans option (but this does not affect the C library, which may or may not accept signaling NaNs). - check the constants mpfr_set_emin (-16382-63) and mpfr_set_emax (16383) in get_ld.c and the other constants, and provide a testcase for large and small numbers. - from Kevin Ryde : Also for pi.c, a pre-calculated compiled-in pi to a few thousand digits would be good value I think. After all, say 10000 bits using 1250 bytes would still be small compared to the code size! Store pi in round to zero mode (to recover other modes). - add other prototypes for round to nearest-away (mpfr_round_nearest_away only deals with the prototypes of say mpfr_sin) or implement it as a native rounding mode - add a new roundind mode: round to odd. If the result is not exactly representable, then round to the odd mantissa. This rounding has the nice property that for k > 1, if: y = round(x, p+k, TO_ODD) z = round(y, p, TO_NEAREST_EVEN), then z = round(x, p, TO_NEAREST_EVEN) so it avoids the double-rounding problem. VL: I prefer the (original?) term "sticky rounding", as used in J Strother Moore, Tom Lynch, Matt Kaufmann. A Mechanically Checked Proof of the Correctness of the Kernel of the AMD5K86 Floating-Point Division Algorithm. IEEE Transactions on Computers, 1996. and http://www.russinoff.com/libman/text/node26.html - new rounding mode MPFR_RNDE when the result is known to be exact? * In normal mode, this would allow MPFR to optimize using this information. * In debug mode, MPFR would check that the result is exact (i.e. that the ternary value is 0). - add tests of the ternary value for constants - When doing Extensive Check (--enable-assert=full), since all the functions use a similar use of MACROS (ZivLoop, ROUND_P), it should be possible to do such a scheme: For the first call to ROUND_P when we can round. Mark it as such and save the approximated rounding value in a temporary variable. Then after, if the mark is set, check if: - we still can round. - The rounded value is the same. It should be a complement to tgeneric tests. - in div.c, try to find a case for which cy != 0 after the line cy = mpn_sub_1 (sp + k, sp + k, qsize, cy); (which should be added to the tests), e.g. by having {vp, k} = 0, or prove that this cannot happen. - add a configure test for --enable-logging to ignore the option if it cannot be supported. Modify the "configure --help" description to say "on systems that support it". - add generic bad cases for functions that don't have an inverse function that is implemented (use a single Newton iteration). - add bad cases for the internal error bound (by using a dichotomy between a bad case for the correct rounding and some input value with fewer Ziv iterations?). - add an option to use a 32-bit exponent type (int) on LP64 machines, mainly for developers, in order to be able to test the case where the extended exponent range is the same as the default exponent range, on such platforms. Tests can be done with the exp-int branch (added on 2010-12-17, and many tests fail at this time). - test underflow/overflow detection of various functions (in particular mpfr_exp) in reduced exponent ranges, including ranges that do not contain 0. - add an internal macro that does the equivalent of the following? MPFR_IS_ZERO(x) || MPFR_GET_EXP(x) <= value - check whether __gmpfr_emin and __gmpfr_emax could be replaced by a constant (see README.dev). Also check the use of MPFR_EMIN_MIN and MPFR_EMAX_MAX. - add a test checking that no mpfr.h macros depend on mpfr-impl.h (the current tests cannot check that since mpfr-impl.h is always included). - move some macro definitions from acinclude.m4 to the m4 directory as suggested by the Automake manual? The reason is that the acinclude.m4 file is big and a bit difficult to read. - use symbol versioning. - check whether mpz_t caching (pool) is necessary. Timings with -static with details about the C / C library implementation should be put somewhere as a comment in the source or in the doc. Using -static is important because otherwise the cache saves the dynamic call to mpz_init and mpz_clear; so, what we're measuring is not clear. See thread: https://gmplib.org/list-archives/gmp-devel/2015-September/004147.html Summary: It will not be integrated in GMP because 1) This yields problems with threading (in MPFR, we have TLS variables, but this is not the case of GMP). 2) The gain (if confirmed with -static) would be due to a poor malloc implementation (timings would depend on the platform). 3) Applications would use more RAM. Additional notes [VL]: the major differences in the timings given by Patrick in 2014-01 under Linux were: Before: arccos(x) took 0.054689 ms (32767 eval in 1792 ms) arctan(x) took 0.042116 ms (32767 eval in 1380 ms) After: arccos(x) took 0.043580 ms (32767 eval in 1428 ms) arctan(x) took 0.035401 ms (32767 eval in 1160 ms) mpfr_acos doesn't use mpz, but calls mpfr_atan, so that the issue comes from mpfr_atan, which uses mpz a lot. The problem mainly comes from the reallocations in GMP because mpz_init is used instead of mpz_init2 with the estimated maximum size. Other places in the code that uses mpz_init may be concerned. Issues with mpz_t caching: * The pool can take much memory, which may no longer be useful. For instance: mpfr_init2 (x, 10000000); mpfr_log_ui (x, 17, MPFR_RNDN); /* ... */ mpfr_clear (x); /* followed by code using only small precision */ while contrary to real caches, they contain no data. This is not valuable memory: freeing/allocating a large block of memory is much faster than the actual computations, so that mpz_t caching has no impact on the performance in such cases. A pool with large blocks also potentially destroys the data locality. * It assumes that the real GMP functions are __gmpz_init and __gmpz_clear, which are not part of the official GMP API, thus is based on GMP internals, which may change in the future or may be different in forks / compatible libraries / etc. This can be solved if MPFR code calls mpfr_mpz_init / mpfr_mpz_clear directly, avoiding the #define's. Questions that need to be answered: * What about the comparisons with other memory allocators? * Shouldn't the pool be part of the memory allocator? For the default memory allocator (malloc): RFE? If it is decided to keep some form of mpz_t caching, a possible solution for both issues: define mpfr_mpz_init2 and mpfr_mpz_clear2, which both take 2 arguments like mpz_init2, where mpfr_mpz_init2 behaves in a way similar to mpz_init2, and mpfr_mpz_clear2 behaves in a way similar to mpz_clear but where the size argument is a hint for the pool; if it is too large, then the mpz_t should not be pushed back to the pool. The size argument of mpfr_mpz_init2 could also be a hint to decide which element to pull from the pool. - in tsum, add testcases for mpfr_sum triggering the bug fixed in r9722, that is, with a large error during the computation of the secondary term (when the TMD occurs). - use the keyword "static" in array indices of parameter declarations with C99 compilers (6.7.5.3p7) when the pointer is expected not to be null? For instance, if mpfr.h is changed to have: __MPFR_DECLSPEC void mpfr_dump (const __mpfr_struct [static 1]); and one calls mpfr_dump (NULL); one gets a warning with Clang. This is just an example; this needs to be done in a clean way. See: http://stackoverflow.com/a/3430353/3782797 https://hamberg.no/erlend/posts/2013-02-18-static-array-indices.html - change most mpfr_urandomb occurrences to mpfr_urandom in the tests? (The one done in r10573 allowed us to find a bug even without assertion checking.) - tzeta has been much slower since r9848 (which increases the precision of the input for the low output precisions), at least with the x86 32-bit ABI. This seems to come from the fact that the working precision in the mpfr_zeta implementation depends on the precision of the input. Once mpfr_zeta has improved, change the last argument of test_generic in tzeta.c back to 5 (as it was before r10667). - check the small-precision tables in the tests? This may require to export some pointer to the tables, but this could be done only if some debug macro is defined. - optionally use malloc() for the caches? See mpfr_mp_memory_cleanup. Note: This can be implemented by adding a TLS flag saying whether we are under cache generation or not, and by making the MPFR allocation functions consider this flag. Moreover, this can only work for mpfr_t caching (floating-point constants), not for mpz_t caching (Bernoulli constants) because we do not have the control of memory allocation for mpz_init. - use GCC's nonnull attribute (available since GCC 4.0) where applicable. - avoid the use of MPFR_MANT(x) as an lvalue; use other (more high level) internal macros if possible, such as MPFR_TMP_INIT1, MPFR_TMP_INIT and MPFR_ALIAS. ############################################################################## 7. Portability ############################################################################## - add a web page with results of builds on different architectures - [Kevin about texp.c long strings] For strings longer than c99 guarantees, it might be cleaner to introduce a "tests_strdupcat" or something to concatenate literal strings into newly allocated memory. I thought I'd done that in a couple of places already. Arrays of chars are not much fun. - use https://gcc.gnu.org/viewcvs/gcc/trunk/config/stdint.m4 for mpfr-gmp.h - By default, GNU Automake adds -I options to local directories, with the side effect that these directories have the precedence to search for system headers (#include <...>). This may make the build fail if a C implementation includes a file that has the same name as one used in such a directory. For instance, if one adds an empty file "src/bits/types.h", then the MPFR build fails under Linux because /usr/include/stdio.h has #include Possible workaround: * disable the default -I options with nostdinc as documented in the Automake manual; * have a rule that copies the needed files ("mpfr.h" or they should be prefixed with "mpfr-") to $(top_builddir)/include; * use "-I$(top_builddir)/include".