The Museum of HP Calculators

HP Forum Archive 19

[ Return to Index | Top of Index ]

Squaring the circle: all-MCODE Bessel functions
Message #1 Posted by Ángel Martin on 15 Feb 2010, 2:05 p.m.

It's done at last: an All-MCODE implementation of the Bessel functions IBS and JBS, for any real integer and any real argument.

Not an easy feat by my standards, I had to go around numerous issues and limitations to coax the venerable coconut to pull this off! (no doubt my poor programming skills :)

So if you're intrigued by it or interested in 30-year old technology warping its sensible limits, then make sure you check out the forthcoming version of the SandMath module, soon at a TOS near you :)

Also a new mini-white paper on the Bessel functions is ready, including source code listings and even a FOCAL version for complex arguments and complex indexes.

Can you spell J[(1+i),(-1-i)]?

Cheers, "AM

      
Re: Squaring the circle: all-MCODE Bessel functions
Message #2 Posted by PeterP on 15 Feb 2010, 4:55 p.m.,
in response to message #1 by Ángel Martin

do you ever sleep Angel? Your output these last few months has been nothing short of astonishing!

Congratulations

Peter

            
Re: Squaring the circle: all-MCODE Bessel functions
Message #3 Posted by JMBaillard on 16 Feb 2010, 4:58 a.m.,
in response to message #2 by PeterP

Hi,

my HP48 returns:

J1+i(-1-i) = -8.88878423739 + 2.2945020794 i

Best regards, Jean-Marc

                  
Re: Squaring the circle: all-MCODE Bessel functions
Message #4 Posted by Ángel Martin on 16 Feb 2010, 7:29 a.m.,
in response to message #3 by JMBaillard

Excellent answer and correct result.

Here's validation of the ZJBS output then: 1, ENTER^, ZENTER^, ZNEG, ZJBS
-> -8,888784241+J2,294502076

and here's the answer from Wolfram Alfa:

-8.88878423734867397936687003371852995586164928241049073400... + 2.29450207936319254530496624987136389187429554091103480761... i

I was never entirely comfortable using the 48, I guess it shows ;-)

Best wishes,
ÁM

Edited: 16 Feb 2010, 7:35 a.m.

                        
Re: Squaring the circle: all-MCODE Bessel functions
Message #5 Posted by JMBaillard on 17 Feb 2010, 4:54 p.m.,
in response to message #4 by Ángel Martin

Hi,

I've also written a "JNZ" program for the HP-41 and it yields:

-8.888784232 + 2.294502069 i

Roundoff-errors are more important, but my preference still goes to the HP41!

All the Best, Jean-Marc.

                              
Re: Squaring the circle: all-MCODE Bessel functions
Message #6 Posted by Ángel Martin on 18 Feb 2010, 2:29 a.m.,
in response to message #5 by JMBaillard

Let it be publicy said that your 41 Math library is nothing short of sensational, Jean-Marc. You have single-handedly built an impressibe body of programs, all posted here in the museum and obligued reference.

I checked your Bessel functions entry but didn't see the complex treatment. Did I miss it somehow or is not posted yet?

http://www.hpmuseum.org/software/41/41besslj.htm

Yes rounding errors become more apparent when there's only 9 significant digits instead of ... 13, or 15? I'm not even sure what's the precission in the 48X anymore.

Incidentally I found this problem to be less of an issue when doing complex numbers calculations. I'm probably wrong but it somehow seems to be a more forgiving environment. I now realize I should have used the 13-digit precision mainframe routines in some of my 41Z functions, and will revise them eventually.

The exception to this "benign rule" has been the Lambert W function. Using the ZSOLVE program I can't always find the right root for the functional equation that defines W(z). I also tried with your program "ROOTZ", but the same issue occurs.

http://www.hpmuseum.org/software/41/41cmpxf.htm

Have you tried using ROOTZ to calculate the Lamber W function with it?
The equation to solve is Exp[W(z)]*W(z)-Z0 = 0, and I'm using the following code to define it:

LBL "*W"
ZEXP
LASTZ
Z*
ZRCL
4
Z-
END

where the function's complex argument Z0 is stored in complex register ZR4, and the evaluated guesses of the equation z are in the complex stack register "Z" upon entry.

It fares well for "small" Z0 values (say 1+i) but starts acting up when the argument "grows" - say (5+5i).

In truth, I wonder to what extent one can use the secant method in the complex plane - where the mere concept doesn't hold up very gracefully.

Cheers, ÁM

