Arity icon indicating copy to clipboard operation
Arity copied to clipboard

Precision issues

Open ChinhQT opened this issue 2 years ago • 26 comments

Hi, I don't know why this one is wrong. The Arity's result is tanh^-1(2%^9) = 5.5511151E-16 But the correct result is tanh^-1(2%^9) = 5.12E-16.

ChinhQT avatar May 18 '23 05:05 ChinhQT

are you sure? 2%^9 equals 5.12E-16 Should atanh of this value be the same?

woheller69 avatar May 18 '23 06:05 woheller69

tanh and atanh (2%^9) = 5.12E-16. They have the same result. I've just tested on other app. It still has the same result. You can check them again.

ChinhQT avatar May 18 '23 06:05 ChinhQT

I checked the formula and it is correct. So I guess these are precision issues of Java's Math library.

woheller69 avatar May 18 '23 06:05 woheller69

Yes, looks like Java's library has some general issues with trigonometry. For example, sin(11π) gives 4.899E-15, sin(2000π) is evaluated as -6.428E-13 (while both expressions obviously should be equal to 0)

Fontan030 avatar May 18 '23 10:05 Fontan030

Java double seems to be good for about 15 digits...

woheller69 avatar May 18 '23 10:05 woheller69

Are you sure ? I don't think it's an error of Java's library. The ClevCal on Android, it works well. The result is correct. I don't why.

ChinhQT avatar May 18 '23 14:05 ChinhQT

Maybe ClevCal uses some other arithmetic library.

The formula used by ArityEngine is: 0.5 * Math.log(1. + (x + x) / (1. - x)); which equals 0.5 * Math.log(1. + x) / (1. - x));

https://en.wikipedia.org/wiki/Inverse_hyperbolic_functions

woheller69 avatar May 18 '23 15:05 woheller69

btw, Java has already problems subtracting 56.3 - 55.7. It gives 0.59999999 I recently fixed that by using BigDecimal for subtractions

woheller69 avatar May 18 '23 15:05 woheller69

With the formula 0.5 * Math.log(1. + x) / (1. - x)) x = 2%^9. The result is wrong. When you use x = 2%^8. The Arity works well. The result is correct. I don't know why.

ChinhQT avatar May 19 '23 01:05 ChinhQT

well I guess it is the 1 + 2%^9 or 1 - 2%^9. One or both of them probably exceeds the possibilty of "double"

woheller69 avatar May 19 '23 03:05 woheller69

Do you have any suggestion to resolve this issue, please ?

ChinhQT avatar May 20 '23 03:05 ChinhQT

No. For high-precision stuff use another calculator.

woheller69 avatar May 20 '23 03:05 woheller69

next version will have a better check for pi multiples, so sin(11PI) and sin(2000PI) should be 0 But that is all I can do.

woheller69 avatar May 20 '23 05:05 woheller69

well I guess it is the 1 + 2%^9 or 1 - 2%^9. One or both of them probably exceeds the possibilty of "double"

Yep, I looked into this and it's actually a known issue.
libc (llvm) gets around this by using a taylor series approximation for small values, the code boils down to 1/11 x^11 + 1/9 x^9 + 1/7 x^7 + 1/5 x^5 + 1/3 x^3 + x. glibc (gnu) uses log1p, which calculates ln(1 + x) with greater precision around x=0.
Neither clang (llvm), gcc (gnu), nor msvc (microsoft) optimized 0.5 * log(1 + (x + x) / (1 - x)) to 0.5 * log1p((x + x) / (1 - x)), so I doubt java would as well.

If java has log1p, it should be an easy replacement for small input values, and if not, the taylor approximation also works.

These methods worked with the given value, even on 32-bit float.

supsm avatar Oct 13 '23 02:10 supsm

I am using my own fork https://github.com/woheller69/ArityEngine Will apply the fix there :-) Maybe we should add a test case also

woheller69 avatar Oct 13 '23 04:10 woheller69

if 0.5 * log1p((x + x) / (1 - x)) should be used:

Why are you proposing

 0.5 * Math.log1p(twox + twox * x / (1 - x));

Shouldn't it be simply

 0.5 * Math.log1p(twox / (1 - x));

woheller69 avatar Oct 13 '23 05:10 woheller69

They're mathematically the same thing. I'm not exactly sure why the gcc wizards chose this, but it must have been purposeful somehow as they had the other option for another case. Besides, they're some of the world's smartest people working on critical infrastructure, so there's plenty of reason for it to be as optimized as possible.

I think it might be because for exceedingly small values of x, twox + twox * x / (1 - x) guarantees at least twox (x + x is just an increase in floating point exponent with no effect on precision, and operator precedence dictates that the addition is done last; twox * x / (1 - x) can never be negative), which is an alright approximation at those values. twox / (1 - x), on the other hand, has no such guarantee and could easily be rounded down. This is just my personal opinion, perhaps the authors had a different rationale.
Either way, their approach should translate to java pretty effectively as the code is already "generic" and not tailored to a specific architecture or platform.

supsm avatar Oct 13 '23 06:10 supsm

Are there similar tweaks for things like

sin(11Pi) or sin(2000Pi) ? (See above)

For exact multiples of PI I fixed it already in my fork. But maybe there are better solutions to improve precision

Modulo causes issues with complex numbers...

woheller69 avatar Oct 13 '23 06:10 woheller69

I don't believe so. Neither gcc nor clang fixes this issue (they give the same result as java), and in fact the input value is simply not representable very well using double-precision floating point. The closest representable value is 6283.18530718, and its sin value is pretty close to what java and C output (2.6666141e-13). This value is slightly different from reported above, so I'm not sure if there was some input rounding somewhere. I used an input value of 6283.1853071795864 as 2000pi in my testing.

supsm avatar Oct 13 '23 06:10 supsm

Maybe this library is an option?

https://github.com/eobermuhlner/big-math

woheller69 avatar Oct 13 '23 06:10 woheller69

But probably that does not help either, as the input value does not come as BigDecimal. The whole Arity Engine would have to be changed...

woheller69 avatar Oct 13 '23 06:10 woheller69

@supsm I added your modifications and some test cases to my fork of Arity engine: https://github.com/woheller69/ArityEngine

Here is an updated version of Arity. tanh and atanh (2%^9) now yield 5.12E-16

-uninstall -remove .zip and install

app-release.apk.zip

woheller69 avatar Oct 14 '23 13:10 woheller69

This should fix some more issues with sin and cos calculations for multiples of PI or PI/2 https://github.com/woheller69/ArityEngine/commit/fab02d62e38fe5af723891b6c7d273cc1b7efca7

app-release.apk.zip

woheller69 avatar Oct 14 '23 19:10 woheller69

I found a similar problem:

4÷3×π×(((6371+0.1)×1000)^3−(6371×1000)^3×2÷1000. Dividing this formula by 1000:

4÷3×π×(((6371+0.1)×1000)^3−(6371×1000)^3×2÷1000÷1000 does not result in a lower exponent, only the mantissa changes

Murmele avatar Jun 05 '24 09:06 Murmele

If you want to divide the whole formula by 1000 you should put it in brackets

(4÷3×π×(((6371+0.1)×1000)^3−(6371×1000)^3×2÷1000)÷1000

woheller69 avatar Jun 05 '24 09:06 woheller69

Sorry I completely overseen that one bracket is not closed

Murmele avatar Jun 05 '24 11:06 Murmele