Squaring the circle: all-MCODE Bessel functions

[Ángel Martin]

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

[PeterP]

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

Congratulations

Peter

[Jean-Michel]

Hi,

my HP48 returns:

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

Best regards,
Jean-Marc

[Ángel Martin]

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.

[Jean-Michel]

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.

[Ángel Martin]

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.

[Paul Dale]

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.

[Ángel Martin]

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

[Paul Dale]

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

[Jean-Michel]

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.

[Ángel Martin]

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.

[Paul Dale]

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

[Jean-Michel]

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.

[Ángel Martin]

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