Re: BUG #16281: LN() function inaccurate at 1000th fractional digit

On Fri, 28 Feb 2020 at 15:39, Tom Lane <tgl(at)sss(dot)pgh(dot)pa(dot)us> wrote:
>
> Dean Rasheed <dean(dot)a(dot)rasheed(at)gmail(dot)com> writes:
>
> > Yeah, that's not at all surprising. All the transcendental numeric
> > functions are basically +/-1 in the final digit.
>
> TBH, I'd be ecstatic if I thought that was the maximum possible error ;-).
> I don't think anyone's really done any error propagation analysis on
> the numeric transcendentals.
>

There's a good paper [1] on that subject by the authors of MPFR. A
prerequisite for keeping errors under control is having appropriately
chosen rscales throughout, so that errors grow in a controlled way
with each operation, which [2] went a long way towards ensuring.
Barring oversights in [2], I think it's entirely plausible that the we
do have that kind of accuracy, given the fairly generous local rscales
chosen in these functions, compared with the number of internal
operations they perform (e.g., I think that ln_var() never needs more
than around 400 terms in its Taylor series), but it would take a much
more thorough analysis to prove it. Speaking of oversights though...

Looking more closely at ln_var(), it seems that there was an oversight
related to the way that it computes the local_rscale for the Taylor
series expansion --- it fails to account for the fact that the result
is multiplied by fact (2^(nsqrt+1), where nsqrt is the number of
square roots performed in the range reduction phase, which in practice
is at most 22).

Since 2^22 has 7 decimal digits, multiplying by 2^22 almost entirely
wipes out the 8-digit safety margin used in the Taylor series
expansion. The attached patch corrects that.

With this patch, all the examples originally posted return the correct
results (calculated with bc). I'd be interested to know how the OP
constructed these examples, and whether there were any that were off
by more than 1 ULP.

Regards,
Dean

[1] https://www.mpfr.org/algorithms.pdf
[2] https://www.postgresql.org/message-id/E1Zxgva-0000tn-Ox%40gemulon.postgresql.org