Re: Free42 - A Complex glitch?



#25

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


#26

perhaps it was arcsin(2.4)?


#27

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.


#28

Ah ok!!

#29

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


#30

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.)


#31

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.


#32

Quote:
I uploaded 1.4.70 to my web site

Many thanks for the update, Thomas! :)
#33

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.


#34

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

#35

Hi Ángel,

And a little OT : WP-34s give a Domaine Error with this example.

Regards,

Miguel


#36

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


#37

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.

#38

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


#39

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


#40

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.


#41

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

#42

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

#43

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.


#44

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.

#45

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. ;-)


#46

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!


#47

Since Free42 just happens to be programmable itself, you could write all the regression tests using plain old HP-42S user code...

#48

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.


Possibly Related Threads…
Thread Author Replies Views Last Post
  HP Prime: complex numbers in CAS. Alberto Candel 1 1,941 12-06-2013, 02:36 PM
Last Post: parisse
  [HP Prime] Plots containing complex numbers bug? Chris Pem10 7 3,703 12-05-2013, 07:40 AM
Last Post: cyrille de Brébisson
  Complex Number Entry on Prime Jeff O. 19 5,305 11-16-2013, 12:34 PM
Last Post: Jeff O.
  HP Prime complex results Javier Goizueta 0 1,008 10-06-2013, 12:59 PM
Last Post: Javier Goizueta
  HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 4,343 09-30-2013, 03:44 AM
Last Post: From Hong Kong
  [HP-Prime xcas] operations with complex numbers + BUGs + Request CompSystems 9 3,560 09-08-2013, 10:40 PM
Last Post: CompSystems
  Elliptic integrals of 1st and 2nd kind calculated by complex agm Gjermund Skailand 3 1,521 06-29-2013, 03:39 PM
Last Post: Gjermund Skailand
  HP-11C and complex numbers Antlab 5 2,248 06-28-2013, 08:59 AM
Last Post: Antlab
  71B Complex User Defined Functions in Calc Mode? Michael Burr 0 960 03-18-2013, 09:38 PM
Last Post: Michael Burr
  HP 34S complex # problem Richard Berler 2 1,241 02-17-2013, 11:01 PM
Last Post: Richard Berler

Forum Jump: