Post Reply 
(Free42) accuracy of LN1+X
06-12-2020, 03:23 PM
Post: #1
(Free42) accuracy of LN1+X
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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
Find all posts by this user
Quote this message in a reply
06-12-2020, 06:51 PM
Post: #2
RE: (Free42) accuracy of LN1+X
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
Find all posts by this user
Quote this message in a reply
06-12-2020, 08:21 PM
Post: #3
RE: (Free42) accuracy of LN1+X
(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);
}

There are only 10 types of people in this world. Those who understand binary and those who don't.
Find all posts by this user
Quote this message in a reply
06-12-2020, 10:07 PM
Post: #4
RE: (Free42) accuracy of LN1+X
(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).
Find all posts by this user
Quote this message in a reply
06-12-2020, 10:45 PM (This post was last modified: 06-12-2020 11:24 PM by rprosperi.)
Post: #5
RE: (Free42) accuracy of LN1+X
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.

--Bob Prosperi
Find all posts by this user
Quote this message in a reply
06-13-2020, 12:45 AM (This post was last modified: 06-13-2020 10:46 AM by Albert Chan.)
Post: #6
RE: (Free42) accuracy of LN1+X
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 ...) ε
Find all posts by this user
Quote this message in a reply
06-13-2020, 11:07 AM
Post: #7
RE: (Free42) accuracy of LN1+X
(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.

— Ian Abbott
Find all posts by this user
Quote this message in a reply
06-14-2020, 01:09 AM
Post: #8
RE: (Free42) accuracy of LN1+X
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
Find all posts by this user
Quote this message in a reply
06-14-2020, 12:48 PM
Post: #9
RE: (Free42) accuracy of LN1+X
(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'
Find all posts by this user
Quote this message in a reply
06-14-2020, 10:39 PM
Post: #10
RE: (Free42) accuracy of LN1+X
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.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)