(Free42) accuracy of LN1+X

[Werner]

I would expect LN1+X(a) to return *exactly* a for a < 1e-34, but that is not the case.
eg. for a=1e-36, LN1+X returns 9.999...e-37, likewise for a=1e-37.
a=1e-38 is correct etc. I gather LN1+X is not a native function of the Intel library?

Cheers, Werner

[Albert Chan]

Perhaps it is implemented with binary128 version log1p ?
Example:

>>> from gmpy2 import *
>>> mpfr('1e-99',113)       # binary128 precision = 113 bits
mpfr('9.99999999999999999999999999999999925e-100',113)

Free42:

1e-99 LN1+X     → 9.999999999999999999999999999999999e-100

[grsbanks]

(06-12-2020 06:51 PM)Albert Chan Wrote:  Perhaps it is implemented with binary128 version log1p ?

Looks like it is.

Code:

Phloat log1p(Phloat p) {
    BID_UINT128 res;
    bid128_log1p(&res, &p.val);
    return Phloat(res);
}

[johanw]

(06-12-2020 03:23 PM)Werner Wrote:  I would expect LN1+X(a) to return *exactly* a for a < 1e-34, but that is not the case.
eg. for a=1e-36, LN1+X returns 9.999...e-37, likewise for a=1e-37.
a=1e-38 is correct etc. I gather LN1+X is not a native function of the Intel library?

I have been trying to understand this post but I failed. I probably overlook something trivial, but what is X(a)?

When I enter LN(1 + 1e-34) in Free42 I get 0 as an answer (1e-34 enter 1 + LN).

[rprosperi]

The function Werner is referring to is "LN1+X( )", so LN1+X(A) means evaluating that function for the value A.

This is a special version of the LN function for argument (meaning X) values very close to 0, which otherwise would not calculate as accurately.

[Albert Chan]

With such a tiny x, log1p(x) is just going to return x

Error is due to decimal128 unable to round-trip thru binary128
You need ceil(34 * log2(10)) + 1 = 114 bits, see Number of Digits required for round-trip conversions

Trivia: LN1+X(ε) will not return (1.000 ... 1) ε

Assuming you do get a result of (1.000 ... 1) ε
This implied decimal128 to binary128 returns ≥ (10^33 + ½) / (ε/10^33)
(if ε = 1/10^n, replace "≥" with ">", due to round-to-even rule)

That meant it required relative error ≥ ½ / 10^33 = 500e-36

binary128 (113 bits precision), relative error ≤ (½ ulp) / (2^112 ulp) ≈ 96e-36

Assumption were wrong, you will not see this. QED

However, returning (0.999 ...) ε is possible.
Required relative error dropped by a factor of 10, to 50e-36 < 96e-36

Update:
(0.999 ... 8) ε required relative error ≥ 3 * 50e-36 = 150e-36, thus not possible.
Based on relative errors, only 2 patterns possible, (1.000 ...) ε, or (0.999 ...) ε

[ijabbott]

(06-12-2020 08:21 PM)grsbanks Wrote:  

(06-12-2020 06:51 PM)Albert Chan Wrote:  Perhaps it is implemented with binary128 version log1p ?

Looks like it is.

Code:

Phloat log1p(Phloat p) {
    BID_UINT128 res;
    bid128_log1p(&res, &p.val);
    return Phloat(res);
}

That code snippet by itself says nothing about the implementation of bid128_log1p, although looking at the source for the function itself (e.g. this copy on GitHub), most cases are wrapped around calls to bid128_to_binary128 and binary128_to_bid128 with the real work done in binary.

[Paul Dale]

I've not checked the Intel code for this function, but the library does make use of binary floating point functions in several places.

This function is actually fairly easy to implement once its problems are realised. The code on the 34S is:

Code:

decNumber *decNumberLn1p(decNumber *r, const decNumber *x) {
    decNumber u, v, w;

    if (decNumberIsSpecial(x) || dn_eq0(x)) {
        return decNumberCopy(r, x);
    }
    dn_p1(&u, x);
    dn_m1(&v, &u);
    if (dn_eq0(&v)) {
        return decNumberCopy(r, x);
    }
    dn_divide(&w, x, &v);
    dn_ln(&v, &u);
    return dn_multiply(r, &v, &w);
}

Pauli

[Albert Chan]

(06-14-2020 01:09 AM)Paul Dale Wrote:  This function is actually fairly easy to implement once its problems are realised.

Probably not.

decimal128 unable to round-trip thru binary128 occurs at all ranges, not just tiny x
We just don't noticed it when x is bigger ...

For example, try this:

1e-12 XEQ "LN1+X"     → 9.999999999995000000000003333333334e-13

We expected result should be rounded(x - x²/2 + x³/3), with last digit = 3

x = 1e-12 converted to binary128, we got x' ≈ (10^33+0.08)/10^45, over-estimated 0.45 ULP

log1p(x') ≈ 9.99999999999500000000000333333333413e-13

This example is lucky that roundings compensated the x'-x error.
You can potientially get into double-rounding errors, all on the same side.

However, even if conversion is exact, we still cannot get correctly rounded result.

>>> from gmpy2 import *
>>> get_context().precision = 113 # binary128
>>> y = log1p(mpfr('1e-12'))
>>> format(y, '.33e')
'9.999999999995000000000003333333334e-13'
>>> format(next_below(y), '.33e')
'9.999999999995000000000003333333332e-13'

[Paul Dale]

I wouldn't say using a round trip decimal to binary and back is a sensible way to implement this exactly.... Fast yes, accurate not really.