1 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
2 Contributed by the Arenaire and Cacao projects, INRIA.
4 This file is part of the GNU MPFR Library.
6 The GNU MPFR Library is free software; you can redistribute it and/or modify
7 it under the terms of the GNU Lesser General Public License as published by
8 the Free Software Foundation; either version 3 of the License, or (at your
9 option) any later version.
11 The GNU MPFR Library is distributed in the hope that it will be useful, but
12 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
13 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14 License for more details.
16 You should have received a copy of the GNU Lesser General Public License
17 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
18 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
19 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
24 3. Changes in existing functions
25 4. New functions to implement
30 ##############################################################################
32 ##############################################################################
34 - add a description of the algorithms used + proof of correctness
36 ##############################################################################
38 ##############################################################################
40 - if we want to distinguish GMP and MPIR, we can check at configure time
41 the following symbols which are only defined in MPIR:
43 #define __MPIR_VERSION 0
44 #define __MPIR_VERSION_MINOR 9
45 #define __MPIR_VERSION_PATCHLEVEL 0
47 There is also a library symbol mpir_version, which should match VERSION, set
48 by configure, for example 0.9.0.
50 ##############################################################################
51 3. Changes in existing functions
52 ##############################################################################
54 - many functions currently taking into account the precision of the *input*
55 variable to set the initial working precison (acosh, asinh, cosh, ...).
56 This is nonsense since the "average" working precision should only depend
57 on the precision of the *output* variable (and maybe on the *value* of
58 the input in case of cancellation).
59 -> remove those dependencies from the input precision.
61 - mpfr_get_str should support base up to 62 too.
64 change the meaning of the 2nd argument (err). Currently the error is
65 at most 2^(MPFR_EXP(b)-err), i.e. err is the relative shift wrt the
66 most significant bit of the approximation. I propose that the error
67 is now at most 2^err ulps of the approximation, i.e.
68 2^(MPFR_EXP(b)-MPFR_PREC(b)+err).
70 - mpfr_set_q first tries to convert the numerator and the denominator
71 to mpfr_t. But this convertion may fail even if the correctly rounded
72 result is representable. New way to implement:
73 Function q = a/b. nq = PREC(q) na = PREC(a) nb = PREC(b)
76 n <- na-nb+ (HIGH(a,nb) >= b)
79 a = q*bb+r --> q has exactly n bits.
82 aa = q*b+r --> q has exaclty n bits.
83 If RNDN, takes nq+1 bits. (See also the new division function).
86 ##############################################################################
87 4. New functions to implement
88 ##############################################################################
90 - implement mpfr_z_sub, mpfr_q_sub, mpfr_z_div, mpfr_q_div?
91 - implement functions for random distributions, see for example
92 http://websympa.loria.fr/wwsympa/arc/mpfr/2010-01/msg00034.html
93 (suggested by Charles Karney <ckarney@Sarnoff.com>, 18 Jan 2010):
94 * a Bernoulli distribution with prob p/q (exact)
95 * a general discrete distribution (i with prob w[i]/sum(w[i]) (Walker
96 algorithm, but make it exact)
97 * a uniform distribution in (a,b)
98 * exponential distribution (mean lambda) (von Neumann's method?)
99 * normal distribution (mean m, s.d. sigma) (ratio method?)
100 - wanted for Magma [John Cannon <john@maths.usyd.edu.au>, Tue, 19 Apr 2005]:
101 HypergeometricU(a,b,s) = 1/gamma(a)*int(exp(-su)*u^(a-1)*(1+u)^(b-a-1),
104 PolylogP, PolylogD, PolylogDold: see http://arxiv.org/abs/math.CA/0702243
105 and the references herein.
106 JBessel(n, x) = BesselJ(n+1/2, x)
107 IncompleteGamma [also wanted by <keith.briggs@bt.com> 4 Feb 2008: Gamma(a,x),
108 gamma(a,x), P(a,x), Q(a,x); see A&S 6.5, ref. [Smith01] in algorithms.bib]
109 KBessel, KBessel2 [2nd kind]
112 ExponentialIntegralE1
113 E1(z) = int(exp(-t)/t, t=z..infinity), |arg z| < Pi
114 mpfr_eint1: implement E1(x) for x > 0, and Ei(-x) for x < 0
121 GammaD(x) = Gamma(x+1/2)
122 - functions defined in the LIA-2 standard
123 + minimum and maximum (5.2.2): max, min, max_seq, min_seq, mmax_seq
124 and mmin_seq (mpfr_min and mpfr_max correspond to mmin and mmax);
125 + rounding_rest, floor_rest, ceiling_rest (5.2.4);
126 + remr (5.2.5): x - round(x/y) y;
127 + error functions from 5.2.7 (if useful in MPFR);
128 + power1pm1 (5.3.6.7): (1 + x)^y - 1;
129 + logbase (5.3.6.12): \log_x(y);
130 + logbase1p1p (5.3.6.13): \log_{1+x}(1+y);
131 + rad (5.3.9.1): x - round(x / (2 pi)) 2 pi = remr(x, 2 pi);
132 + axis_rad (5.3.9.1) if useful in MPFR;
133 + cycle (5.3.10.1): rad(2 pi x / u) u / (2 pi) = remr(x, u);
134 + axis_cycle (5.3.10.1) if useful in MPFR;
135 + sinu, cosu, tanu, cotu, secu, cscu, cossinu, arcsinu, arccosu,
136 arctanu, arccotu, arcsecu, arccscu (5.3.10.{2..14}):
137 sin(x 2 pi / u), etc.;
138 [from which sinpi(x) = sin(Pi*x), ... are trivial to implement, with u=2.]
139 + arcu (5.3.10.15): arctan2(y,x) u / (2 pi);
140 + rad_to_cycle, cycle_to_rad, cycle_to_cycle (5.3.11.{1..3}).
141 - From GSL, missing special functions (if useful in MPFR):
142 (cf http://www.gnu.org/software/gsl/manual/gsl-ref.html#Special-Functions)
143 + The Airy functions Ai(x) and Bi(x) defined by the integral representations:
144 * Ai(x) = (1/\pi) \int_0^\infty \cos((1/3) t^3 + xt) dt
145 * Bi(x) = (1/\pi) \int_0^\infty (e^(-(1/3) t^3) + \sin((1/3) t^3 + xt)) dt
146 * Derivatives of Airy Functions
147 + The Bessel functions for n integer and n fractional:
148 * Regular Modified Cylindrical Bessel Functions I_n
149 * Irregular Modified Cylindrical Bessel Functions K_n
150 * Regular Spherical Bessel Functions j_n: j_0(x) = \sin(x)/x,
151 j_1(x)= (\sin(x)/x-\cos(x))/x & j_2(x)= ((3/x^2-1)\sin(x)-3\cos(x)/x)/x
152 Note: the "spherical" Bessel functions are solutions of
153 x^2 y'' + 2 x y' + [x^2 - n (n+1)] y = 0 and satisfy
154 j_n(x) = sqrt(Pi/(2x)) J_{n+1/2}(x). They should not be mixed with the
155 classical Bessel Functions, also noted j0, j1, jn, y0, y1, yn in C99
157 Cf http://en.wikipedia.org/wiki/Bessel_function#Spherical_Bessel_functions
158 *Irregular Spherical Bessel Functions y_n: y_0(x) = -\cos(x)/x,
159 y_1(x)= -(\cos(x)/x+\sin(x))/x &
160 y_2(x)= (-3/x^3+1/x)\cos(x)-(3/x^2)\sin(x)
161 * Regular Modified Spherical Bessel Functions i_n:
162 i_l(x) = \sqrt{\pi/(2x)} I_{l+1/2}(x)
163 * Irregular Modified Spherical Bessel Functions:
164 k_l(x) = \sqrt{\pi/(2x)} K_{l+1/2}(x).
166 Cl_2(x) = - \int_0^x dt \log(2 \sin(t/2))
167 Cl_2(\theta) = \Im Li_2(\exp(i \theta)) (dilogarithm).
168 + Dawson Function: \exp(-x^2) \int_0^x dt \exp(t^2).
169 + Debye Functions: D_n(x) = n/x^n \int_0^x dt (t^n/(e^t - 1))
170 + Elliptic Integrals:
171 * Definition of Legendre Forms:
172 F(\phi,k) = \int_0^\phi dt 1/\sqrt((1 - k^2 \sin^2(t)))
173 E(\phi,k) = \int_0^\phi dt \sqrt((1 - k^2 \sin^2(t)))
174 P(\phi,k,n) = \int_0^\phi dt 1/((1 + n \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))
175 * Complete Legendre forms are denoted by
178 * Definition of Carlson Forms
179 RC(x,y) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1)
180 RD(x,y,z) = 3/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-3/2)
181 RF(x,y,z) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2)
182 RJ(x,y,z,p) = 3/2 \int_0^\infty dt
183 (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2) (t+p)^(-1)
184 + Elliptic Functions (Jacobi)
185 + N-relative exponential:
186 exprel_N(x) = N!/x^N (\exp(x) - \sum_{k=0}^{N-1} x^k/k!)
187 + exponential integral:
188 E_2(x) := \Re \int_1^\infty dt \exp(-xt)/t^2.
189 Ei_3(x) = \int_0^x dt \exp(-t^3) for x >= 0.
190 Ei(x) := - PV(\int_{-x}^\infty dt \exp(-t)/t)
191 + Hyperbolic/Trigonometric Integrals
192 Shi(x) = \int_0^x dt \sinh(t)/t
193 Chi(x) := Re[ \gamma_E + \log(x) + \int_0^x dt (\cosh[t]-1)/t]
194 Si(x) = \int_0^x dt \sin(t)/t
195 Ci(x) = -\int_x^\infty dt \cos(t)/t for x > 0
196 AtanInt(x) = \int_0^x dt \arctan(t)/t
197 [ \gamma_E is the Euler constant ]
198 + Fermi-Dirac Function:
199 F_j(x) := (1/r\Gamma(j+1)) \int_0^\infty dt (t^j / (\exp(t-x) + 1))
200 + Pochhammer symbol (a)_x := \Gamma(a + x)/\Gamma(a) : see [Smith01] in
202 logarithm of the Pochhammer symbol
203 + Gegenbauer Functions
205 + Eta Function: \eta(s) = (1-2^{1-s}) \zeta(s)
206 Hurwitz zeta function: \zeta(s,q) = \sum_0^\infty (k+q)^{-s}.
207 + Lambert W Functions, W(x) are defined to be solutions of the equation:
209 This function has multiple branches for x < 0 (2 funcs W0(x) and Wm1(x))
210 + Trigamma Function psi'(x).
211 and Polygamma Function: psi^{(m)}(x) for m >= 0, x > 0.
213 - from gnumeric (www.gnome.org/projects/gnumeric/doc/function-reference.html):
220 - mpfr_frexp(mpfr_t rop, mpfr_exp_t *n, mpfr_t op, mpfr_rnd_t rnd) suggested
221 by Steve Kargl <sgk@troutmask.apl.washington.edu> Sun, 7 Aug 2005
222 - mpfr_inp_raw, mpfr_out_raw (cf mail "Serialization of mpfr_t" from Alexey
223 and answer from Granlund on mpfr list, May 2007)
224 - [maybe useful for SAGE] implement companion frac_* functions to the rint_*
225 functions. For example mpfr_frac_floor(x) = x - floor(x). (The current
226 mpfr_frac function corresponds to mpfr_rint_trunc.)
227 - scaled erfc (http://websympa.loria.fr/wwsympa/arc/mpfr/2009-05/msg00054.html)
228 - asec, acsc, acot, asech, acsch and acoth (mail from Björn Terelius on mpfr
231 ##############################################################################
233 ##############################################################################
235 - compute exp by using the series for cosh or sinh, which has half the terms
236 (see Exercise 4.11 from Modern Computer Arithmetic, version 0.3)
237 The same method can be used for log, using the series for atanh, i.e.,
238 atanh(x) = 1/2*log((1+x)/(1-x)).
239 - improve mpfr_gamma (see http://code.google.com/p/fastfunlib/). A possible
240 idea is to implement a fast algorithm for the argument reconstruction
241 gamma(x+k). One could also use the series for 1/gamma(x), see for example
242 http://dlmf.nist.gov/5/7/ or formula (36) from
243 http://mathworld.wolfram.com/GammaFunction.html
244 - fix regression with mpfr_mpz_root (from Keith Briggs, 5 July 2006), for
245 example on 3Ghz P4 with gmp-4.2, x=12.345:
246 prec=50000 k=2 k=3 k=10 k=100
247 mpz_root 0.036 0.072 0.476 7.628
248 mpfr_mpz_root 0.004 0.004 0.036 12.20
249 See also mail from Carl Witty on mpfr list, 09 Oct 2007.
250 - implement Mulders algorithm for squaring and division
251 - for sparse input (say x=1 with 2 bits), mpfr_exp is not faster than for
252 full precision when precision <= MPFR_EXP_THRESHOLD. The reason is
253 that argument reduction kills sparsity. Maybe avoid argument reduction
255 - speed up const_euler for large precision [for x=1.1, prec=16610, it takes
256 75% of the total time of eint(x)!]
257 - speed up mpfr_atan for large arguments (to speed up mpc_log)
258 [from Mark Watkins on Fri, 18 Mar 2005]
259 Also mpfr_atan(x) seems slower (by a factor of 2) for x near from 1.
260 Example on a Athlon for 10^5 bits: x=1.1 takes 3s, whereas 2.1 takes 1.8s.
261 The current implementation does not give monotonous timing for the following:
262 mpfr_random (x); for (i = 0; i < k; i++) mpfr_atan (y, x, MPFR_RNDN);
263 for precision 300 and k=1000, we get 1070ms, and 500ms only for p=400!
264 - improve mpfr_sin on values like ~pi (do not compute sin from cos, because
265 of the cancellation). For instance, reduce the input modulo pi/2 in
266 [-pi/4,pi/4], and define auxiliary functions for which the argument is
267 assumed to be already reduced (so that the sin function can avoid
268 unnecessary computations by calling the auxiliary cos function instead of
269 the full cos function). This will require a native code for sin, for
270 example using the reduction sin(3x)=3sin(x)-4sin(x)^3.
271 See http://websympa.loria.fr/wwsympa/arc/mpfr/2007-08/msg00001.html and
272 the following messages.
273 - improve generic.c to work for number of terms <> 2^k
274 - rewrite mpfr_greater_p... as native code.
275 - inline mpfr_neg? Problems with NAN flags:
276 #define mpfr_neg(_d,_x,_r) \
277 (__builtin_constant_p ((_d)==(_x)) && (_d)==(_x) ? \
278 ((_d)->_mpfr_sign = -(_d)->_mpfr_sign, 0) : \
279 mpfr_neg ((_d), (_x), (_r))) */
281 - mpf_t uses a scheme where the number of limbs actually present can
282 be less than the selected precision, thereby allowing low precision
283 values (for instance small integers) to be stored and manipulated in
284 an mpf_t efficiently.
286 Perhaps mpfr should get something similar, especially if looking to
287 replace mpf with mpfr, though it'd be a major change. Alternately
288 perhaps those mpfr routines like mpfr_mul where optimizations are
289 possible through stripping low zero bits or limbs could check for
290 that (this would be less efficient but easier).
292 - try the idea of the paper "Reduced Cancellation in the Evaluation of Entire
293 Functions and Applications to the Error Function" by W. Gawronski, J. Mueller
294 and M. Reinhard, to be published in SIAM Journal on Numerical Analysis: to
295 avoid cancellation in say erfc(x) for x large, they compute the Taylor
296 expansion of erfc(x)*exp(x^2/2) instead (which has less cancellation),
297 and then divide by exp(x^2/2) (which is simpler to compute).
299 - replace the *_THRESHOLD macros by global (TLS) variables that can be
300 changed at run time (via a function, like other variables)? One benefit
301 is that users could use a single MPFR binary on several machines (e.g.,
302 a library provided by binary packages or shared via NFS) with different
303 thresholds. On the default values, this would be a bit less efficient
304 than the current code, but this isn't probably noticeable (this should
305 be tested). Something like:
306 long *mpfr_tune_get(void) to get the current values (the first value
307 is the size of the array).
308 int mpfr_tune_set(long *array) to set the tune values.
309 int mpfr_tune_run(long level) to find the best values (the support
310 for this feature is optional, this can also be done with an
313 - better distinguish different processors (for example Opteron and Core 2)
314 and use corresponding default tuning parameters (as in GMP). This could be
315 done in configure.in to avoid hacking config.guess, for example define
317 Note (VL): the effect on cross-compilation (that can be a processor
318 with the same architecture, e.g. compilation on a Core 2 for an
319 Opteron) is not clear. The choice should be consistent with the
320 build target (e.g. -march or -mtune value with gcc).
321 Also choose better default values. For instance, the default value of
322 MPFR_MUL_THRESHOLD is 40, while the best values that have been found
323 are between 11 and 19 for 32 bits and between 4 and 10 for 64 bits!
325 - during the Many Digits competition, we noticed that (our implantation of)
326 Mulders short product was slower than a full product for large sizes.
327 This should be precisely analyzed and fixed if needed.
329 ##############################################################################
331 ##############################################################################
333 - Once the double inclusion of mpfr.h is fully supported, add tstdint
334 to check_PROGRAMS in the tests/Makefile.am file.
336 - [suggested by Tobias Burnus <burnus(at)net-b.de> and
337 Asher Langton <langton(at)gcc.gnu.org>, Wed, 01 Aug 2007]
338 support quiet and signaling NaNs in mpfr:
339 * functions to set/test a quiet/signaling NaN: mpfr_set_snan, mpfr_snan_p,
340 mpfr_set_qnan, mpfr_qnan_p
341 * correctly convert to/from double (if encoding of s/qNaN is fixed in 754R)
343 - check again coverage: on July 27, Patrick Pelissier reports that the
344 following files are not tested at 100%: add1.c, atan.c, atan2.c,
345 cache.c, cmp2.c, const_catalan.c, const_euler.c, const_log2.c, cos.c,
346 gen_inverse.h, div_ui.c, eint.c, exp3.c, exp_2.c, expm1.c, fma.c, fms.c,
347 lngamma.c, gamma.c, get_d.c, get_f.c, get_ld.c, get_str.c, get_z.c,
348 inp_str.c, jn.c, jyn_asympt.c, lngamma.c, mpfr-gmp.c, mul.c, mul_ui.c,
349 mulders.c, out_str.c, pow.c, print_raw.c, rint.c, root.c, round_near_x.c,
350 round_raw_generic.c, set_d.c, set_ld.c, set_q.c, set_uj.c, set_z.c, sin.c,
351 sin_cos.c, sinh.c, sqr.c, stack_interface.c, sub1.c, sub1sp.c, subnormal.c,
352 uceil_exp2.c, uceil_log2.c, ui_pow_ui.c, urandomb.c, yn.c, zeta.c, zeta_ui.c.
354 - check the constants mpfr_set_emin (-16382-63) and mpfr_set_emax (16383) in
355 get_ld.c and the other constants, and provide a testcase for large and
358 - from Kevin Ryde <user42@zip.com.au>:
359 Also for pi.c, a pre-calculated compiled-in pi to a few thousand
360 digits would be good value I think. After all, say 10000 bits using
361 1250 bytes would still be small compared to the code size!
362 Store pi in round to zero mode (to recover other modes).
364 - add a new rounding mode: round to nearest, with ties away from zero
365 (this is roundTiesToAway in 754-2008, could be used by mpfr_round)
366 - add a new roundind mode: round to odd. If the result is not exactly
367 representable, then round to the odd mantissa. This rounding
368 has the nice property that for k > 1, if:
369 y = round(x, p+k, TO_ODD)
370 z = round(y, p, TO_NEAREST_EVEN), then
371 z = round(x, p, TO_NEAREST_EVEN)
372 so it avoids the double-rounding problem.
374 - add tests of the ternary value for constants
376 - When doing Extensive Check (--enable-assert=full), since all the
377 functions use a similar use of MACROS (ZivLoop, ROUND_P), it should
378 be possible to do such a scheme:
379 For the first call to ROUND_P when we can round.
380 Mark it as such and save the approximated rounding value in
381 a temporary variable.
382 Then after, if the mark is set, check if:
383 - we still can round.
384 - The rounded value is the same.
385 It should be a complement to tgeneric tests.
387 - add a new exception "division by zero" (IEEE-754 terminology) / "infinitary"
388 (LIA-2 terminology). In IEEE 754R (2006 February 14 8:00):
389 "The division by zero exception shall be signaled iff an exact
390 infinite result is defined for an operation on finite operands.
391 [such as a pole or logarithmic singularity.] In particular, the
392 division by zero exception shall be signaled if the divisor is
393 zero and the dividend is a finite nonzero number."
395 - in div.c, try to find a case for which cy != 0 after the line
396 cy = mpn_sub_1 (sp + k, sp + k, qsize, cy);
397 (which should be added to the tests), e.g. by having {vp, k} = 0, or
398 prove that this cannot happen.
400 - add a configure test for --enable-logging to ignore the option if
401 it cannot be supported. Modify the "configure --help" description
402 to say "on systems that support it".
404 - allow generic tests to run with a restricted exponent range.
406 - add generic bad cases for functions that don't have an inverse
407 function that is implemented (use a single Newton iteration).
409 - add bad cases for the internal error bound (by using a dichotomy
410 between a bad case for the correct rounding and some input value
411 with fewer Ziv iterations?).
413 - add an option to use a 32-bit exponent type (int) on LP64 machines,
414 mainly for developers, in order to be able to test the case where the
415 extended exponent range is the same as the default exponent range, on
418 - test underflow/overflow detection of various functions (in particular
419 mpfr_exp) in reduced exponent ranges, including ranges that do not
423 ##############################################################################
425 ##############################################################################
427 - support the decimal64 function without requiring --with-gmp-build
429 - [Kevin about texp.c long strings]
430 For strings longer than c99 guarantees, it might be cleaner to
431 introduce a "tests_strdupcat" or something to concatenate literal
432 strings into newly allocated memory. I thought I'd done that in a
433 couple of places already. Arrays of chars are not much fun.
435 - use http://gcc.gnu.org/viewcvs/trunk/config/stdint.m4 for mpfr-gmp.h
437 - rename configure.in to configure.ac