SUBROUTINE XSNRMP(NU,MU1,MU2,SARG,MODE,SPN,IPN,ISIG)
C***BEGIN PROLOGUE  XSNRMP
C***DATE WRITTEN   820712   (YYMMDD)
C***REVISION DATE  871110   (YYMMDD)
C***CATEGORY NO.  C3a2,C9
C***KEYWORDS  LEGENDRE FUNCTIONS
C***AUTHOR  LOZIER, DANIEL W. (NATIONAL BUREAU OF STANDARDS)
C           SMITH, JOHN M. (NBS AND GEORGE MASON UNIVERSITY)
C***PURPOSE  TO COMPUTE THE NORMALIZED LEGENDRE POLYNOMIAL
C***DESCRIPTION
C
C        SUBROUTINE TO CALCULATE NORMALIZED LEGENDRE POLYNOMIALS
C        (XDNRMP is double-precision version)
C        XSNRMP calculates normalized Legendre polynomials of varying
C        order and fixed argument and degree. The order MU and degree
C        NU are nonegative 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 XSSET.
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 XSNRMP 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 XSNRMP 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 event
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        XSCON. 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 XSCON(SPN(I), IPN(I))
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 XSCON, (SPN(I),
C        IPN(I)) for I = J+1, J+2, ... cannot be used in subsequent
C        extended-range computations.
C
C***REFERENCES  (SEE DESCRIPTION ABOVE)
C***ROUTINES CALLED  XERROR, XSADD, XSADJ, XSRED, XSSET
C***END PROLOGUE  XSNRMP