relibc_openlibm/slatec/xnrmp.f

*DECK XNRMP
      SUBROUTINE XNRMP (NU, MU1, MU2, SARG, MODE, SPN, IPN, ISIG,
     1   IERROR)
C***BEGIN PROLOGUE  XNRMP
C***PURPOSE  Compute normalized Legendre polynomials.
C***LIBRARY   SLATEC
C***CATEGORY  C3A2, C9
C***TYPE      SINGLE PRECISION (XNRMP-S, DXNRMP-D)
C***KEYWORDS  LEGENDRE FUNCTIONS
C***AUTHOR  Lozier, Daniel W., (National Bureau of Standards)
C           Smith, John M., (NBS and George Mason University)
C***DESCRIPTION
C
C        SUBROUTINE TO CALCULATE NORMALIZED LEGENDRE POLYNOMIALS
C        (DXNRMP is double-precision version)
C        XNRMP calculates normalized Legendre polynomials of varying
C        order and fixed argument and degree. The order MU and degree
C        NU are non-negative integers and the argument is real. Because
C        the algorithm requires the use of numbers outside the normal
C        machine range, this subroutine employs a special arithmetic
C        called extended-range arithmetic. See J.M. Smith, F.W.J. Olver,
C        and D.W. Lozier, Extended-Range Arithmetic and Normalized
C        Legendre Polynomials, ACM Transactions on Mathematical Soft-
C        ware, 93-105, March 1981, for a complete description of the
C        algorithm and special arithmetic. Also see program comments
C        in XSET.
C
C        The normalized Legendre polynomials are multiples of the
C        associated Legendre polynomials of the first kind where the
C        normalizing coefficients are chosen so as to make the integral
C        from -1 to 1 of the square of each function equal to 1. See
C        E. Jahnke, F. Emde and F. Losch, Tables of Higher Functions,
C        McGraw-Hill, New York, 1960, p. 121.
C
C        The input values to XNRMP are NU, MU1, MU2, SARG, and MODE.
C        These must satisfy
C          1. NU .GE. 0 specifies the degree of the normalized Legendre
C             polynomial that is wanted.
C          2. MU1 .GE. 0 specifies the lowest-order normalized Legendre
C             polynomial that is wanted.
C          3. MU2 .GE. MU1 specifies the highest-order normalized Leg-
C             endre polynomial that is wanted.
C         4a. MODE = 1 and -1.0 .LE. SARG .LE. 1.0 specifies that
C             Normalized Legendre(NU, MU, SARG) is wanted for MU = MU1,
C             MU1 + 1, ..., MU2.
C         4b. MODE = 2 and -3.14159... .LT. SARG .LT. 3.14159... spec-
C             ifies that Normalized Legendre(NU, MU, COS(SARG)) is want-
C             ed for MU = MU1, MU1 + 1, ..., MU2.
C
C        The output of XNRMP consists of the two vectors SPN and IPN
C        and the error estimate ISIG. The computed values are stored as
C        extended-range numbers such that
C             (SPN(1),IPN(1))=NORMALIZED LEGENDRE(NU,MU1,X)
C             (SPN(2),IPN(2))=NORMALIZED LEGENDRE(NU,MU1+1,X)
C                .
C                .
C             (SPN(K),IPN(K))=NORMALIZED LEGENDRE(NU,MU2,X)
C        where K = MU2 - MU1 + 1 and X = SARG or COS(SARG) according
C        to whether MODE = 1 or 2. Finally, ISIG is an estimate of the
C        number of decimal digits lost through rounding errors in the
C        computation. For example if SARG is accurate to 12 significant
C        decimals, then the computed function values are accurate to
C        12 - ISIG significant decimals (except in neighborhoods of
C        zeros).
C
C        The interpretation of (SPN(I),IPN(I)) is SPN(I)*(IR**IPN(I))
C        where IR is the internal radix of the computer arithmetic. When
C        IPN(I) = 0 the value of the normalized Legendre polynomial is
C        contained entirely in SPN(I) and subsequent single-precision
C        computations can be performed without further consideration of
C        extended-range arithmetic. However, if IPN(I) .NE. 0 the corre-
C        sponding value of the normalized Legendre polynomial cannot be
C        represented in single-precision because of overflow or under-
C        flow. THE USER MUST TEST IPN(I) IN HIS/HER PROGRAM. In the case
C        that IPN(I) is nonzero, the user should try using double pre-
C        cision if it has a wider exponent range. If double precision
C        fails, the user could rewrite his/her program to use extended-
C        range arithmetic.
C
C        The interpretation of (SPN(I),IPN(I)) can be changed to
C        SPN(I)*(10**IPN(I)) by calling the extended-range subroutine
C        XCON. This should be done before printing the computed values.
C        As an example of usage, the Fortran coding
C              J = K
C              DO 20 I = 1, K
C              CALL XCON(SPN(I), IPN(I),IERROR)
C              IF (IERROR.NE.0) RETURN
C              PRINT 10, SPN(I), IPN(I)
C           10 FORMAT(1X, E30.8 , I15)
C              IF ((IPN(I) .EQ. 0) .OR. (J .LT. K)) GO TO 20
C              J = I - 1
C           20 CONTINUE
C        will print all computed values and determine the largest J
C        such that IPN(1) = IPN(2) = ... = IPN(J) = 0. Because of the
C        change of representation caused by calling XCON, (SPN(I),
C        IPN(I)) for I = J+1, J+2, ... cannot be used in subsequent
C        extended-range computations.
C
C        IERROR is an error indicator. If no errors are detected,
C        IERROR=0 when control returns to the calling routine. If
C        an error is detected, IERROR is returned as nonzero. The
C        calling routine must check the value of IERROR.
C
C        If IERROR=112 or 113, invalid input was provided to XNRMP.
C        If IERROR=101,102,103, or 104, invalid input was provided
C        to XSET.
C        If IERROR=105 or 106, an internal consistency error occurred
C        in XSET (probably due to a software malfunction in the
C        library routine I1MACH).
C        If IERROR=107, an overflow or underflow of an extended-range
C        number was detected in XADJ.
C        If IERROR=108, an overflow or underflow of an extended-range
C        number was detected in XC210.
C
C***SEE ALSO  XSET
C***REFERENCES  Smith, Olver and Lozier, Extended-Range Arithmetic and
C                 Normalized Legendre Polynomials, ACM Trans on Math
C                 Softw, v 7, n 1, March 1981, pp 93--105.
C***ROUTINES CALLED  XADD, XADJ, XERMSG, XRED, XSET
C***REVISION HISTORY  (YYMMDD)
C   820712  DATE WRITTEN
C   881020  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
C   901019  Revisions to prologue.  (DWL and WRB)
C   901106  Changed all specific intrinsics to generic.  (WRB)
C           Corrected order of sections in prologue and added TYPE
C           section.  (WRB)
C           CALLs to XERROR changed to CALLs to XERMSG.  (WRB)
C   920127  Revised PURPOSE section of prologue.  (DWL)
C***END PROLOGUE  XNRMP
      INTEGER NU, MU1, MU2, MODE, IPN, ISIG
      REAL SARG, SPN
      DIMENSION SPN(*), IPN(*)
      REAL C1,C2,P,P1,P2,P3,S,SX,T,TX,X,RK
