path: root/TODO

Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 Free Software Foundation, Inc.
Contributed by the Spaces project, INRIA Lorraine.

This file is part of the MPFR Library.

The MPFR Library is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License (either version 2.1
of the License, or, at your option, any later version) and the GNU General
Public License as published by the Free Software Foundation (most of MPFR is
under the former, some under the latter).

The 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 MPFR Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
02110-1301, USA.


##############################################################################
Documentation:
##############################################################################

- add a description of the algorithms used + proof of correctness

- mpfr_set_prec: add an explanation of how to speed up calculations
  which increase their precision at each step.


##############################################################################
Installation:
##############################################################################

- nothing to do currently :-)


##############################################################################
Changes in existing functions:
##############################################################################

- many functions currently taking into account the precision of the *input*
  variable to set the initial working precison (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_get_str should support base up to 62 too.

- 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 convertion 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 exaclty n bits.
  If RNDN, takes nq+1 bits. (See also the new division function).

- random functions: get rid of _gmp_rand.


##############################################################################
New functions to implement:
##############################################################################

- document MPFR_SIGN and/or provide mpfr_get_sign / mpfr_set_sign.
- function mpfr_lgamma, with the prototype
    int mpfr_lgamma (mpfr_ptr, int *, mpfr_srcptr, mpfr_rnd_t)
  that also returns the value of signgam; see thread
    http://sympa.loria.fr/wwsympa/arc/mpfr/2006-10/msg00029.html
- functions operating on mpfr_t and double: mpfr_add_d, mpfr_sub_d,
  mpfr_d_sub, mpfr_mul_d, mpfr_div_d, mpfr_d_div
  [suggested by Keith Briggs, 3 Jan 2006]
- modf (to extract integer and fractional parts), suggested by
        Dmitry Antipov <dmitry.antipov@mail.ru> Thu, 13 Jun 2002
- mpfr_fmod (mpfr_t, mpfr_srcptr, mpfr_srcptr, mp_rnd_t)
  [suggested by Tomas Zahradnicky <tomas@24uSoftware.com>, 29 Nov 2003]
  Kevin: Might want to be called mpfr_mod, to match mpz_mod.
  -> we probably want to allow both mpfr_fmod and mpfr_mod.
  Proposed implementation (apart from special cases):
  int mpfr_fmod (mpfr_t r, mpfr_t x, mpfr_t y, mpfr_rnd_t r)
  {
     mpfr_t q;
     int inexact;

     /* if 2^(ex-1) <= |x| < 2^ex, and 2^(ey-1) <= |y| < 2^ey, then
	|x/y| < 2^(ex-ey+1) */
     mpfr_init2 (q, MAX(MPFR_PREC_MIN, mpfr_get_exp (x)-mpfr_get_exp (y) + 1));
     mpfr_div (q, x, y, GMP_RNDZ);
     mpfr_trunc (q, q);
     mpfr_prec_round (q, mpfr_get_prec (q) + mpfr_get_prec (y), GMP_RNDZ);
     mpfr_mul (q, q, y, GMP_RNDZ); /* exact */
     inexact = mpfr_sub (r, x, q, r);
     mpfr_clear (q);
     return inexact;
  }
- mpfr_fms for a-b*c
  [suggested by Tomas Zahradnicky <tomas@24uSoftware.com>, 29 Nov 2003]
- 1/sqrt(x) [Regis Dupont <dupont@lix.polytechnique.fr>, 15 Sep 2004]
- dilog() [the dilogarithm: dilog(x) =  int(ln(t)/(1-t), t=1..x)]
- mpfr_printf [Sisyphus <kalinabears@iinet.net.au> Tue, 04 Jan 2005]
  for example mpfr_printf ("%.2Ff\n", x)
- 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)
  IncompleteGamma
  KBessel, KBessel2 [2nd kind]
  JacobiTheta
  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
  DawsonIntegral
  Psi = LogDerivative
  GammaD(x) = Gamma(x+1/2)
- 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;
  + rec_sqrt (5.2.6): 1 / sqrt(x);
  + 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.;
  + 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 http://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 http://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)
    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))
  + Trigamma Function psi'(x).
    and Polygamma Function: psi^{(m)}(x) for m >= 0, x > 0.

- from gnumeric (www.gnome.org/projects/gnumeric/doc/function-reference.html):
  - beta
  - betaln
  - degrees
  - radians
  - sqrtpi

- mpfr_frexp(mpfr_t rop, mp_exp_t *n, mpfr_t op, mp_rnd_t rnd) suggested by
  Steve Kargl <sgk@troutmask.apl.washington.edu> Sun, 7 Aug 2005


##############################################################################
Efficiency:
##############################################################################

- 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
- 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, GMP_RNDN);
  for precision 300 and k=1000, we get 1070ms, and 500ms only for p=400!
