(Free42) accuracy of LN1+X

06122020, 03:23 PM
Post: #1




(Free42) accuracy of LN1+X
I would expect LN1+X(a) to return *exactly* a for a < 1e34, but that is not the case.
eg. for a=1e36, LN1+X returns 9.999...e37, likewise for a=1e37. a=1e38 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 

06122020, 06:51 PM
Post: #2




RE: (Free42) accuracy of LN1+X
Perhaps it is implemented with binary128 version log1p ?
Example: >>> from gmpy2 import * >>> mpfr('1e99',113) # binary128 precision = 113 bits mpfr('9.99999999999999999999999999999999925e100',113) Free42: 1e99 LN1+X → 9.999999999999999999999999999999999e100 

06122020, 08:21 PM
Post: #3




RE: (Free42) accuracy of LN1+X
(06122020 06:51 PM)Albert Chan Wrote: Perhaps it is implemented with binary128 version log1p ? Looks like it is. Code: Phloat log1p(Phloat p) { There are only 10 types of people in this world. Those who understand binary and those who don't. 

06122020, 10:07 PM
Post: #4




RE: (Free42) accuracy of LN1+X
(06122020 03:23 PM)Werner Wrote: I would expect LN1+X(a) to return *exactly* a for a < 1e34, 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)? When I enter LN(1 + 1e34) in Free42 I get 0 as an answer (1e34 enter 1 + LN). 

06122020, 10:45 PM
(This post was last modified: 06122020 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 

06132020, 12:45 AM
(This post was last modified: 06132020 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 roundtrip thru binary128 You need ceil(34 * log2(10)) + 1 = 114 bits, see Number of Digits required for roundtrip 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 roundtoeven rule) That meant it required relative error ≥ ½ / 10^33 = 500e36 binary128 (113 bits precision), relative error ≤ (½ ulp) / (2^112 ulp) ≈ 96e36 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 50e36 < 96e36 Update: (0.999 ... 8) ε required relative error ≥ 3 * 50e36 = 150e36, thus not possible. Based on relative errors, only 2 patterns possible, (1.000 ...) ε, or (0.999 ...) ε 

06132020, 11:07 AM
Post: #7




RE: (Free42) accuracy of LN1+X
(06122020 08:21 PM)grsbanks Wrote:(06122020 06:51 PM)Albert Chan Wrote: Perhaps it is implemented with binary128 version log1p ? 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 

06142020, 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) { Pauli 

06142020, 12:48 PM
Post: #9




RE: (Free42) accuracy of LN1+X
(06142020 01:09 AM)Paul Dale Wrote: This function is actually fairly easy to implement once its problems are realised. Probably not. decimal128 unable to roundtrip thru binary128 occurs at all ranges, not just tiny x We just don't noticed it when x is bigger ... For example, try this: 1e12 XEQ "LN1+X" → 9.999999999995000000000003333333334e13 We expected result should be rounded(x  x²/2 + x³/3), with last digit = 3 x = 1e12 converted to binary128, we got x' ≈ (10^33+0.08)/10^45, overestimated 0.45 ULP log1p(x') ≈ 9.99999999999500000000000333333333413e13 This example is lucky that roundings compensated the x'x error. You can potientially get into doublerounding 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('1e12')) >>> format(y, '.33e') '9.999999999995000000000003333333334e13' >>> format(next_below(y), '.33e') '9.999999999995000000000003333333332e13' 

06142020, 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.


« Next Oldest  Next Newest »

User(s) browsing this thread: