summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/TODO
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/TODO')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/TODO316
1 files changed, 287 insertions, 29 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/TODO b/Build/source/libs/mpfr/mpfr-src/TODO
index d4628514f00..251ac589752 100644
--- a/Build/source/libs/mpfr/mpfr-src/TODO
+++ b/Build/source/libs/mpfr/mpfr-src/TODO
@@ -20,7 +20,7 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
Table of contents:
1. Documentation
-2. Installation
+2. Compiler/library detection
3. Changes in existing functions
4. New functions to implement
5. Efficiency
@@ -31,10 +31,10 @@ Table of contents:
1. Documentation
##############################################################################
-- add a description of the algorithms used + proof of correctness
+- add a description of the algorithms used and a proof of correctness
##############################################################################
-2. Installation
+2. Compiler/library detection
##############################################################################
- if we want to distinguish GMP and MPIR, we can check at configure time
@@ -47,6 +47,11 @@ Table of contents:
There is also a library symbol mpir_version, which should match VERSION, set
by configure, for example 0.9.0.
+- update ICC detection.
+ * Use only __INTEL_COMPILER instead of the obsolete macro __ICC?
+ * Possibly enable some features with ICC. For instance, removing
+ "&& !defined(__ICC)" from src/mpfr.h works with ICC 15.0.0 20140723.
+
##############################################################################
3. Changes in existing functions
##############################################################################
@@ -82,12 +87,58 @@ Table of contents:
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
##############################################################################
+- 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_get_decimal128 and mpfr_set_decimal128
+- implement mpfr_log_ui to compute log(n) for an unsigned long n.
+ We can write for argument reduction n = 2^k * n/2^k, where
+ 2/3 <= n/2^k < 4/3, i.e., k = floor(log2(3n))-1, thus
+ log(n) = k*log(2) + log(n/2^k), and we can use binary splitting on the
+ Taylor expansion of log(1+x) to compute log(n/2^k), where at most
+ p*log(2)/log(3) terms are needed for precision p.
+ Other idea (from Fredrik Johansson): compute log(m) + log(n/m) where
+ m=2^a*3^b*5^c*7^d and m is close to n.
- 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):
@@ -104,21 +155,19 @@ Table of contents:
PolylogP, PolylogD, PolylogDold: see http://arxiv.org/abs/math.CA/0702243
and the references herein.
JBessel(n, x) = BesselJ(n+1/2, x)
- IncompleteGamma [also wanted by <keith.briggs@bt.com> 4 Feb 2008: Gamma(a,x),
- gamma(a,x), P(a,x), Q(a,x); see A&S 6.5, ref. [Smith01] in algorithms.bib]
KBessel, KBessel2 [2nd kind]
JacobiTheta
+ (see http://www.ams.org/journals/mcom/0000-000-00/S0025-5718-2017-03245-2/home.html)
LogIntegral
- ExponentialIntegralE1
- E1(z) = int(exp(-t)/t, t=z..infinity), |arg z| < Pi
- mpfr_eint1: implement E1(x) for x > 0, and Ei(-x) for x < 0
- E1(NaN) = NaN
- E1(+Inf) = +0
- E1(-Inf) = -Inf
- E1(+0) = +Inf
- E1(-0) = -Inf
+ 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);
@@ -207,11 +256,30 @@ Table of contents:
+ 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):
- - beta
+ - 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
@@ -230,6 +298,17 @@ Table of contents:
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
+- use the src/x86_64/corei5/mparam.h file once GMP recognizes correctly the
+ Core i5 processors (note that gcc -mtune=native gives __tune_corei7__
+ and not __tune_corei5__ on those processors)
- 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
@@ -238,35 +317,33 @@ Table of contents:
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). One could also use the series for 1/gamma(x), see for example
+ 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.
-- implement Mulders algorithm for squaring and division
- 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 const_euler for large precision [for x=1.1, prec=16610, it takes
- 75% of the total time of eint(x)!]
-- speed up mpfr_atan for large arguments (to speed up mpc_log)
- [from Mark Watkins on Fri, 18 Mar 2005]
- Also mpfr_atan(x) seems slower (by a factor of 2) for x near from 1.
- Example on a Athlon for 10^5 bits: x=1.1 takes 3s, whereas 2.1 takes 1.8s.
- The current implementation does not give monotonous timing for the following:
- mpfr_random (x); for (i = 0; i < k; i++) mpfr_atan (y, x, MPFR_RNDN);
- for precision 300 and k=1000, we get 1070ms, and 500ms only for p=400!
+- 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
@@ -327,6 +404,37 @@ Table of contents:
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
##############################################################################
@@ -337,6 +445,10 @@ Table of contents:
* 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 again coverage: on 2007-07-27, Patrick Pelissier reports that the
following files are not tested at 100%: add1.c, atan.c, atan2.c,
@@ -359,8 +471,9 @@ Table of contents:
1250 bytes would still be small compared to the code size!
Store pi in round to zero mode (to recover other modes).
-- add a new rounding mode: round to nearest, with ties away from zero
- (this is roundTiesToAway in 754-2008, could be used by mpfr_round)
+- 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:
@@ -368,6 +481,36 @@ Table of contents:
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).
+
+- new "rounding mode" MPFR_RNDF (faithful rounding)?
+ This is not really a rounding mode since it is non-deterministic. The
+ goal is to avoid the Table Maker's Dilemma in internal computations.
+ The definition of faithful rounding of a real number x is: return either
+ RNDD(x) or RNDU(x). This means that if x is exactly representable, one
+ returns x exactly. In MPFR, the ternary value should be unspecified for
+ efficiency reasons.
+ Note: One typically implements faithful rounding by computing an
+ approximation to the result with some adequately chosen error bound,
+ then by rounding this approximation to nearest.
+ Concerning the choice of the error bound, if the result x is equal to
+ 1 + t, where t is a very small positive number, then the error bound
+ needs to be at most ulp(x)/4 + t. Since t can be arbitrarily small,
+ the error bound needs to be at most ulp(x)/4. And this error bound
+ is sufficient in all cases. Note that with the even rounding rule or
+ rounding away from zero, it is not needed to relax the condition when
+ x is exactly representable.
- add tests of the ternary value for constants
@@ -416,6 +559,121 @@ Table of contents:
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).
+
+- add internal or public variants of some basic functions (+, -, *) with
+ mpz_t as the exponent for correctly rounded polynomials (like fmma),
+ in order to avoid internal overflow and underflow in extreme cases?
+ Alternatively, for fmma, modify add*.c and sub*.c code to define
+ mpfr_add_raw, which takes arrays of limbs and their precision, and the
+ positive exponent delta (as mpfr_uexp_t to be always representable).
+ The exponent delta should be sufficient for now since it can be bounded
+ by MPFR_PREC_MAX+1 if need be.
+
+- 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.
+
##############################################################################
7. Portability
@@ -423,8 +681,6 @@ Table of contents:
- add a web page with results of builds on different architectures
-- support the decimal64 function without requiring --with-gmp-build
-
- [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
@@ -432,3 +688,5 @@ Table of contents:
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
+
+- with MinGW, build with -D__USE_MINGW_ANSI_STDIO by default?