Squaring the circle: all-MCODE Bessel functions



#15

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


#16

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

Congratulations

Peter


#17

Hi,

my HP48 returns:

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

Best regards,
Jean-Marc


#18

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.


#19

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.


#20

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.


#21

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.


#22

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


#23

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

#24

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.


#25

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.


#26

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

#27

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.

#28

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


Possibly Related Threads...
Thread Author Replies Views Last Post
  Converting Great Circle Navigation from 41C to 42S Bill Triplett 11 2,170 11-13-2013, 07:24 AM
Last Post: Kimberly Thompson
  [41CL] New Extra Functions version Monte Dalrymple 0 683 11-08-2013, 04:32 PM
Last Post: Monte Dalrymple
  HP-41 MCODE: The Last Function - at last! Ángel Martin 0 633 11-08-2013, 05:11 AM
Last Post: Ángel Martin
  HP Prime - Drawing a circle from a program Jean-Michel 6 1,500 11-07-2013, 09:20 AM
Last Post: Tim Wessman
  HP Prime: in need of help with defining functions Alberto Candel 14 2,527 10-27-2013, 10:48 AM
Last Post: Alberto Candel
  HP Prime spreadsheet functions SanS 0 1,289 10-04-2013, 04:23 AM
Last Post: SanS
  prime programming / great circle formulae Geoff Quickfall 15 2,380 09-27-2013, 07:16 PM
Last Post: Kimberly Thompson
  Stats functions on the HP34S Nicholas van Stigt 5 1,266 09-24-2013, 02:45 AM
Last Post: Nick_S
  Trig Functions Howard Owen 11 2,208 09-16-2013, 02:53 PM
Last Post: Fred Lusk
  50g piecewise functions Kurt Fankhauser 6 1,483 09-15-2013, 08:01 PM
Last Post: Kurt Fankhauser

Forum Jump: