Re: What to do about tgammal?

On Sat, Dec 18, 2021 at 10:41:14AM +0000, Mark Murray wrote:
> 
> Hmm. I think my understanding of ULP is missing something?
> 
> I thought that ULP could not be greater than the mantissa size
> in bits?
> 
> I.e., I thought it represents average rounding error (compared with
> "perfect rounding"), not truncation error, as the above very large
> ULPs suggest.
> 

The definition of ULP differs according which expert you
choose to follow. :-)  For me (a non-expert), ULP is measured
in the system of the "accurate answer", which is assumed to
have many more bits of precision than the "approximate answer".
From a very old das@ email and for long double I have 

/* From a das email:
 *
 *   ulps = err * (2**(LDBL_MANT_DIG - e)), where e = ilogb(accurate).
 *
 * I use:
 *
 *   mpfr_sub(err, accurate, approximate, RNDN);
 *   mpfr_abs(err, err, RNDN);
 *   mpfr_mul_2si(ulpx, err, -mpfr_get_exp(accurate) + LDBL_MANT_DIG, RNDN);
 *   ulp = mpfr_get_d(ulpx, RNDN);
 *   if (ulp  0.506 && mpfr_cmpabs(err, ldbl_minx)  0) {
 *     print warning...;
 *   }
 */

Here, "err = |accurate - approximate|" is done in the precision
of the "accurate answer", and RNDN is round-to-nearest.  The
line 'mpfr_mul_2si(...)' is doing the scaling to manipulate the
exponent of the radix 2.

This is agreement with Goldberg.  IIRC, Jean-Michel Muller et al
have much more detailed discussion about error and ULPs in their
book "Handbook of Floating-Point Arithmetic".   

https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

So, on ld80 and tgamma, I have 

% tlibm -l -x 8 -X 1755 -n 100000 -P tgamma
Interval tested for tgammal: [8,1755]
       ulp <= 0.5: 90.207%     90207 |  90.207%     90207
0.5 <  ulp <  0.6:  7.085%      7085 |  97.292%     97292
0.6 <  ulp <  0.7:  2.383%      2383 |  99.675%     99675
0.7 <  ulp <  0.8:  0.314%       314 |  99.989%     99989
0.8 <  ulp <  0.9:  0.011%        11 | 100.000%    100000
Max ulp: 0.853190 at 9.236118561185611856579e+02

The "ulp <= 0.5" line indicates the number of correctly
rounded case.  If the value of ULP exceeds 1, then you are
starting to lose more than 1 bit in the significant.  Consider
individual points.  Here's a correctly rounded case:

% tlibm -l -a 9.236118561185000000 tgamma
   xu = LD80C(0x93c72441a44a57a9,     3,  9.23611856118499999994e+00L),
libmu = LD80C(0x82f85b3fe4b36f9b,    16,  6.70567128873706962722e+04L),
mpfru = LD80C(0x82f85b3fe4b36f9b,    16,  6.70567128873706962722e+04L),
 ULP = 0.42318

Notice, ULP < 0.5 and the bits in libm and mpfr answer rounded to 
long double are all the same as is the decimal representation. Now
consider the max ULP case:

% tlibm -l -a  9.236118561185611856579e+02 tgamma
   xu = LD80C(0xe6e728a690c4fa87,     9,  9.23611856118561185658e+02L),
libmu = LD80C(0xba6bea661a7eda3c,  7762,  5.72944398777756327202e+2336L),
mpfru = LD80C(0xba6bea661a7eda3d,  7762,  5.72944398777756327245e+2336L),
 ULP = 0.85319

ULP is < 1, which is desireable.  The bits agree except for the
last place where we're off by 1, ie, 0xc vs 0xd.  The decimal
representation is of course off.

So, now on grimoirei without my patch, I have 

% ./tlibm_libm -l -a 1.59255925592559255925592559255925594e+01 tgamma
   x =  1.59255925592559255925592559255925594e+01L
libm =  1.06660136733166064453125000000000000e+12L
mpfr =  1.06660136733166211967839038834033459e+12L,
 ULP = 13932370826141802496.00000

You have 16 correct decimal digits out of 36, which is the best
you can expect from mapping tgammal to tgamma.  ULP is significant
larger than 1.  The ulp is almost guaranteed to be larger than
2**(113-53).

Now with my patch and worse case for 10000 x in [10:16]

% ./tlibm_lmath -l -a 1.59255925592559255925592559255925594e+01 tgamma
  x =  1.59255925592559255925592559255925594e+01L
libm =  1.06660136733166211967839038834033897e+12
mpfr =  1.06660136733166211967839038834033459e+12L,
 ULP = 41.40877

I don't print out the hex representation in ld128, but you see
the number of correct decimal digits is 33 digits compared to
36.

-- 
Steve