▼
Posts: 727
Threads: 43
Joined: Jul 2005
Quote: You may be right. In fact I remember well that fix since it was me who pointed it out to Thomas. I ran into it during my 41Z testing, as I was using Free42 as a reference.
Ángel, do you remember what the test case was? I confirmed that changing the arcsin algorithm back to -i*arcsinh(i*z) fixes the branch cut bug that started the recent thread, but I forget the exact problem that I tried to fix in 1.4.67.
Thanks,
- Thomas
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
perhaps it was arcsin(2.4)?
▼
Posts: 727
Threads: 43
Joined: Jul 2005
No, that's the bug that started *this* thread. I'm asking about the older bug, that Ángel reported to me, that I tried to address in Free42 release 1.4.67. It had something to do with large arguments, but I don't remember what the problem was, exactly.
▼
Posts: 2,247
Threads: 200
Joined: Jun 2005
Posts: 1,253
Threads: 117
Joined: Nov 2005
Yes, it was "large arguments" - meaning large Re and Im parts of course. I think it all started with 50+50i, and from then on up to the precision limit of the calc.
Here's one test case from our emails back then:
Quote:
z=200+200i
asin(sin(z)) is not correct.
Free42 returns 0-162i, whilst the correct answer is -1,62+200i (as returned by the real 42S, and double-checked on the Wolfram site
http://www.wolframalpha.com/input/?i=asin%28sin%28200%2B200i%29%29
Hope this refreshes your memory.
Best,'AM
▼
Posts: 727
Threads: 43
Joined: Jul 2005
Ah, I got it. The problem occurs when z is large while Re(z)<0; in that case, sqrt(z^2+1) approaches -z because of the limited precision, and so ln(z+sqrt(z^2+1)) becomes noise. It seems that the fix is to simply avoid Re(z)<0 altogether, instead computing asinh(z)=-asinh(-z) whenever Re(z)<0. (I guess I should do something similar for acosh.)
▼
Posts: 727
Threads: 43
Joined: Jul 2005
OK, I fixed complex ASINH and ACOSH. ASINH now calculates -arcsinh(-x) if Re(x)<0, and in the ACOSH implementation, I replaced the subexpression sqrt(z^2-1) with sqrt(z+1)*sqrt(z-1).
This fixes complex ASIN and ACOS, too, since they use ASINH and ACOSH underneath.
I uploaded 1.4.70 to my web site and to Android Market. It hasn't appeared on the Market yet; that can take a few hours. I'll update Free42 on iTunes shortly; that will probably take about a week to appear.
▼
Posts: 1,216
Threads: 75
Joined: Jun 2011
Quote:
I uploaded 1.4.70 to my web site
Many thanks for the update, Thomas! :)
Posts: 1,253
Threads: 117
Joined: Nov 2005
Hi THomas, that's basically what I also did in the 41Z code, for ACOSH. I like your method for ASINH more than what I concocted, will look into it as it'll probably save me some bytes.
Free42 Decimal now works great, however I'm still getting different results on Free42 Binary 1.4.70 for:
asin(2.404) = 1.571 - i1.524
asin(2.404+0i) = 1.571 + i1.524
Did you apply the fix to both?
Edited: 19 Sept 2011, 2:51 a.m.
▼
Posts: 727
Threads: 43
Joined: Jul 2005
The fix is automatically applied to both, because it's the same code in both cases, just compiled with a different definition of the 'phloat' type depending on which version is built.
The bug in this case is in the code that calculates the complex square root in the log(x+sqrt(x^2+1)) evaluation; it returns the wrong root because the imaginary part of its argument is -0.
In a word: aargh! I'm sure there are good reasons why IEEE-754 has negative zero, but in this case it (or maybe I should say, the way atan2() treats it) is biting me in the ass. It's easy to fix, but I'll have to check all occurrences of atan2() to see if they have the same potential problem. Fortunately, there aren't all that many...
- Thomas
Posts: 239
Threads: 55
Joined: Sep 2006
Hi Ángel,
And a little OT : WP-34s give a Domaine Error with this example.
Regards,
Miguel
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
The
asin(sin(200+200i))
example?
My 34S is not generating an error -- it isn't getting the correct result though :-( Something to look into.
- Pauli
▼
Posts: 1,253
Threads: 117
Joined: Nov 2005
ok, so how about adding a check in the code for the case:
if Re(z^2+1) =~ Re(z^2) - within the precision of the CPU, that is
and if so, simply use Ln(2z) instead?
This works well for me on the 41Z, can't believe a humble 10-digit mantissa CPU handles it better ?
Cheers,
ÁM.
Posts: 1,216
Threads: 75
Joined: Jun 2011
Quote:
The
asin(sin(200+200i))
example?
My 34S is not generating an error -- it isn't getting the correct result though :-( Something to look into.
Well, at least the emulator gives "Domain Error" on ASIN for this example in DEG/360 mode.
Franz
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
Interesting. This doesn't error in the same place on the real device or the console emulator, however the results are incorrect.
Anyway, I think I've fixed it now. Commit coming soon.
As a bonus, I believe I've also fixed the numeric issues as well.
Not what other problems have I just introduced :-(
Thanks for the clarification.
- Pauli
▼
Posts: 1,216
Threads: 75
Joined: Jun 2011
Quote:
Anyway, I think I've fixed it now. Commit coming soon.
Ok, tried it now with SVN1619 and the error has gone!
Now for something different:
After a long time I've now also downloaded wp34stest.exe again and tried it, but other than a few months ago it just gives a bunch of error messages instead of calculated values.
Here's the output (continuously running):
NaN MAPM Warning: 'm_apm_set_string', Non-digit char found in parse
MAPM Warning: Text =
MAPM Warning: NaN
1.000000000000000000000000000000E+0
Something broken?
Edited: 20 Sept 2011, 7:45 a.m.
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
No idea sorry. Hugh wrote this program to help him test the numeric accuracy of my routines. Not having ready access to any Windows machine, I can't easily run it.
- Pauli
Posts: 3,229
Threads: 42
Joined: Jul 2006
The domain error looks like it was a problem with the hyperbolic sinh of large negative reals. Fixed now.
Introduced when I increased the accuracy of this function for near zero arguments.
- Pauli
Posts: 52
Threads: 2
Joined: Sep 2011
BTW, is there a way to interact with the Free42 from an external program and read back the LCD output? This way you could implement some kind of unit tests to assure correct calculation behavior.
▼
Posts: 185
Threads: 20
Joined: Apr 2007
Quote:
BTW, is there a way to interact with the Free42 from an external program and read back the LCD output? This way you could implement some kind of unit tests to assure correct calculation behavior.
The Free42s output can be sent to the emulated printer, which can be output to a text file. (the PC version)
Then one could have an external program, like Excel, read and parse the text output.
Posts: 727
Threads: 43
Joined: Jul 2005
Any regression testing tool or scripting language that is capable of sending keystrokes to an application can be used to drive Free42. Reading the display can be accomplished by sending Control-C and then reading the clipboard.
Writing all those unit tests will be a bear of a job, though. Besides, once I fix ASIN/ASINH/ACOS/ACOSH, there will be no bugs left in Free42 anyway. ;-)
▼
Posts: 52
Threads: 2
Joined: Sep 2011
You are right, but you don't need to write all regression tests up-front, just add them over time :-)
I used to work for a telecom company and we only started using regression-tests, when it was really easy to write them. So I guess the approach of "GUI-scripting" and emulating mouse-clicks on keys could be quite cumbersome. It would be easier, if you could just supply a text-file with input-keys and the result is validated agains a certain expected output, e.g.:
SIN 3 COS 5 +
1.04853065433
(so each key-press could be separated by a space, expected result is on following line)
Of course "someone" could build a wrapper, to parse such a file, does "GUI-Scripting" and reads back the result.
BTW, I like your implementation on the iPhone, it just rocks!
▼
Posts: 727
Threads: 43
Joined: Jul 2005
Since Free42 just happens to be programmable itself, you could write all the regression tests using plain old HP-42S user code...
Posts: 727
Threads: 43
Joined: Jul 2005
To summarize the whole Free42 complex ASIN saga:
Earlier this year, Ángel Martin reported to me that Free42 returned an incorrect result for arcsin(sin(200+200i)).
I thought I'd fixed that bug by changing the implementation of complex arcsin from -i*arcsinh(i*x) to conj(-i*arcsinh(conj(-i*x))).
I'm kicking myself now for ever making that change. The original algorithm and the revised one are equivalent, except for the edge case of |Re(x)|>1 and Im(x)=0, and the problem that Ángel reported had nothing to do with that edge case. The 1.4.67 change didn't fix the problem with arcsin(sin(200+200i)); it merely moved it to a different part of the complex plane: instead of returning wrong results for 200+200i and -200+200i, it now returned wrong results for 200-200i and -200-200i.
The problem with arcsin(2.404+0i) that two people reported to me via email, and that also started this thread (or, rather, the original thread that got moved to the archive a couple of days ago!), was caused by my non-fix to the arcsin(sin(200+200i)) bug. Undoing that non-fix made the arcsin(2.404+0i) bug go away, but of course it didn't do anything about the arcsin(sin(200+200i)) bug, other than moving it back to its original spot.
The problem with arcsin(sin(200+200i)) turned out to be that my naive implementation of arcsinh(x)=ln(x+sqrt(x^2+1)) fails when |x| is large while Re(x)<0. In that case, the sqrt(x^2+1) part becomes -x+epsilon, and so x+sqrt(x^2+1) becomes noise.
Changing the arcsinh algorithm to -arcsinh(-x) when Re(x)<0 causes the sqrt(x^2+1) part to approach x instead of -x for large |x|, and this avoids the loss of precision in x+sqrt(x^2+1).
So, I fixed the arcsin(2.404+0i) bug I introduced in 1.4.67 by simply reverting back to the 1.4.66 code, and I fixed the arcsin(sin(200+200i)) bug by changing the arcsinh algorithm so as to avoid the loss of precision in x+sqrt(x^2+1) for Re(x)<0.
Note that the reason that Free42 1.4.(67,68,69) did return the correct result for arcsin(2.404) but not for arcsin(2.404+0i) is that the case for real X in ASIN is handled by completely different code than the case for complex X. The "real X" code was right all along.
Also note that there was a bug in ACOS and ACOSH, caused by the same loss of precision issue, in the arccosh code, in the x+sqrt(x^2-1) subexpression for large |x| and Re(x)<0. Nobody reported that bug, and it wasn't affected by the 1.4.67 changes, but that code should be OK now, too.
Whew. :-)
Thank you to Franz Huber and Miguel Toro-Alvarez for contacting me about the arcsin(2.404+0i) bug, and to Ángel Martin for bringing the arcsin(sin(200+200i)) case to my attention, and to everyone who weighed in on the HP Forum!
- Thomas
Edited: 18 Sept 2011, 10:24 p.m.
|