### Program Code for Calculation of Laguerre Polynomials(specially adapted for large indices)

```/********************************************************************************
*  This PL/1 subroutine calculates of the Laguerre polynomials, based on the  recursive
*  relationships (Abramowitz and Stegun, Handbook of Mathematical Functions):
*
*  L0q(x)=1
*  L1q(x)=q+1-x
*  Lpq(x)=1/p.[(2p+q-1-x).Lp-1q(x) -(p-1+q).Lp-2q(x)]
*
*
*  Input parameters are:
*    X: argument of the Laguerre polynomial
*    N0: principal quantum number
*    L0: angular momentum quantum number
*     ( note that  p=N0-L0-1;   q=2L0+1 )
*
*  Output parameters are:
*    L: mantissa of Laguerre-polynomial
*    LEX: exponent of Laguerre-polynomial
*
*  Required subroutines:
*        (to avoid numerical overflow or underflow for large indices).
*
*
*   This program has worked well and correct for me in the past (using IBM PL/1 compilers),
*   but I can't guarantee this for all compilers, especially if the code is translated
*   into a different language.
*
*   If there are any problems, please let me know and I mention the issue here.
*   (for contact details see  http://www.plasmaphysics.org.uk/feedback.htm ).
*
*                                               Thomas Smid, May 2007
*****************************************************************************/

LAGUERRE : PROCEDURE(X,N,NLL) ;
DEFAULT RANGE(A:Z) DEC FLOAT VALUE (DEC FLOAT(16)) ;
DCL (L(*,*),NLL(*)) CTL DEC FLOAT(16) ;
DCL (N,I,J,IJ,J1,J2) DEC FIXED(7,0) ;

ON UNDERFLOW ;

ALLOCATE L(N,N) ;

L=0 ;

DO I=1 TO N ;
L(I,I)=1 ;
END ;

DO I=1 TO N-1 ;
L(I,I+1)=2*I-X ;
END ;

DO I=1 TO N-2 ;
DO J=I+2 TO N;
J1=J-1; J2=J-2; IJ=I+J-2 ;
L(I,J)=((2*J-2-X)*L(I,J1)-IJ*L(I,J2))/(J-I) ;
END ;
END ;

DO I=1 TO N  ;
NLL(I)=L(I,N) ;
END ;

FREE L ;

END LAGUERRE ;

/***********************************************************************************
*  The subroutines NUM, ADD, MULT, DIV enable basic arithmetic operations of numbers
*  with exponents up to ±1015.
*  NUM separates usual floating point numbers (X) into a mantissa (X) and exponent (XEX),
*  and ADD, MULT, DIV provide the basic arithmetic operations of these representations,
*  where XM,XEX,YM,YEX are the two input numbers and XYM, XYEX the result.
***********************************************************************************/

NUM: PROCEDURE(X,XEX) ;
DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
DCL X DEC FLOAT(16) ;
DCL XEX DEC FIXED(15,0) ;
LAX=LOG10(ABS(X)) ;
XEX=FLOOR(LAX) ;
X=10**(LAX-XEX)*SIGN(X) ;
END NUM ;

DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
DCL (XM,YM,XYM) DEC FLOAT(16) ;
DCL (XEX,YEX,XYEX) DEC FIXED(15,0) ;
IF XM=0&YM=0 THEN DO;
XYM=0; XYEX=0;
RETURN;
END ;
IF XM=0 THEN DO;
XYM=YM; XYEX=YEX;
RETURN;
END ;
IF YM=0 THEN DO;
XYM=XM; XYEX=XEX;
RETURN;
END ;

IF XEX>=YEX THEN DO ;
XYEX=XEX ;
XYM=YEX-XEX; XYME=XYM;
IF XYM<-60 THEN XYM=0; ELSE XYM=10**XYME ;
XYM=XM+YM*XYM ;
END ;
ELSE DO ;
XYEX=YEX;
XYM=XEX-YEX; XYME=XYM;
IF XYM<-60 THEN XYM=0; ELSE XYM=10**XYME ;
XYM=XM*XYM+YM ;
END ;
IF ABS(XYM)>=10 THEN DO ;
XYM=XYM/10; XYEX=XYEX+1 ;
END ;
IF ABS(XYM)<1 THEN DO;
XYM=XYM*10; XYEX=XYEX-1;
END ;

MULT: PROCEDURE(XM,XEX,YM,YEX,XYM,XYEX) ;
DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
DCL (XM,YM,XYM) DEC FLOAT(16) ;
DCL (XEX,YEX,XYEX) DEC FIXED(15,0) ;
IF XM=0|YM=0 THEN DO ;
XYM=0; XYEX=0;
RETURN;
END ;
XYM=XM*YM ;
XYEX=XEX+YEX ;
IF ABS(XYM)>=10 THEN DO ;
XYM=XYM/10; XYEX=XYEX+1 ;
END ;
END MULT ;

DIV: PROCEDURE(XM,XEX,YM,YEX,XYM,XYEX) ;
DEFAULT RANGE(A:Z) DEC FLOAT VALUE(DEC FLOAT(16)) ;
DCL (XM,YM,XYM) DEC FLOAT(16) ;
DCL (XEX,YEX,XYEX) DEC FIXED(15,0) ;
IF XM=0 THEN DO;
XYM=0; XYEX=0;
RETURN;
END ;
XYM=XM/YM ;
XYEX=XEX-YEX ;
IF ABS(XYM)<1 THEN DO ;
XYM=XYM*10; XYEX=XYEX-1 ;
END ;
END DIV ;

```