C CALL XSET TO INITIALIZE EXTENDED-RANGE ARITHMETIC (SEE XSET
C LISTING FOR DETAILS)
C***FIRST EXECUTABLE STATEMENT  XNRMP
      IERROR=0
      CALL XSET (0, 0, 0.0, 0,IERROR)
      IF (IERROR.NE.0) RETURN
C
C        TEST FOR PROPER INPUT VALUES.
C
      IF (NU.LT.0) GO TO 110
      IF (MU1.LT.0) GO TO 110
      IF (MU1.GT.MU2) GO TO 110
      IF (NU.EQ.0) GO TO 90
      IF (MODE.LT.1 .OR. MODE.GT.2) GO TO 110
      GO TO (10, 20), MODE
   10 IF (ABS(SARG).GT.1.0) GO TO 120
      IF (ABS(SARG).EQ.1.0) GO TO 90
      X = SARG
      SX = SQRT((1.0+ABS(X))*((0.5-ABS(X))+0.5))
      TX = X/SX
      ISIG = LOG10(2.0*NU*(5.0+TX**2))
      GO TO 30
   20 IF (ABS(SARG).GT.4.0*ATAN(1.0)) GO TO 120
      IF (SARG.EQ.0.0) GO TO 90
      X = COS(SARG)
      SX = ABS(SIN(SARG))
      TX = X/SX
      ISIG = LOG10(2.0*NU*(5.0+ABS(SARG*TX)))
