Re: Use of C99 extra long double math functions after r236148

On Sun, Jul 08, 2012 at 02:06:46PM -0500, Stephen Montgomery-Smith wrote:
> Here is a technical question.  I know that people always talk about 
> ulp's in the context of how good a function implementation is.  I think 
> the ulp is the number of base 2 digits at the end of the mantissa that 
> we cannot rely on.

There are a few definitions of ULP (units in the last place).
Google 'Muller ULP' and check goldberg's paper (google
'goldberg floating point'.

Here's how I define ULP for my testing and the MPFR code.

/* From a das email:
     ulps = err * (2**(LDBL_MANT_DIG - e)), where e = ilogb(z).

 I use:

     mpfr_sub(err, acc_out, approx_out, GMP_RNDN);
     mpfr_abs(err, err, GMP_RNDN);
     mpfr_mul_2si(ulpx, err, -mpfr_get_exp(acc_out) + LDBL_MANT_DIG, GMP_RNDN);
     ulp = mpfr_get_d(ulpx, GMP_RNDN);
     if (ulp  0.506 && mpfr_cmpabs(err, ldbl_minx)  0) {
       print warning...;
     }
*/
#include "gmp.h"
#include "mpfr.h"
#include "sgk.h"

static mpfr_t mp_ulp_err;

/*
 * Set the precision of the ULP computation.
 */
void
mp_ulp_init(int prec)
{
	mpfr_init2(mp_ulp_err, (mpfr_prec_t)prec);
}

/*
 * Given the 'approximate' value of a function and
 * an 'accurate' value compute the ULP.
 */
double
mp_ulp(mpfr_t approximate, mpfr_t accurate, int dig)
{
        double ulp;
        mpfr_sub(mp_ulp_err, accurate, approximate, RD);
        mpfr_abs(mp_ulp_err, mp_ulp_err, RD);
        mpfr_mul_2si(mp_ulp_err, mp_ulp_err, - mpfr_get_exp(accurate) + dig,
	    RD);
        ulp = mpfr_get_d(mp_ulp_err, RD);
        return (ulp);   
}

> So if one were to write a naive implementation of lexp(x) that used 
> Taylor's series if x is positive, and 1/lexp(-x) is x is negative - well 
> one could fairly easily estimate an upper bound on the ulp, but it 
> wouldn't be low (like ulp=1 or 2), but probably rather higher (ulp of 
> the order of 10 or 20).

I've already written an ld80 expl().  I have tested billions if not
trillions of values.  Don't recall the max ULP I saw, but I know it
was less than 1.  I don't have an ld128 version, so I won't submit
a patch for libm.

> So do people really work hard to get that last drop of ulp out of their 
> calculations?

I know very few scientist who work hard to reduce the ULP.  Most
have little understanding of floating point.  

>  Would a ulp=10 be considered unacceptable?

Yes, it is unacceptable for the math library.  Remember ULPs
propagate and can quickly grow if the initial ulp for a 
result is large. 

> Also, looking through the source code for the FreeBSD implementation of 
> exp, I saw that they used some rather smart rational function.  (I don't 
> know how they came up with it.)

See Steven Moshier's book, which is the basis for CEPHES.
Instead of using a minimax polynomial approximation for
the function in some interval, one uses a minimax approximation
with a rational function.

-- 
Steve