Lommel Functions for the HP-41

Overview
 

-"LOML1" & "LOML2" calculate the Lommel functions of the 1st & 2nd kind.
-They use the formulae:
 

    s⁽¹⁾_m,n(x)  = x^m+1 / [ (m+1)² - n² ]  ₁F₂ ( 1 ; (m-n+3)/2 , (m+n+3)/2 ; -x²/4 )

    s⁽²⁾_m,n(x)  = x^m+1 / [ (m+1)² - n² ]  ₁F₂ ( 1 ; (m-n+3)/2 , (m+n+3)/2 ; x²/4 )
                      + 2^m+n-1 Gam(n) Gam((m+n+1)/2) x ^-n / Gam((-m+n+1)/2)  ₀F₁ ( ; 1-n ; -x²/4 )
                      + 2^m-n-1 Gam(-n) Gam((m-n+1)/2) xⁿ / Gam((-m-n+1)/2)  ₀F₁ ( ; 1+n ; -x²/4 )

   where  _pF_q is the generalized hypergeometric function and Gam is Euler Gamma function.
 

Program Listing
 

Data Registers:   R00 thru R09: temp
Flag:  F01
Subroutines: /

-The M-code routine "HGF+" may be replaced by XEQ "GHGF" ( or "0F1" and "1F2" ) but in this case, replace R00 by another unused register ( R10 )
-Likewise, 1/G+ may be replaced by XEQ "1/G+" but store the intermediate result in another register ( XEQ "1/G+" does not save Y-register )
 ( cf "Hypergeometric Functions " & "Gamma Function" )
 
 

01  LBL "LOML2"   02  CF 01   03  XEQ 01   04  STO 09   05  XEQ 00   06  STO 08             07  RCL 05   08  CHS   09  STO 05   10  XEQ 00   11  RCL 08   12  +   13  2   14  ST/ Y   15  RCL 07   16  Y^X   17  *   18  RCL 09   19  +   20  RTN   21  LBL 00

22  1   23  RCL 05   24  -   25  STO 01             26  CLST   27  SIGN   28  RCL 06   29  CHS   30  HGF+   31  RCL 05   32  1   33  +   34  RCL 07   35  -   36  2   37  /   38  1/G+   39  *   40  2   41  RCL 05   42  Y^X

43  *   44  RCL 05   45  1/G+   46  /   47  RCL 04   48  RCL 05   49  +   50  2   51  /   52  1/G+   53  /   54  RCL 00   55  RCL 05   56  Y^X   57  /   58  RTN   59  LBL "LOML1"   60  SF 01   61  LBL 01   62  STO 00             63  RDN

64  STO 03   65  STO 05   66  CHS   67  X<>Y   68  STO 04   69  STO 07   70  ST+ 03   71  +   72  3   73  ST+ 03   74  +   75  2   76  ST/ 03   77  /   78  STO 02             79  1   80  STO 01   81  ST+ 04   82  LASTX   83  RCL 00   84  LASTX

85  /   86  X^2   87  STO 06   88  FS? 01   89  CHS   90  HGF+   91  RCL 00   92  RCL 04   93  Y^X   94  *   95  RCL 04   96  X^2   97  RCL 05             98  X^2   99  - 100  / 101  END

 
    ( 138 bytes / SIZE 010 )
 
 

STACK	INPUTS	OUTPUTS
Z	m	/
Y	n	/
X	x	f(x)

   With   f(x) =  s⁽¹⁾_m,n(x)  or  s⁽²⁾_m,n(x)

Examples:

   2   SQRT
   3   SQRT
  PI  XEQ "LOML1"  >>>>   s⁽¹⁾_{sqrt(2),sqrt(3)}(PI) = 3.003060384            ---Execution time =5s---

   2   SQRT
   3   SQRT
  PI  XEQ "LOML2"  >>>>   s⁽²⁾_{sqrt(2),sqrt(3)}(PI) = 9.048798662            ---Execution time =23s---

Variant:

-Replace line 88 by FC? 02
-Delete lines 59-60
-Replace liines 01-02  by  LBL "LOML"   FC? 01  GTO 01

-Then,   CF 02  XEQ "LOML"  will return  s⁽¹⁾_m,n(x)
  and      SF 02  XEQ "LOML"  will return  s⁽²⁾_m,n(x)

Notes:

-There will be a DATA ERROR if  m+1 = +/-n
-Since the formulas involve  Gam(n) & Gam(-n),  "LOML2" will not work if  n is an integer.
-There are other similar limitations.
 

Reference:

[1]  http://mathworld.wolfram.com/LommelFunction.html