Edited: 18 Feb 2010, 2:53 a.m.

                                    
Re: Squaring the circle: all-MCODE Bessel functions
Message #7 Posted by Paul Dale on 18 Feb 2010, 4:22 a.m.,
in response to message #6 by Ángel Martin

I've done complex Lambert's W for my 20b scientific firmware :)

My references to the algorithm sources are:

	http://en.wikipedia.org/wiki/Lambert_W_function
	http://ioannis.virtualcomposer2000.com/math/LWExpansion.html
	http://ioannis.virtualcomposer2000.com/math/LWCalculating.html

These are quite probably not up to date.

I used the algorithms from the wikipedia page.

Start with the approximation under 7.2 and then refine it using Halley's method under 7.1 (I limited my maximum iterations to 20).

- Pauli

Edited: 21 Feb 2010, 4:34 a.m.

                                          
Re: Squaring the circle: all-MCODE Bessel functions
Message #8 Posted by Ángel Martin on 21 Feb 2010, 1:58 p.m.,
in response to message #7 by Paul Dale

Thanks for the references Paul. So you wrote the firmware extensions for the 20B in machine code as well? Hmm, that makes me consider purchasing the machine, despite I'm sure I won;t understand half of it (never been a financial guy).

Cheers, 'AM

                                                
Re: Squaring the circle: all-MCODE Bessel functions
Message #9 Posted by Paul Dale on 21 Feb 2010, 3:53 p.m.,
in response to message #8 by Ángel Martin

I've written everything so far in C. The 20b uses an ARM CPU and compilers for that are easy to get. I'm using HP-GCC.

- Pauli

                                    
Re: Squaring the circle: all-MCODE Bessel functions
Message #10 Posted by JMBaillard on 18 Feb 2010, 4:30 p.m.,
in response to message #6 by Ángel Martin

Hi Angel,

Many thanks for your appreciations. The complex version "JNZ" is missing in the museum, all the more that I wrote it after your post on J1+i(-1-i)

-Here is the listing:

Data Registers: R00 thru R10: temp

Flag: F01 CF 01 to compute Jn(z) SF 01 to compute In(z)

Subroutines: "GAMZ" or "GAMZ+" ( cf "Gamma Function For the HP-41" § 1°)-h) ) "Z^2" "Z*Z" "Z^Z" "Z/Z" ( cf "Complex Functions for the HP-41" )

"Z^2" "Z*Z" ... and so on ... may replaced by M-Code routines ( cf "A few M-Code Routines for the HP-41" )

01 LBL "JNZ" 02 STO 01 03 RDN 04 STO 02 05 RDN 06 STO 03 07 X<>Y 08 STO 04 09 2 10 ST/ 01 11 ST/ 02 12 RCL 02 13 RCL 01 14 XEQ "Z^2" 15 STO 09 16 X<>Y 17 STO 10 18 CLX 19 STO 06 20 STO 08 21 SIGN 22 STO 00 23 STO 05 24 STO 07 25 LBL 01 26 RCL 10 27 RCL 09 28 RCL 06 29 RCL 05 30 XEQ "Z*Z" 31 RCL 00 32 FC? 01 33 CHS 34 ST/ Z 35 / 36 RCL 03 37 RCL 00 38 + 39 RCL 04 40 X<>Y 41 XEQ "Z/Z" 42 STO 05 43 X<>Y 44 STO 06 45 RCL 08 46 + 47 STO 08 48 LATSX 49 - 50 X<>Y 51 RCL 07 52 + 53 STO 07 54 LASTX 55 - 56 R-P 57 ISG 00 58 CLX 59 X#0? 60 GTO 01 61 RCL 02 62 RCL 01 63 RCL 04 64 RCL 03 65 XEQ "Z^Z" 66 RCL 08 67 RCL 07 68 XEQ "Z*Z" 69 STO 07 70 X<>Y 71 STO 08 72 RCL 04 73 RCL 03 74 1 75 + 76 XEQ "GAMZ" or XEQ "GAMZ+" 77 RCL 08 78 RCL 07 79 R^ 80 R^ 81 XEQ "Z/Z" 82 END

( 125 bytes / SIZE 011 )

STACK INPUTS OUTPUTS

T b / Z a / Y y Im [ f(z) ] X x Re [ f(z) ]

n = a + i b where z = x + i y and f(z) = Jn(z) if flag F01 is clear, f(z) = In(z) if flag F01 is set.

Examples: n = 1 + 4 i , z = 2 + 3 i

• CF 01

4 ENTER^ 1 ENTER^ 3 ENTER^ 2 XEQ "JNZ" >>>> J1+4i ( 2+3.i ) = 0.342267397 - 0.424349408 i ( Execution time = 33s )

