+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: Not HP Calculators (/forum-7.html)
+--- Forum: Not quite HP Calculators - but related (/forum-8.html)
+--- Thread: (Free42) accuracy of LN1+X (/thread-15180.html)
(Free42) accuracy of LN1+X - Werner - 06-12-2020 03:23 PM
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
RE: (Free42) accuracy of LN1+X - Albert Chan - 06-12-2020 06:51 PM
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
RE: (Free42) accuracy of LN1+X - grsbanks - 06-12-2020 08:21 PM
(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);
}RE: (Free42) accuracy of LN1+X - johanw - 06-12-2020 10:07 PM
(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.I have been trying to understand this post but I failed. I probably overlook something trivial, but what is X(a)?
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?
When I enter LN(1 + 1e-34) in Free42 I get 0 as an answer (1e-34 enter 1 + LN).
RE: (Free42) accuracy of LN1+X - rprosperi - 06-12-2020 10:45 PM
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.
RE: (Free42) accuracy of LN1+X - Albert Chan - 06-13-2020 12:45 AM
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 ...) ε
RE: (Free42) accuracy of LN1+X - ijabbott - 06-13-2020 11:07 AM
(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.
RE: (Free42) accuracy of LN1+X - Paul Dale - 06-14-2020 01:09 AM
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
RE: (Free42) accuracy of LN1+X - Albert Chan - 06-14-2020 12:48 PM
(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'
RE: (Free42) accuracy of LN1+X - Paul Dale - 06-14-2020 10:39 PM
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.