diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/TODO')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/TODO | 725 |
1 files changed, 0 insertions, 725 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/TODO b/Build/source/libs/mpfr/mpfr-src/TODO deleted file mode 100644 index 684a442533d..00000000000 --- a/Build/source/libs/mpfr/mpfr-src/TODO +++ /dev/null @@ -1,725 +0,0 @@ -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 <math.h>? - 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 <ckarney@Sarnoff.com>, 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 <john@maths.usyd.edu.au>, 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 - <maechler@stat.math.ethz.ch> 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 <burnus(at)net-b.de> and - Asher Langton <langton(at)gcc.gnu.org>, 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 <user42@zip.com.au>: - 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 <bits/types.h> - 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". |