• SF 01

4 ENTER^ 1 ENTER^ 3 ENTER^ 2 XEQ "JNZ" >>>> I1+4i ( 2+3.i ) = 1.391608538 + 0.315041667 i ( Execution time = 33s )

With respect to the Lambert W functions, "ROOTZ" seems to work well provided we choose 2 initial guesses close to Ln(z) [ if z # 0 ]

For example, if z = 5 + 5i , try z0 = 1.9 + 0.8 i and z1 = 2 + 0.9 i

With the subroutine

LBL "W" STO 11 X<>Y STO 12 X<>Y XEQ "E^Z" RCL 12 RCL 11 XEQ "Z*Z" RCL 10 RCL 09 XEQ "Z-Z" RTN

"ROOTZ" returns w(5+5i) = 1.501424530 + 0.477488825 i after 6 or 7 iterations.

But Newton's method has a better convergence, and it is used in the following routine:

LBL "LAMBZ" STO O1 X<>Y STO 02 R-P X=0? GTO 01 RCL 02 RCL 01 XEQ "LNZ" GTO 02 LBL 01 CLST LBL 02 STO 03 X<>Y STO 04 LBL 03 RCL 04 RCL 03 XEQ "E^Z" STO 05 X<>Y STO 06 X<>Y RCL 04 RCL 03 XEQ "Z*Z" ST+ 05 X<>Y ST+ 06 X<>Y RCL 02 RCL 01 XEQ "Z-Z" RCL 06 RCL 05 XEQ "Z/Z" CHS RCL 03 + STO 03 LASTX - X<>Y CHS RCL 04 + STO 04 LASTX - R-P VIEW 03 E-9 X<Y? GTO 03 RCL 04 RCL 03 END

and for example, 5 ENTER^ XEQ "LAMBZ" gives

w(5+5i) = 1.501424530 + 0.477488825 i

( the same result as "ROOTZ" with one iteration less )

Likewise w(7+9i) = 1.796464032 + 0.591616466 i

"RTZ4" should give even faster results...

My contribution to the HP-41 programs has not stopped: Since 2007, I've sent about 80 ( yes, eighty ) new or updated pages and I'm surprised not to see them in the MoHPC!

Perhaps will I create my own website, but it's impossible for the moment: my access provider has many problems since the end of October and it still does not work normally!

Best Regards, Jean-Marc.

                                          
Re: Squaring the circle: all-MCODE Bessel functions
Message #11 Posted by Ángel Martin on 21 Feb 2010, 4:04 a.m.,
in response to message #10 by JMBaillard

Hi Jean-Marc,

Excellent news, thanks for looking into this. Yes, ZSOLVE now works fine for Lambert W using the better initial estimation as you suggest, using the logarithm instead of my crude choice.

I'm impressed by the short execution times of your programs in general, also JNZ. Not bad at all considering you're not using MCODE routines. Are you programming the standard definition doing the summation, or other approximation perhaps?

Shame to hear that your newer contributions aren't posted yet - I hope that gets corrected soon! BTW, I didn't find the reference you gave on the "Some MCODE routines for the 41" - is that another one not posted as well? I was curious to see if you used the standard mainframe routines or the 13-bit extended precision to program the complex functions.

Thanks again for your inputs.

Edited: 21 Feb 2010, 4:17 a.m.

                                                
Re: Squaring the circle: all-MCODE Bessel functions
Message #12 Posted by Paul Dale on 21 Feb 2010, 4:36 a.m.,
in response to message #11 by Ángel Martin

I've also noticed the software library hasn't been updated for ages. I've only submitted a few programs but they were a long time ago.

- Pauli

                                                
Re: Squaring the circle: all-MCODE Bessel functions
Message #13 Posted by JMBaillard on 21 Feb 2010, 5:39 p.m.,
in response to message #11 by Ángel Martin

Hi Angel,

I've just sent you the page on M-code routines
and since I don't know how to use 13-bit precision
I would need some tutorial on this point...

Best Regards,
Jean-Marc.

            
Re: Squaring the circle: all-MCODE Bessel functions
Message #14 Posted by Ángel Martin on 17 Feb 2010, 1:48 a.m.,
in response to message #2 by PeterP

Thanks Peter, it was just rounding off some loose ends in the SandMath - but I admit it's taken me farther (and longer) than what I originally thought.

Now I'm preparing some mini-papers to document the highlights of the module. As addictive as 41 MCODE is, it'll be good to rest for a while!

Cheers, 'AM


[ Return to Index | Top of Index ]

Go back to the main exhibit hall