- implement range reduction in sin/cos/tan for large arguments
        (currently too slow for 2^1024), and/or implement Mark Watkins's
        algorithm (see also remquo):
  My suggestion for a long-term fix to these sin/cos problems when one
  of them is close to zero is the following:
  Give input inp:
  1) Reduce inp to the interval 0 to Pi/2
  2) If inp is less than Pi/4, compute sin(inp) via a power series around 0.
  3) If inp is more than Pi/4, compute s=Pi/2-inp to the working precision
     (most of the possible cancellation will occur here --- or in step 1)
     and calculate cos(inp) as a power series (in s) around Pi/2 --- thus
     the output is near zero precisely when the s-input is near zero,
     avoiding most cancellation problems.
  4) Having sin or cos, compute the other via the sqrt. This will not cause
     cancellation problems because the computed sin or cos is not near 1.
- combined mpfr_sinh_cosh() [Geoff Bailey, 20 Apr 2005,
  and Kaveh R. Ghazi, 17 Jan 2007]
- improve generic.c to work for number of terms <> 2^k
- rewrite mpfr_greater_p... as native code.
- inline mpfr_neg? Problems with NAN flags:
   #define mpfr_neg(_d,_x,_r)                          \
    (__builtin_constant_p ((_d)==(_x)) && (_d)==(_x) ? \
     ((_d)->_mpfr_sign = -(_d)->_mpfr_sign, 0)       : \
      mpfr_neg ((_d), (_x), (_r)))  */

- 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).


##############################################################################
Miscellaneous:
##############################################################################

- In addition to major/minor/patchlevel version information, add a macro
  and/or function to the MPFR interface, that would tell what patches
  have been applied and other possible information (some information is
  already available as a suffix in the MPFR_VERSION_STRING macro, but it
  should be short, thus it is not suitable if one wishes to mention all
  the individual patches that have been applied); this can include the
  features that have been enabled or disabled in configure. For instance,
  see what Mutt and Subversion do concerning their CVS/Subversion builds.
  Of course, patches must include some tag (one should not rely on a
  wrapper to apply them), that could be added to a PATCHES file. Mutt's
  doc/patch-notes.txt file mentions an interesting scheme: the following
  chunk should be included in the patch, after replacing <your-id-here>
  by the patch id.
  ------------------------------snip------------------------------
  --- PATCHES~    Tue Nov  6 19:59:33 2001
  +++ PATCHES     Tue Nov  6 19:59:42 2001
  @@ -1,0 +1 @@
  +<your-id-here>
  ------------------------------snip------------------------------
  References:
    http://sympa.loria.fr/wwsympa/arc/mpfr/2005-12/msg00049.html
    http://gcc.gnu.org/bugzilla/show_bug.cgi?id=29405
    http://svn.collab.net/repos/svn/trunk/subversion/include/svn_version.h

- 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.

- rename mpf2mpfr.h to gmp-mpf2mpfr.h?
  (will wait until mpfr is fully integrated into gmp :-)

- 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 a new rounding mode: rounding away from 0. This can be easily
  implemented as follows: round to zero, and if the result is inexact,
        add one ulp to the mantissa.
- add a new rounding mode: round to nearest, with ties away from zero
  (will be in 754r, could be used by mpfr_round)
- 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.

- check/define the sign of infinity for gamma(-integer)

- 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.

- add a new exception "division by zero" (IEEE-754 terminology) / "pole"
  (LIA-2 terminology). In IEEE 754R (2006 February 14 8:00):
    "The division by zero exception shall be signaled iff an exact
    infinite result is defined for an operation on finite operands.
    [such as a pole or logarithmic singularity.] In particular, the
    division by zero exception shall be signaled if the divisor is
    zero and the dividend is a finite nonzero number."

- 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".


##############################################################################
Portability:
##############################################################################

- [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 http://gcc.gnu.org/viewcvs/trunk/config/stdint.m4 for mpfr-gmp.h


##############################################################################
Possible future MPF / MPFR integration:
##############################################################################

- mpf routines can become "extern inline"s calling mpfr equivalents,
  probably just with GMP_RNDZ hard coded, since that's what mpf has
  always done.

- Want to preserve the mpf_t structure size, for binary compatibility.
  Breaking compatibility would cause lots of pain and potential subtle
  breakage for users.  If the fields in mpf_t are not enough then
  extra space under _mp_d can be used.

- mpf_sgn has been a macro directly accessing the _mp_size field, so a
  compatible representation would be required.  At worst that field
  could be maintained for mpf_sgn, but not otherwise used internally.

  mpf_sgn should probably throw an exception if called with NaN, since
  there's no useful value it can return, so it might want to become a
  function.  Inlined copies in existing binaries would hopefully never
  see a NaN, if they only do old-style mpf things.

- mpfr routines replacing mpf routines must be reentrant and thread
  safe, since of course that's what has been documented for mpf.

- mpfr_random will not be wanted since there's no corresponding
  mpf_random and new routines should not use the old style global
  random state.