C
C        BEGIN CALCULATION
C
   30 MU = MU2
      I = MU2 - MU1 + 1
C
C        IF MU.GT.NU, NORMALIZED LEGENDRE(NU,MU,X)=0.
C
   40 IF (MU.LE.NU) GO TO 50
      SPN(I) = 0.0
      IPN(I) = 0
      I = I - 1
      MU = MU - 1
      IF (I .GT. 0) GO TO 40
      ISIG = 0
      GO TO 160
   50 MU = NU
C
C        P1 = 0. = NORMALIZED LEGENDRE(NU,NU+1,X)
C
      P1 = 0.0
      IP1 = 0
C
C        CALCULATE P2 = NORMALIZED LEGENDRE(NU,NU,X)
C
      P2 = 1.0
      IP2 = 0
      P3 = 0.5
      RK = 2.0
      DO 60 J=1,NU
        P3 = ((RK+1.0)/RK)*P3
        P2 = P2*SX
        CALL XADJ(P2, IP2,IERROR)
        IF (IERROR.NE.0) RETURN
        RK = RK + 2.0
   60 CONTINUE
      P2 = P2*SQRT(P3)
      CALL XADJ(P2, IP2,IERROR)
      IF (IERROR.NE.0) RETURN
      S = 2.0*TX
      T = 1.0/NU
      IF (MU2.LT.NU) GO TO 70
      SPN(I) = P2
      IPN(I) = IP2
      I = I - 1
      IF (I .EQ. 0) GO TO 140
C
C        RECURRENCE PROCESS
C
   70 P = MU*T
      C1 = 1.0/SQRT((1.0-P+T)*(1.0+P))
      C2 = S*P*C1*P2
      C1 = -SQRT((1.0+P+T)*(1.0-P))*C1*P1
      CALL XADD(C2, IP2, C1, IP1, P, IP,IERROR)
      IF (IERROR.NE.0) RETURN
      MU = MU - 1
      IF (MU.GT.MU2) GO TO 80
C
C        STORE IN ARRAY SPN FOR RETURN TO CALLING ROUTINE.
C
      SPN(I) = P
      IPN(I) = IP
      I = I - 1
      IF (I .EQ. 0) GO TO 140
   80 P1 = P2
      IP1 = IP2
      P2 = P
      IP2 = IP
      IF (MU.LE.MU1) GO TO 140
      GO TO 70
C
C        SPECIAL CASE WHEN X=-1 OR +1, OR NU=0.
C
   90 K = MU2 - MU1 + 1
      DO 100 I=1,K
        SPN(I) = 0.0
        IPN(I) = 0
  100 CONTINUE
      ISIG = 0
      IF (MU1.GT.0) GO TO 160
      ISIG = 1
      SPN(1) = SQRT(NU+0.5)
      IPN(1) = 0
      IF (MOD(NU,2).EQ.0) GO TO 160
      IF (MODE.EQ.1 .AND. SARG.EQ.1.0) GO TO 160
      IF (MODE.EQ.2) GO TO 160
      SPN(1) = -SPN(1)
      GO TO 160
C
C          ERROR PRINTOUTS AND TERMINATION.
C
  110 CALL XERMSG ('SLATEC', 'XNRMP', 'NU, MU1, MU2 or MODE not valid',
     +             112, 1)
      IERROR=112
      RETURN
  120 CALL XERMSG ('SLATEC', 'XNRMP', 'SARG out of range', 113, 1)
      IERROR=113
      RETURN
C
C        RETURN TO CALLING PROGRAM
C
  140 K = MU2 - MU1 + 1
      DO 150 I=1,K
        CALL XRED(SPN(I),IPN(I),IERROR)
        IF (IERROR.NE.0) RETURN
  150 CONTINUE
  160 RETURN
      END