1*
2*  original code from the Slatec library
3*
4*  slight modifications by Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr> :
5*
6*     1/ some (minor modifications) so that the "enter" is
7*        now X and not THETA (X=cos(THETA)). This leads to
8*        better accuracy for x near 0 and seems more
9*        natural but may be there is some drawback ?
10*     2/ remove parts which send warning messages to the
11*        Slatec XERMSG routine (nevertheless all the errors
12*        flags are communicated throw IERROR).
13*        Normaly the scilab interface verify the validity of the
14*        input arguments but the verifications in this code are
15*        still here.
16*     3/ substitute calls to I1MACH by calls to dlamch
17*        (Scilab uses dlamch and not I1MACH to get machine
18*        parameter so it seems more logical).
19*
20*    IERROR values :
21*          210 :  DNU1, NUDIFF, MU1, MU2, or ID not valid
22*          211 :  X out of range (must be in [0,1)
23*          201, 202, 203, 204 : invalid input was provided to DXSET
24*                       (should not occurred in IEEE floating point)
25*          205, 206 : internal consistency error occurred in DXSET
26*                     (probably due to a software malfunction in the
27*                      library routine I1MACH) Should not occurred
28*                      in IEEE floating point, if dlamch works well.
29*          207 : an overflow or underflow of an extended-range number
30*                was detected in DXADJ.
31*          208 : an error which may occur in DXC210 but this one is not
32*                call from DXLEGF (don't know why it is given below).
33*
34*     Normally on the return to scilab, only 207 may be present.
35
36
37*DECK DXLEGF
38      SUBROUTINE DXLEGF(DNU1, NUDIFF, MU1, MU2, X, ID, PQA, IPQA,
39     1   IERROR)
40C***BEGIN PROLOGUE  DXLEGF
41C***PURPOSE  Compute normalized Legendre polynomials and associated
42C            Legendre functions.
43C***LIBRARY   SLATEC
44C***CATEGORY  C3A2, C9
45C***TYPE      DOUBLE PRECISION (XLEGF-S, DXLEGF-D)
46C***KEYWORDS  LEGENDRE FUNCTIONS
47C***AUTHOR  Smith, John M., (NBS and George Mason University)
48C***DESCRIPTION
49C
50C   DXLEGF: Extended-range Double-precision Legendre Functions
51C
52C   A feature of the DXLEGF subroutine for Legendre functions is
53C the use of extended-range arithmetic, a software extension of
54C ordinary floating-point arithmetic that greatly increases the
55C exponent range of the representable numbers. This avoids the
56C need for scaling the solutions to lie within the exponent range
57C of the most restrictive manufacturer's hardware. The increased
58C exponent range is achieved by allocating an integer storage
59C location together with each floating-point storage location.
60C
61C   The interpretation of the pair (X,I) where X is floating-point
62C and I is integer is X*(IR**I) where IR is the internal radix of
63C the computer arithmetic.
64C
65C   This subroutine computes one of the following vectors:
66C
67C 1. Legendre function of the first kind of negative order, either
68C    a. P(-MU1,NU,X), P(-MU1-1,NU,X), ..., P(-MU2,NU,X) or
69C    b. P(-MU,NU1,X), P(-MU,NU1+1,X), ..., P(-MU,NU2,X)
70C 2. Legendre function of the second kind, either
71C    a. Q(MU1,NU,X), Q(MU1+1,NU,X), ..., Q(MU2,NU,X) or
72C    b. Q(MU,NU1,X), Q(MU,NU1+1,X), ..., Q(MU,NU2,X)
73C 3. Legendre function of the first kind of positive order, either
74C    a. P(MU1,NU,X), P(MU1+1,NU,X), ..., P(MU2,NU,X) or
75C    b. P(MU,NU1,X), P(MU,NU1+1,X), ..., P(MU,NU2,X)
76C 4. Normalized Legendre polynomials, either
77C    a. PN(MU1,NU,X), PN(MU1+1,NU,X), ..., PN(MU2,NU,X) or
78C    b. PN(MU,NU1,X), PN(MU,NU1+1,X), ..., PN(MU,NU2,X)
79C
80C where X = COS(THETA).
81C
82C   The input values to DXLEGF are DNU1, NUDIFF, MU1, MU2, THETA (now X),
83C and ID. These must satisfy
84C
85C    DNU1 is DOUBLE PRECISION and greater than or equal to -0.5;
86C    NUDIFF is INTEGER and non-negative;
87C    MU1 is INTEGER and non-negative;
88C    MU2 is INTEGER and greater than or equal to MU1;
89
90C    THETA is DOUBLE PRECISION and in the half-open interval (0,PI/2];
91*    modification :  X is given (and not THETA) X must be in [0,1)
92
93C    ID is INTEGER and equal to 1, 2, 3 or 4;
94C
95C and  additionally either NUDIFF = 0 or MU2 = MU1.
96C
97C   If ID=1 and NUDIFF=0, a vector of type 1a above is computed
98C with NU=DNU1.
99C
100C   If ID=1 and MU1=MU2, a vector of type 1b above is computed
101C with NU1=DNU1, NU2=DNU1+NUDIFF and MU=MU1.
102C
103C   If ID=2 and NUDIFF=0, a vector of type 2a above is computed
104C with NU=DNU1.
105C
106C   If ID=2 and MU1=MU2, a vector of type 2b above is computed
107C with NU1=DNU1, NU2=DNU1+NUDIFF and MU=MU1.
108C
109C   If ID=3 and NUDIFF=0, a vector of type 3a above is computed
110C with NU=DNU1.
111C
112C   If ID=3 and MU1=MU2, a vector of type 3b above is computed
113C with NU1=DNU1, NU2=DNU1+NUDIFF and MU=MU1.
114C
115C   If ID=4 and NUDIFF=0, a vector of type 4a above is computed
116C with NU=DNU1.
117C
118C   If ID=4 and MU1=MU2, a vector of type 4b above is computed
119C with NU1=DNU1, NU2=DNU1+NUDIFF and MU=MU1.
120C
121C   In each case the vector of computed Legendre function values
122C is returned in the extended-range vector (PQA(I),IPQA(I)). The
123C length of this vector is either MU2-MU1+1 or NUDIFF+1.
124C
125C   Where possible, DXLEGF returns IPQA(I) as zero. In this case the
126C value of the Legendre function is contained entirely in PQA(I),
127C so it can be used in subsequent computations without further
128C consideration of extended-range arithmetic. If IPQA(I) is nonzero,
129C then the value of the Legendre function is not representable in
130C floating-point because of underflow or overflow. The program that
131C calls DXLEGF must test IPQA(I) to ensure correct usage.
132C
133C   IERROR is an error indicator. If no errors are detected, IERROR=0
134C when control returns to the calling routine. If an error is detected,
135C IERROR is returned as nonzero. The calling routine must check the
136C value of IERROR.
137C
138C   If IERROR=210 or 211, invalid input was provided to DXLEGF.
139C   If IERROR=201,202,203, or 204, invalid input was provided to DXSET.
140C   If IERROR=205 or 206, an internal consistency error occurred in
141C DXSET (probably due to a software malfunction in the library routine
142C I1MACH).
143C   If IERROR=207, an overflow or underflow of an extended-range number
144C was detected in DXADJ.
145C   If IERROR=208, an overflow or underflow of an extended-range number
146C was detected in DXC210.
147C
148C***SEE ALSO  DXSET
149C***REFERENCES  Olver and Smith, Associated Legendre Functions on the
150C                 Cut, J Comp Phys, v 51, n 3, Sept 1983, pp 502--518.
151C               Smith, Olver and Lozier, Extended-Range Arithmetic and
152C                 Normalized Legendre Polynomials, ACM Trans on Math
153C                 Softw, v 7, n 1, March 1981, pp 93--105.
154C***ROUTINES CALLED  DXPMU, DXPMUP, DXPNRM, DXPQNU, DXQMU, DXQNU, DXRED,
155C                    DXSET, XERMSG
156C***REVISION HISTORY  (YYMMDD)
157C   820728  DATE WRITTEN
158C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
159C   901019  Revisions to prologue.  (DWL and WRB)
160C   901106  Changed all specific intrinsics to generic.  (WRB)
161C           Corrected order of sections in prologue and added TYPE
162C           section.  (WRB)
163C           CALLs to XERROR changed to CALLs to XERMSG.  (WRB)
164C   920127  Revised PURPOSE section of prologue.  (DWL)
165C***END PROLOGUE  DXLEGF
166      DOUBLE PRECISION PQA,DNU1,DNU2,SX,X,PI2
167      DIMENSION PQA(*),IPQA(*)
168C
169C***FIRST EXECUTABLE STATEMENT  DXLEGF
170      IERROR=0
171      CALL DXSET (0, 0, 0.0D0, 0,IERROR)
172      IF (IERROR.NE.0) RETURN
173      PI2=2.D0*ATAN(1.D0)
174C
175C        ZERO OUTPUT ARRAYS
176C
177      L=(MU2-MU1)+NUDIFF+1
178      DO 290 I=1,L
179      PQA(I)=0.D0
180  290 IPQA(I)=0
181C
182C        CHECK FOR VALID INPUT VALUES
183C
184*** normally all these are verified by the scilab interface
185      IF(NUDIFF.LT.0) GO TO 400
186      IF(DNU1.LT.-.5D0) GO TO 400
187      IF(MU2.LT.MU1) GO TO 400
188      IF(MU1.LT.0) GO TO 400
189*      IF(THETA.LE.0.D0.OR.THETA.GT.PI2) GO TO 420
190      IF(X.LT.0.D0.OR.X.GT.1.d0) GO TO 420
191      IF(ID.LT.1.OR.ID.GT.4) GO TO 400
192      IF((MU1.NE.MU2).AND.(NUDIFF.GT.0)) GO TO 400
193C
194C        IF DNU1 IS NOT AN INTEGER, NORMALIZED P(MU,DNU,X)
195C        CANNOT BE CALCULATED.  IF DNU1 IS AN INTEGER AND
196C        MU1.GT.DNU2 THEN ALL VALUES OF P(+MU,DNU,X) AND
197C        NORMALIZED P(MU,NU,X) WILL BE ZERO.
198C
199      DNU2=DNU1+NUDIFF
200      IF((ID.EQ.3).AND.(MOD(DNU1,1.D0).NE.0.D0)) GO TO 295
201      IF((ID.EQ.4).AND.(MOD(DNU1,1.D0).NE.0.D0)) GO TO 400
202      IF((ID.EQ.3.OR.ID.EQ.4).AND.MU1.GT.DNU2) RETURN
203  295 CONTINUE
204C
205*      X=COS(THETA)
206*      SX=1.D0/SIN(THETA)
207      SX=1.D0/SQRT((1.d0-X)*(1.d0+X))
208      IF(ID.EQ.2) GO TO 300
209      IF(MU2-MU1.LE.0) GO TO 360
210C
211C        FIXED NU, VARIABLE MU
212C        CALL DXPMU TO CALCULATE P(-MU1,NU,X),....,P(-MU2,NU,X)
213C
214      CALL DXPMU(DNU1,DNU2,MU1,MU2,X,SX,ID,PQA,IPQA,IERROR)
215      IF (IERROR.NE.0) RETURN
216      GO TO 380
217C
218  300 IF(MU2.EQ.MU1) GO TO 320
219C
220C        FIXED NU, VARIABLE MU
221C        CALL DXQMU TO CALCULATE Q(MU1,NU,X),....,Q(MU2,NU,X)
222C
223      CALL DXQMU(DNU1,DNU2,MU1,MU2,X,SX,ID,PQA,IPQA,IERROR)
224      IF (IERROR.NE.0) RETURN
225      GO TO 390
226C
227C        FIXED MU, VARIABLE NU
228C        CALL DXQNU TO CALCULATE Q(MU,DNU1,X),....,Q(MU,DNU2,X)
229C
230  320 CALL DXQNU(DNU1,DNU2,MU1,X,SX,ID,PQA,IPQA,IERROR)
231      IF (IERROR.NE.0) RETURN
232      GO TO 390
233C
234C        FIXED MU, VARIABLE NU
235C        CALL DXPQNU TO CALCULATE P(-MU,DNU1,X),....,P(-MU,DNU2,X)
236C
237  360 CALL DXPQNU(DNU1,DNU2,MU1,X,SX,ID,PQA,IPQA,IERROR)
238      IF (IERROR.NE.0) RETURN
239C
240C        IF ID = 3, TRANSFORM P(-MU,NU,X) VECTOR INTO
241C        P(MU,NU,X) VECTOR.
242C
243  380 IF(ID.EQ.3) CALL DXPMUP(DNU1,DNU2,MU1,MU2,PQA,IPQA,IERROR)
244      IF (IERROR.NE.0) RETURN
245C
246C        IF ID = 4, TRANSFORM P(-MU,NU,X) VECTOR INTO
247C        NORMALIZED P(MU,NU,X) VECTOR.
248C
249      IF(ID.EQ.4) CALL DXPNRM(DNU1,DNU2,MU1,MU2,PQA,IPQA,IERROR)
250      IF (IERROR.NE.0) RETURN
251C
252C        PLACE RESULTS IN REDUCED FORM IF POSSIBLE
253C        AND RETURN TO MAIN PROGRAM.
254C
255  390 DO 395 I=1,L
256      CALL DXRED(PQA(I),IPQA(I),IERROR)
257      IF (IERROR.NE.0) RETURN
258  395 CONTINUE
259      RETURN
260C
261C        *****     ERROR TERMINATION     *****
262C
263*  400 CALL XERMSG ('SLATEC', 'DXLEGF',
264*     +             'DNU1, NUDIFF, MU1, MU2, or ID not valid', 210, 1)
265 400  continue
266      IERROR=210
267      RETURN
268*  420 CALL XERMSG ('SLATEC', 'DXLEGF', 'THETA out of range', 211, 1)
269 420  continue
270      IERROR=211
271      RETURN
272      END
273*DECK DXPMU
274      SUBROUTINE DXPMU (NU1, NU2, MU1, MU2, X, SX, ID, PQA, IPQA,
275     1   IERROR)
276C***BEGIN PROLOGUE  DXPMU
277C***SUBSIDIARY
278C***PURPOSE  To compute the values of Legendre functions for DXLEGF.
279C            Method: backward mu-wise recurrence for P(-MU,NU,X) for
280C            fixed nu to obtain P(-MU2,NU1,X), P(-(MU2-1),NU1,X), ...,
281C            P(-MU1,NU1,X) and store in ascending mu order.
282C***LIBRARY   SLATEC
283C***CATEGORY  C3A2, C9
284C***TYPE      DOUBLE PRECISION (XPMU-S, DXPMU-D)
285C***KEYWORDS  LEGENDRE FUNCTIONS
286C***AUTHOR  Smith, John M., (NBS and George Mason University)
287C***ROUTINES CALLED  DXADD, DXADJ, DXPQNU
288C***REVISION HISTORY  (YYMMDD)
289C   820728  DATE WRITTEN
290C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
291C   901019  Revisions to prologue.  (DWL and WRB)
292C   901106  Changed all specific intrinsics to generic.  (WRB)
293C           Corrected order of sections in prologue and added TYPE
294C           section.  (WRB)
295C   920127  Revised PURPOSE section of prologue.  (DWL)
296C***END PROLOGUE  DXPMU
297      DOUBLE PRECISION PQA,NU1,NU2,P0,X,SX,X1,X2
298      DIMENSION PQA(*),IPQA(*)
299C
300C        CALL DXPQNU TO OBTAIN P(-MU2,NU,X)
301C
302C***FIRST EXECUTABLE STATEMENT  DXPMU
303      IERROR=0
304      CALL DXPQNU(NU1,NU2,MU2,X,SX,ID,PQA,IPQA,IERROR)
305      IF (IERROR.NE.0) RETURN
306      P0=PQA(1)
307      IP0=IPQA(1)
308      MU=MU2-1
309C
310C        CALL DXPQNU TO OBTAIN P(-MU2-1,NU,X)
311C
312      CALL DXPQNU(NU1,NU2,MU,X,SX,ID,PQA,IPQA,IERROR)
313      IF (IERROR.NE.0) RETURN
314      N=MU2-MU1+1
315      PQA(N)=P0
316      IPQA(N)=IP0
317      IF(N.EQ.1) GO TO 300
318      PQA(N-1)=PQA(1)
319      IPQA(N-1)=IPQA(1)
320      IF(N.EQ.2) GO TO 300
321      J=N-2
322  290 CONTINUE
323C
324C        BACKWARD RECURRENCE IN MU TO OBTAIN
325C              P(-MU2,NU1,X),P(-(MU2-1),NU1,X),....P(-MU1,NU1,X)
326C              USING
327C              (NU-MU)*(NU+MU+1.)*P(-(MU+1),NU,X)=
328C                2.*MU*X*SQRT((1./(1.-X**2))*P(-MU,NU,X)-P(-(MU-1),NU,X)
329C
330      X1=2.D0*MU*X*SX*PQA(J+1)
331      X2=-(NU1-MU)*(NU1+MU+1.D0)*PQA(J+2)
332      CALL DXADD(X1,IPQA(J+1),X2,IPQA(J+2),PQA(J),IPQA(J),IERROR)
333      IF (IERROR.NE.0) RETURN
334      CALL DXADJ(PQA(J),IPQA(J),IERROR)
335      IF (IERROR.NE.0) RETURN
336      IF(J.EQ.1) GO TO 300
337      J=J-1
338      MU=MU-1
339      GO TO 290
340  300 RETURN
341      END
342*DECK DXPMUP
343      SUBROUTINE DXPMUP (NU1, NU2, MU1, MU2, PQA, IPQA, IERROR)
344C***BEGIN PROLOGUE  DXPMUP
345C***SUBSIDIARY
346C***PURPOSE  To compute the values of Legendre functions for DXLEGF.
347C            This subroutine transforms an array of Legendre functions
348C            of the first kind of negative order stored in array PQA
349C            into Legendre functions of the first kind of positive
350C            order stored in array PQA. The original array is destroyed.
351C***LIBRARY   SLATEC
352C***CATEGORY  C3A2, C9
353C***TYPE      DOUBLE PRECISION (XPMUP-S, DXPMUP-D)
354C***KEYWORDS  LEGENDRE FUNCTIONS
355C***AUTHOR  Smith, John M., (NBS and George Mason University)
356C***ROUTINES CALLED  DXADJ
357C***REVISION HISTORY  (YYMMDD)
358C   820728  DATE WRITTEN
359C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
360C   901019  Revisions to prologue.  (DWL and WRB)
361C   901106  Changed all specific intrinsics to generic.  (WRB)
362C           Corrected order of sections in prologue and added TYPE
363C           section.  (WRB)
364C   920127  Revised PURPOSE section of prologue.  (DWL)
365C***END PROLOGUE  DXPMUP
366      DOUBLE PRECISION DMU,NU,NU1,NU2,PQA,PROD
367      DIMENSION PQA(*),IPQA(*)
368C***FIRST EXECUTABLE STATEMENT  DXPMUP
369      IERROR=0
370      NU=NU1
371      MU=MU1
372      DMU=MU
373      N=INT(NU2-NU1+.1D0)+(MU2-MU1)+1
374      J=1
375*      IF(MOD(REAL(NU),1.).NE.0.) GO TO 210
376      IF(MOD(NU,1.d0).NE.0.d0) GO TO 210
377  200 IF(DMU.LT.NU+1.D0) GO TO 210
378      PQA(J)=0.D0
379      IPQA(J)=0
380      J=J+1
381      IF(J.GT.N) RETURN
382C        INCREMENT EITHER MU OR NU AS APPROPRIATE.
383      IF(NU2-NU1.GT..5D0) NU=NU+1.D0
384      IF(MU2.GT.MU1) MU=MU+1
385      GO TO 200
386C
387C        TRANSFORM P(-MU,NU,X) TO P(MU,NU,X) USING
388C        P(MU,NU,X)=(NU-MU+1)*(NU-MU+2)*...*(NU+MU)*P(-MU,NU,X)*(-1)**MU
389C
390  210 PROD=1.D0
391      IPROD=0
392      K=2*MU
393      IF(K.EQ.0) GO TO 222
394      DO 220 L=1,K
395      PROD=PROD*(DMU-NU-L)
396  220 CALL DXADJ(PROD,IPROD,IERROR)
397      IF (IERROR.NE.0) RETURN
398  222 CONTINUE
399      DO 240 I=J,N
400      IF(MU.EQ.0) GO TO 225
401      PQA(I)=PQA(I)*PROD*(-1)**MU
402      IPQA(I)=IPQA(I)+IPROD
403      CALL DXADJ(PQA(I),IPQA(I),IERROR)
404      IF (IERROR.NE.0) RETURN
405  225 IF(NU2-NU1.GT..5D0) GO TO 230
406      PROD=(DMU-NU)*PROD*(-DMU-NU-1.D0)
407      CALL DXADJ(PROD,IPROD,IERROR)
408      IF (IERROR.NE.0) RETURN
409      MU=MU+1
410      DMU=DMU+1.D0
411      GO TO 240
412  230 PROD=PROD*(-DMU-NU-1.D0)/(DMU-NU-1.D0)
413      CALL DXADJ(PROD,IPROD,IERROR)
414      IF (IERROR.NE.0) RETURN
415      NU=NU+1.D0
416  240 CONTINUE
417      RETURN
418      END
419*DECK DXPNRM
420      SUBROUTINE DXPNRM (NU1, NU2, MU1, MU2, PQA, IPQA, IERROR)
421C***BEGIN PROLOGUE  DXPNRM
422C***SUBSIDIARY
423C***PURPOSE  To compute the values of Legendre functions for DXLEGF.
424C            This subroutine transforms an array of Legendre functions
425C            of the first kind of negative order stored in array PQA
426C            into normalized Legendre polynomials stored in array PQA.
427C            The original array is destroyed.
428C***LIBRARY   SLATEC
429C***CATEGORY  C3A2, C9
430C***TYPE      DOUBLE PRECISION (XPNRM-S, DXPNRM-D)
431C***KEYWORDS  LEGENDRE FUNCTIONS
432C***AUTHOR  Smith, John M., (NBS and George Mason University)
433C***ROUTINES CALLED  DXADJ
434C***REVISION HISTORY  (YYMMDD)
435C   820728  DATE WRITTEN
436C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
437C   901019  Revisions to prologue.  (DWL and WRB)
438C   901106  Changed all specific intrinsics to generic.  (WRB)
439C           Corrected order of sections in prologue and added TYPE
440C           section.  (WRB)
441C   920127  Revised PURPOSE section of prologue.  (DWL)
442C***END PROLOGUE  DXPNRM
443      DOUBLE PRECISION C1,DMU,NU,NU1,NU2,PQA,PROD
444      DIMENSION PQA(*),IPQA(*)
445C***FIRST EXECUTABLE STATEMENT  DXPNRM
446      IERROR=0
447      L=(MU2-MU1)+(NU2-NU1+1.5D0)
448      MU=MU1
449      DMU=MU1
450      NU=NU1
451C
452C         IF MU .GT.NU, NORM P =0.
453C
454      J=1
455  500 IF(DMU.LE.NU) GO TO 505
456      PQA(J)=0.D0
457      IPQA(J)=0
458      J=J+1
459      IF(J.GT.L) RETURN
460C
461C        INCREMENT EITHER MU OR NU AS APPROPRIATE.
462C
463      IF(MU2.GT.MU1) DMU=DMU+1.D0
464      IF(NU2-NU1.GT..5D0) NU=NU+1.D0
465      GO TO 500
466C
467C         TRANSFORM P(-MU,NU,X) INTO NORMALIZED P(MU,NU,X) USING
468C              NORM P(MU,NU,X)=
469C                 SQRT((NU+.5)*FACTORIAL(NU+MU)/FACTORIAL(NU-MU))
470C                              *P(-MU,NU,X)
471C
472  505 PROD=1.D0
473      IPROD=0
474      K=2*MU
475      IF(K.LE.0) GO TO 520
476      DO 510 I=1,K
477      PROD=PROD*SQRT(NU+DMU+1.D0-I)
478  510 CALL DXADJ(PROD,IPROD,IERROR)
479      IF (IERROR.NE.0) RETURN
480  520 DO 540 I=J,L
481      C1=PROD*SQRT(NU+.5D0)
482      PQA(I)=PQA(I)*C1
483      IPQA(I)=IPQA(I)+IPROD
484      CALL DXADJ(PQA(I),IPQA(I),IERROR)
485      IF (IERROR.NE.0) RETURN
486      IF(NU2-NU1.GT..5D0) GO TO 530
487      IF(DMU.GE.NU) GO TO 525
488      PROD=SQRT(NU+DMU+1.D0)*PROD
489      IF(NU.GT.DMU) PROD=PROD*SQRT(NU-DMU)
490      CALL DXADJ(PROD,IPROD,IERROR)
491      IF (IERROR.NE.0) RETURN
492      MU=MU+1
493      DMU=DMU+1.D0
494      GO TO 540
495  525 PROD=0.D0
496      IPROD=0
497      MU=MU+1
498      DMU=DMU+1.D0
499      GO TO 540
500  530 PROD=SQRT(NU+DMU+1.D0)*PROD
501      IF(NU.NE.DMU-1.D0) PROD=PROD/SQRT(NU-DMU+1.D0)
502      CALL DXADJ(PROD,IPROD,IERROR)
503      IF (IERROR.NE.0) RETURN
504      NU=NU+1.D0
505  540 CONTINUE
506      RETURN
507      END
508*DECK DXPQNU
509      SUBROUTINE DXPQNU (NU1, NU2, MU, X, SX, ID, PQA, IPQA, IERROR)
510C***BEGIN PROLOGUE  DXPQNU
511C***SUBSIDIARY
512C***PURPOSE  To compute the values of Legendre functions for DXLEGF.
513C            This subroutine calculates initial values of P or Q using
514C            power series, then performs forward nu-wise recurrence to
515C            obtain P(-MU,NU,X), Q(0,NU,X), or Q(1,NU,X). The nu-wise
516C            recurrence is stable for P for all mu and for Q for mu=0,1.
517C***LIBRARY   SLATEC
518C***CATEGORY  C3A2, C9
519C***TYPE      DOUBLE PRECISION (XPQNU-S, DXPQNU-D)
520C***KEYWORDS  LEGENDRE FUNCTIONS
521C***AUTHOR  Smith, John M., (NBS and George Mason University)
522C***ROUTINES CALLED  DXADD, DXADJ, DXPSI
523C***COMMON BLOCKS    DXBLK1
524C***REVISION HISTORY  (YYMMDD)
525C   820728  DATE WRITTEN
526C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
527C   901019  Revisions to prologue.  (DWL and WRB)
528C   901106  Changed all specific intrinsics to generic.  (WRB)
529C           Corrected order of sections in prologue and added TYPE
530C           section.  (WRB)
531C   920127  Revised PURPOSE section of prologue.  (DWL)
532C***END PROLOGUE  DXPQNU
533      DOUBLE PRECISION A,NU,NU1,NU2,PQ,PQA,DXPSI,R,W,X,X1,X2,SX,XS,
534     1 Y,Z
535      DOUBLE PRECISION DI,DMU,PQ1,PQ2,FACTMU,FLOK
536      DIMENSION PQA(*),IPQA(*)
537      COMMON /DXBLK1/ NBITSF
538      SAVE /DXBLK1/
539C
540C        J0, IPSIK, AND IPSIX ARE INITIALIZED IN THIS SUBROUTINE.
541C        J0 IS THE NUMBER OF TERMS USED IN SERIES EXPANSION
542C        IN SUBROUTINE DXPQNU.
543C        IPSIK, IPSIX ARE VALUES OF K AND X RESPECTIVELY
544C        USED IN THE CALCULATION OF THE DXPSI FUNCTION.
545C
546C***FIRST EXECUTABLE STATEMENT  DXPQNU
547      IERROR=0
548      J0=NBITSF
549      IPSIK=1+(NBITSF/10)
550      IPSIX=5*IPSIK
551      IPQ=0
552C        FIND NU IN INTERVAL [-.5,.5) IF ID=2  ( CALCULATION OF Q )
553      NU=MOD(NU1,1.D0)
554      IF(NU.GE..5D0) NU=NU-1.D0
555C        FIND NU IN INTERVAL (-1.5,-.5] IF ID=1,3, OR 4  ( CALC. OF P )
556      IF(ID.NE.2.AND.NU.GT.-.5D0) NU=NU-1.D0
557C        CALCULATE MU FACTORIAL
558      K=MU
559      DMU=MU
560      IF(MU.LE.0) GO TO 60
561      FACTMU=1.D0
562      IF=0
563      DO 50 I=1,K
564      FACTMU=FACTMU*I
565   50 CALL DXADJ(FACTMU,IF,IERROR)
566      IF (IERROR.NE.0) RETURN
567   60 IF(K.EQ.0) FACTMU=1.D0
568      IF(K.EQ.0) IF=0
569C
570C        X=COS(THETA)
571C        Y=SIN(THETA/2)**2=(1-X)/2=.5-.5*X
572C        R=TAN(THETA/2)=SQRT((1-X)/(1+X)
573C
574      Y=0.5d0*(1.d0-X)
575      R=sqrt((1.d0-X)/(1.d0+X))
576C
577C        USE ASCENDING SERIES TO CALCULATE TWO VALUES OF P OR Q
578C        FOR USE AS STARTING VALUES IN RECURRENCE RELATION.
579C
580      PQ2=0.0D0
581      DO 100 J=1,2
582      IPQ1=0
583      IF(ID.EQ.2) GO TO 80
584C
585C        SERIES FOR P ( ID = 1, 3, OR 4 )
586C        P(-MU,NU,X)=1./FACTORIAL(MU)*SQRT(((1.-X)/(1.+X))**MU)
587C                *SUM(FROM 0 TO J0-1)A(J)*(.5-.5*X)**J
588C
589      IPQ=0
590      PQ=1.D0
591      A=1.D0
592      IA=0
593      DO 65 I=2,J0
594      DI=I
595      A=A*Y*(DI-2.D0-NU)*(DI-1.D0+NU)/((DI-1.D0+DMU)*(DI-1.D0))
596      CALL DXADJ(A,IA,IERROR)
597      IF (IERROR.NE.0) RETURN
598      IF(A.EQ.0.D0) GO TO 66
599      CALL DXADD(PQ,IPQ,A,IA,PQ,IPQ,IERROR)
600      IF (IERROR.NE.0) RETURN
601   65 CONTINUE
602   66 CONTINUE
603      IF(MU.LE.0) GO TO 90
604      X2=R
605      X1=PQ
606      K=MU
607      DO 77 I=1,K
608      X1=X1*X2
609   77 CALL DXADJ(X1,IPQ,IERROR)
610      IF (IERROR.NE.0) RETURN
611      PQ=X1/FACTMU
612      IPQ=IPQ-IF
613      CALL DXADJ(PQ,IPQ,IERROR)
614      IF (IERROR.NE.0) RETURN
615      GO TO 90
616C
617C        Z=-LN(R)=.5*LN((1+X)/(1-X))
618C
619   80 Z=-LOG(R)
620      W=DXPSI(NU+1.D0,IPSIK,IPSIX)
621      XS = SX  ! pour le cas ou XS serait modifie par la suite
622*      XS=1.D0/SIN(THETA)
623C
624C        SERIES SUMMATION FOR Q ( ID = 2 )
625C        Q(0,NU,X)=SUM(FROM 0 TO J0-1)((.5*LN((1+X)/(1-X))
626C    +DXPSI(J+1,IPSIK,IPSIX)-DXPSI(NU+1,IPSIK,IPSIX)))*A(J)*(.5-.5*X)**J
627C
628C        Q(1,NU,X)=-SQRT(1./(1.-X**2))+SQRT((1-X)/(1+X))
629C             *SUM(FROM 0 T0 J0-1)(-NU*(NU+1)/2*LN((1+X)/(1-X))
630C                 +(J-NU)*(J+NU+1)/(2*(J+1))+NU*(NU+1)*
631C     (DXPSI(NU+1,IPSIK,IPSIX)-DXPSI(J+1,IPSIK,IPSIX))*A(J)*(.5-.5*X)**J
632C
633C        NOTE, IN THIS LOOP K=J+1
634C
635      PQ=0.D0
636      IPQ=0
637      IA=0
638      A=1.D0
639      DO 85 K=1,J0
640      FLOK=K
641      IF(K.EQ.1) GO TO 81
642      A=A*Y*(FLOK-2.D0-NU)*(FLOK-1.D0+NU)/((FLOK-1.D0+DMU)*(FLOK-1.D0))
643      CALL DXADJ(A,IA,IERROR)
644      IF (IERROR.NE.0) RETURN
645   81 CONTINUE
646      IF(MU.GE.1) GO TO 83
647      X1=(DXPSI(FLOK,IPSIK,IPSIX)-W+Z)*A
648      IX1=IA
649      CALL DXADD(PQ,IPQ,X1,IX1,PQ,IPQ,IERROR)
650      IF (IERROR.NE.0) RETURN
651      GO TO 85
652   83 X1=(NU*(NU+1.D0)*(Z-W+DXPSI(FLOK,IPSIK,IPSIX))+(NU-FLOK+1.D0)
653     1  *(NU+FLOK)/(2.D0*FLOK))*A
654      IX1=IA
655      CALL DXADD(PQ,IPQ,X1,IX1,PQ,IPQ,IERROR)
656      IF (IERROR.NE.0) RETURN
657   85 CONTINUE
658      IF(MU.GE.1) PQ=-R*PQ
659      IXS=0
660      IF(MU.GE.1) CALL DXADD(PQ,IPQ,-XS,IXS,PQ,IPQ,IERROR)
661      IF (IERROR.NE.0) RETURN
662      IF(J.EQ.2) MU=-MU
663      IF(J.EQ.2) DMU=-DMU
664   90 IF(J.EQ.1) PQ2=PQ
665      IF(J.EQ.1) IPQ2=IPQ
666      NU=NU+1.D0
667  100 CONTINUE
668      K=0
669      IF(NU-1.5D0.LT.NU1) GO TO 120
670      K=K+1
671      PQA(K)=PQ2
672      IPQA(K)=IPQ2
673      IF(NU.GT.NU2+.5D0) RETURN
674  120 PQ1=PQ
675      IPQ1=IPQ
676      IF(NU.LT.NU1+.5D0) GO TO 130
677      K=K+1
678      PQA(K)=PQ
679      IPQA(K)=IPQ
680      IF(NU.GT.NU2+.5D0) RETURN
681C
682C        FORWARD NU-WISE RECURRENCE FOR F(MU,NU,X) FOR FIXED MU
683C        USING
684C        (NU+MU+1)*F(MU,NU,X)=(2.*NU+1)*F(MU,NU,X)-(NU-MU)*F(MU,NU-1,X)
685C        WHERE F(MU,NU,X) MAY BE P(-MU,NU,X) OR IF MU IS REPLACED
686C        BY -MU THEN F(MU,NU,X) MAY BE Q(MU,NU,X).
687C        NOTE, IN THIS LOOP, NU=NU+1
688C
689  130 X1=(2.D0*NU-1.D0)/(NU+DMU)*X*PQ1
690      X2=(NU-1.D0-DMU)/(NU+DMU)*PQ2
691      CALL DXADD(X1,IPQ1,-X2,IPQ2,PQ,IPQ,IERROR)
692      IF (IERROR.NE.0) RETURN
693      CALL DXADJ(PQ,IPQ,IERROR)
694      IF (IERROR.NE.0) RETURN
695      NU=NU+1.D0
696      PQ2=PQ1
697      IPQ2=IPQ1
698      GO TO 120
699C
700      END
701*DECK DXPSI
702      DOUBLE PRECISION FUNCTION DXPSI (A, IPSIK, IPSIX)
703C***BEGIN PROLOGUE  DXPSI
704C***SUBSIDIARY
705C***PURPOSE  To compute values of the Psi function for DXLEGF.
706C***LIBRARY   SLATEC
707C***CATEGORY  C7C
708C***TYPE      DOUBLE PRECISION (XPSI-S, DXPSI-D)
709C***KEYWORDS  PSI FUNCTION
710C***AUTHOR  Smith, John M., (NBS and George Mason University)
711C***ROUTINES CALLED  (NONE)
712C***REVISION HISTORY  (YYMMDD)
713C   820728  DATE WRITTEN
714C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
715C   901019  Revisions to prologue.  (DWL and WRB)
716C   901106  Changed all specific intrinsics to generic.  (WRB)
717C           Corrected order of sections in prologue and added TYPE
718C           section.  (WRB)
719C   920127  Revised PURPOSE section of prologue.  (DWL)
720C***END PROLOGUE  DXPSI
721      DOUBLE PRECISION A,B,C,CNUM,CDENOM
722      DIMENSION CNUM(12),CDENOM(12)
723      SAVE CNUM, CDENOM
724C
725C        CNUM(I) AND CDENOM(I) ARE THE ( REDUCED ) NUMERATOR
726C        AND 2*I*DENOMINATOR RESPECTIVELY OF THE 2*I TH BERNOULLI
727C        NUMBER.
728C
729      DATA CNUM(1),CNUM(2),CNUM(3),CNUM(4),CNUM(5),CNUM(6),CNUM(7),
730     1CNUM(8),CNUM(9),CNUM(10),CNUM(11),CNUM(12)
731     2    / 1.D0,     -1.D0,    1.D0,     -1.D0, 1.D0,
732     3   -691.D0,  1.D0,     -3617.D0, 43867.D0, -174611.D0, 77683.D0,
733     4   -236364091.D0/
734      DATA CDENOM(1),CDENOM(2),CDENOM(3),CDENOM(4),CDENOM(5),CDENOM(6),
735     1 CDENOM(7),CDENOM(8),CDENOM(9),CDENOM(10),CDENOM(11),CDENOM(12)
736     2/12.D0,120.D0,   252.D0,   240.D0,132.D0,
737     3  32760.D0, 12.D0,  8160.D0, 14364.D0, 6600.D0, 276.D0, 65520.D0/
738C***FIRST EXECUTABLE STATEMENT  DXPSI
739      N=MAX(0,IPSIX-INT(A))
740      B=N+A
741      K1=IPSIK-1
742C
743C        SERIES EXPANSION FOR A .GT. IPSIX USING IPSIK-1 TERMS.
744C
745      C=0.D0
746      DO 12 I=1,K1
747      K=IPSIK-I
748   12 C=(C+CNUM(K)/CDENOM(K))/B**2
749      DXPSI=LOG(B)-(C+.5D0/B)
750      IF(N.EQ.0) GO TO 20
751      B=0.D0
752C
753C        RECURRENCE FOR A .LE. IPSIX.
754C
755      DO 15 M=1,N
756   15 B=B+1.D0/(N-M+A)
757      DXPSI=DXPSI-B
758   20 RETURN
759      END
760*DECK DXQMU
761      SUBROUTINE DXQMU (NU1, NU2, MU1, MU2, X, SX, ID, PQA, IPQA,
762     1   IERROR)
763C***BEGIN PROLOGUE  DXQMU
764C***SUBSIDIARY
765C***PURPOSE  To compute the values of Legendre functions for DXLEGF.
766C            Method: forward mu-wise recurrence for Q(MU,NU,X) for fixed
767C            nu to obtain Q(MU1,NU,X), Q(MU1+1,NU,X), ..., Q(MU2,NU,X).
768C***LIBRARY   SLATEC
769C***CATEGORY  C3A2, C9
770C***TYPE      DOUBLE PRECISION (XQMU-S, DXQMU-D)
771C***KEYWORDS  LEGENDRE FUNCTIONS
772C***AUTHOR  Smith, John M., (NBS and George Mason University)
773C***ROUTINES CALLED  DXADD, DXADJ, DXPQNU
774C***REVISION HISTORY  (YYMMDD)
775C   820728  DATE WRITTEN
776C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
777C   901019  Revisions to prologue.  (DWL and WRB)
778C   901106  Corrected order of sections in prologue and added TYPE
779C           section.  (WRB)
780C   920127  Revised PURPOSE section of prologue.  (DWL)
781C***END PROLOGUE  DXQMU
782      DIMENSION PQA(*),IPQA(*)
783      DOUBLE PRECISION DMU,NU,NU1,NU2,PQ,PQA,PQ1,PQ2,SX,X,X1,X2
784C***FIRST EXECUTABLE STATEMENT  DXQMU
785      IERROR=0
786      MU=0
787C
788C        CALL DXPQNU TO OBTAIN Q(0.,NU1,X)
789C
790      CALL DXPQNU(NU1,NU2,MU,X,SX,ID,PQA,IPQA,IERROR)
791      IF (IERROR.NE.0) RETURN
792      PQ2=PQA(1)
793      IPQ2=IPQA(1)
794      MU=1
795C
796C        CALL DXPQNU TO OBTAIN Q(1.,NU1,X)
797C
798      CALL DXPQNU(NU1,NU2,MU,X,SX,ID,PQA,IPQA,IERROR)
799      IF (IERROR.NE.0) RETURN
800      NU=NU1
801      K=0
802      MU=1
803      DMU=1.D0
804      PQ1=PQA(1)
805      IPQ1=IPQA(1)
806      IF(MU1.GT.0) GO TO 310
807      K=K+1
808      PQA(K)=PQ2
809      IPQA(K)=IPQ2
810      IF(MU2.LT.1) GO TO 330
811  310 IF(MU1.GT.1) GO TO 320
812      K=K+1
813      PQA(K)=PQ1
814      IPQA(K)=IPQ1
815      IF(MU2.LE.1) GO TO 330
816  320 CONTINUE
817C
818C        FORWARD RECURRENCE IN MU TO OBTAIN
819C                  Q(MU1,NU,X),Q(MU1+1,NU,X),....,Q(MU2,NU,X) USING
820C             Q(MU+1,NU,X)=-2.*MU*X*SQRT(1./(1.-X**2))*Q(MU,NU,X)
821C                               -(NU+MU)*(NU-MU+1.)*Q(MU-1,NU,X)
822C
823      X1=-2.D0*DMU*X*SX*PQ1
824      X2=(NU+DMU)*(NU-DMU+1.D0)*PQ2
825      CALL DXADD(X1,IPQ1,-X2,IPQ2,PQ,IPQ,IERROR)
826      IF (IERROR.NE.0) RETURN
827      CALL DXADJ(PQ,IPQ,IERROR)
828      IF (IERROR.NE.0) RETURN
829      PQ2=PQ1
830      IPQ2=IPQ1
831      PQ1=PQ
832      IPQ1=IPQ
833      MU=MU+1
834      DMU=DMU+1.D0
835      IF(MU.LT.MU1) GO TO 320
836      K=K+1
837      PQA(K)=PQ
838      IPQA(K)=IPQ
839      IF(MU2.GT.MU) GO TO 320
840  330 RETURN
841      END
842*DECK DXQNU
843      SUBROUTINE DXQNU (NU1, NU2, MU1, X, SX, ID, PQA, IPQA,
844     1   IERROR)
845C***BEGIN PROLOGUE  DXQNU
846C***SUBSIDIARY
847C***PURPOSE  To compute the values of Legendre functions for DXLEGF.
848C            Method: backward nu-wise recurrence for Q(MU,NU,X) for
849C            fixed mu to obtain Q(MU1,NU1,X), Q(MU1,NU1+1,X), ...,
850C            Q(MU1,NU2,X).
851C***LIBRARY   SLATEC
852C***CATEGORY  C3A2, C9
853C***TYPE      DOUBLE PRECISION (XQNU-S, DXQNU-D)
854C***KEYWORDS  LEGENDRE FUNCTIONS
855C***AUTHOR  Smith, John M., (NBS and George Mason University)
856C***ROUTINES CALLED  DXADD, DXADJ, DXPQNU
857C***REVISION HISTORY  (YYMMDD)
858C   820728  DATE WRITTEN
859C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
860C   901019  Revisions to prologue.  (DWL and WRB)
861C   901106  Corrected order of sections in prologue and added TYPE
862C           section.  (WRB)
863C   920127  Revised PURPOSE section of prologue.  (DWL)
864C***END PROLOGUE  DXQNU
865      DIMENSION PQA(*),IPQA(*)
866      DOUBLE PRECISION DMU,NU,NU1,NU2,PQ,PQA,PQ1,PQ2,SX,X,X1,X2
867      DOUBLE PRECISION PQL1,PQL2
868C***FIRST EXECUTABLE STATEMENT  DXQNU
869      IERROR=0
870      K=0
871      PQ2=0.0D0
872      IPQ2=0
873      PQL2=0.0D0
874      IPQL2=0
875      IF(MU1.EQ.1) GO TO 290
876      MU=0
877C
878C        CALL DXPQNU TO OBTAIN Q(0.,NU2,X) AND Q(0.,NU2-1,X)
879C
880      CALL DXPQNU(NU1,NU2,MU,X,SX,ID,PQA,IPQA,IERROR)
881      IF (IERROR.NE.0) RETURN
882      IF(MU1.EQ.0) RETURN
883      K=(NU2-NU1+1.5D0)
884      PQ2=PQA(K)
885      IPQ2=IPQA(K)
886      PQL2=PQA(K-1)
887      IPQL2=IPQA(K-1)
888  290 MU=1
889C
890C        CALL DXPQNU TO OBTAIN Q(1.,NU2,X) AND Q(1.,NU2-1,X)
891C
892      CALL DXPQNU(NU1,NU2,MU,X,SX,ID,PQA,IPQA,IERROR)
893      IF (IERROR.NE.0) RETURN
894      IF(MU1.EQ.1) RETURN
895      NU=NU2
896      PQ1=PQA(K)
897      IPQ1=IPQA(K)
898      PQL1=PQA(K-1)
899      IPQL1=IPQA(K-1)
900  300 MU=1
901      DMU=1.D0
902  320 CONTINUE
903C
904C        FORWARD RECURRENCE IN MU TO OBTAIN Q(MU1,NU2,X) AND
905C              Q(MU1,NU2-1,X) USING
906C              Q(MU+1,NU,X)=-2.*MU*X*SQRT(1./(1.-X**2))*Q(MU,NU,X)
907C                   -(NU+MU)*(NU-MU+1.)*Q(MU-1,NU,X)
908C
909C              FIRST FOR NU=NU2
910C
911      X1=-2.D0*DMU*X*SX*PQ1
912      X2=(NU+DMU)*(NU-DMU+1.D0)*PQ2
913      CALL DXADD(X1,IPQ1,-X2,IPQ2,PQ,IPQ,IERROR)
914      IF (IERROR.NE.0) RETURN
915      CALL DXADJ(PQ,IPQ,IERROR)
916      IF (IERROR.NE.0) RETURN
917      PQ2=PQ1
918      IPQ2=IPQ1
919      PQ1=PQ
920      IPQ1=IPQ
921      MU=MU+1
922      DMU=DMU+1.D0
923      IF(MU.LT.MU1) GO TO 320
924      PQA(K)=PQ
925      IPQA(K)=IPQ
926      IF(K.EQ.1) RETURN
927      IF(NU.LT.NU2) GO TO 340
928C
929C              THEN FOR NU=NU2-1
930C
931      NU=NU-1.D0
932      PQ2=PQL2
933      IPQ2=IPQL2
934      PQ1=PQL1
935      IPQ1=IPQL1
936      K=K-1
937      GO TO 300
938C
939C         BACKWARD RECURRENCE IN NU TO OBTAIN
940C              Q(MU1,NU1,X),Q(MU1,NU1+1,X),....,Q(MU1,NU2,X)
941C              USING
942C              (NU-MU+1.)*Q(MU,NU+1,X)=
943C                       (2.*NU+1.)*X*Q(MU,NU,X)-(NU+MU)*Q(MU,NU-1,X)
944C
945  340 PQ1=PQA(K)
946      IPQ1=IPQA(K)
947      PQ2=PQA(K+1)
948      IPQ2=IPQA(K+1)
949  350 IF(NU.LE.NU1) RETURN
950      K=K-1
951      X1=(2.D0*NU+1.D0)*X*PQ1/(NU+DMU)
952      X2=-(NU-DMU+1.D0)*PQ2/(NU+DMU)
953      CALL DXADD(X1,IPQ1,X2,IPQ2,PQ,IPQ,IERROR)
954      IF (IERROR.NE.0) RETURN
955      CALL DXADJ(PQ,IPQ,IERROR)
956      IF (IERROR.NE.0) RETURN
957      PQ2=PQ1
958      IPQ2=IPQ1
959      PQ1=PQ
960      IPQ1=IPQ
961      PQA(K)=PQ
962      IPQA(K)=IPQ
963      NU=NU-1.D0
964      GO TO 350
965      END
966*DECK DXRED
967      SUBROUTINE DXRED (X, IX, IERROR)
968C***BEGIN PROLOGUE  DXRED
969C***PURPOSE  To provide double-precision floating-point arithmetic
970C            with an extended exponent range.
971C***LIBRARY   SLATEC
972C***CATEGORY  A3D
973C***TYPE      DOUBLE PRECISION (XRED-S, DXRED-D)
974C***KEYWORDS  EXTENDED-RANGE DOUBLE-PRECISION ARITHMETIC
975C***AUTHOR  Lozier, Daniel W., (National Bureau of Standards)
976C           Smith, John M., (NBS and George Mason University)
977C***DESCRIPTION
978C     DOUBLE PRECISION X
979C     INTEGER IX
980C
981C                  IF
982C                  RADIX**(-2L) .LE. (ABS(X),IX) .LE. RADIX**(2L)
983C                  THEN DXRED TRANSFORMS (X,IX) SO THAT IX=0.
984C                  IF (X,IX) IS OUTSIDE THE ABOVE RANGE,
985C                  THEN DXRED TAKES NO ACTION.
986C                  THIS SUBROUTINE IS USEFUL IF THE
987C                  RESULTS OF EXTENDED-RANGE CALCULATIONS
988C                  ARE TO BE USED IN SUBSEQUENT ORDINARY
989C                  DOUBLE-PRECISION CALCULATIONS.
990C
991C***SEE ALSO  DXSET
992C***REFERENCES  (NONE)
993C***ROUTINES CALLED  (NONE)
994C***COMMON BLOCKS    DXBLK2
995C***REVISION HISTORY  (YYMMDD)
996C   820712  DATE WRITTEN
997C   881020  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
998C   901019  Revisions to prologue.  (DWL and WRB)
999C   901106  Changed all specific intrinsics to generic.  (WRB)
1000C           Corrected order of sections in prologue and added TYPE
1001C           section.  (WRB)
1002C   920127  Revised PURPOSE section of prologue.  (DWL)
1003C***END PROLOGUE  DXRED
1004      DOUBLE PRECISION X
1005      INTEGER IX
1006      DOUBLE PRECISION RADIX, RADIXL, RAD2L, DLG10R, XA
1007      INTEGER L, L2, KMAX
1008      COMMON /DXBLK2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX
1009      SAVE /DXBLK2/
1010C
1011C***FIRST EXECUTABLE STATEMENT  DXRED
1012      IERROR=0
1013      IF (X.EQ.0.0D0) GO TO 90
1014      XA = ABS(X)
1015      IF (IX.EQ.0) GO TO 70
1016      IXA = ABS(IX)
1017      IXA1 = IXA/L2
1018      IXA2 = MOD(IXA,L2)
1019      IF (IX.GT.0) GO TO 40
1020   10 CONTINUE
1021      IF (XA.GT.1.0D0) GO TO 20
1022      XA = XA*RAD2L
1023      IXA1 = IXA1 + 1
1024      GO TO 10
1025   20 XA = XA/RADIX**IXA2
1026      IF (IXA1.EQ.0) GO TO 70
1027      DO 30 I=1,IXA1
1028        IF (XA.LT.1.0D0) GO TO 100
1029        XA = XA/RAD2L
1030   30 CONTINUE
1031      GO TO 70
1032C
1033   40 CONTINUE
1034      IF (XA.LT.1.0D0) GO TO 50
1035      XA = XA/RAD2L
1036      IXA1 = IXA1 + 1
1037      GO TO 40
1038   50 XA = XA*RADIX**IXA2
1039      IF (IXA1.EQ.0) GO TO 70
1040      DO 60 I=1,IXA1
1041        IF (XA.GT.1.0D0) GO TO 100
1042        XA = XA*RAD2L
1043   60 CONTINUE
1044   70 IF (XA.GT.RAD2L) GO TO 100
1045      IF (XA.GT.1.0D0) GO TO 80
1046      IF (RAD2L*XA.LT.1.0D0) GO TO 100
1047   80 X = SIGN(XA,X)
1048   90 IX = 0
1049  100 RETURN
1050      END
1051*DECK DXSET
1052      SUBROUTINE DXSET (IRAD, NRADPL, DZERO, NBITS, IERROR)
1053C***BEGIN PROLOGUE  DXSET
1054C***PURPOSE  To provide double-precision floating-point arithmetic
1055C            with an extended exponent range.
1056C***LIBRARY   SLATEC
1057C***CATEGORY  A3D
1058C***TYPE      DOUBLE PRECISION (XSET-S, DXSET-D)
1059C***KEYWORDS  EXTENDED-RANGE DOUBLE-PRECISION ARITHMETIC
1060C***AUTHOR  Lozier, Daniel W., (National Bureau of Standards)
1061C           Smith, John M., (NBS and George Mason University)
1062C***DESCRIPTION
1063C
1064C   SUBROUTINE  DXSET  MUST BE CALLED PRIOR TO CALLING ANY OTHER
1065C EXTENDED-RANGE SUBROUTINE. IT CALCULATES AND STORES SEVERAL
1066C MACHINE-DEPENDENT CONSTANTS IN COMMON BLOCKS. THE USER MUST
1067C SUPPLY FOUR CONSTANTS THAT PERTAIN TO HIS PARTICULAR COMPUTER.
1068C THE CONSTANTS ARE
1069C
1070C          IRAD = THE INTERNAL BASE OF DOUBLE-PRECISION
1071C                 ARITHMETIC IN THE COMPUTER.
1072C        NRADPL = THE NUMBER OF RADIX PLACES CARRIED IN
1073C                 THE DOUBLE-PRECISION REPRESENTATION.
1074C         DZERO = THE SMALLEST OF 1/DMIN, DMAX, DMAXLN WHERE
1075C                 DMIN = THE SMALLEST POSITIVE DOUBLE-PRECISION
1076C                 NUMBER OR AN UPPER BOUND TO THIS NUMBER,
1077C                 DMAX = THE LARGEST DOUBLE-PRECISION NUMBER
1078C                 OR A LOWER BOUND TO THIS NUMBER,
1079C                 DMAXLN = THE LARGEST DOUBLE-PRECISION NUMBER
1080C                 SUCH THAT LOG10(DMAXLN) CAN BE COMPUTED BY THE
1081C                 FORTRAN SYSTEM (ON MOST SYSTEMS DMAXLN = DMAX).
1082C         NBITS = THE NUMBER OF BITS (EXCLUSIVE OF SIGN) IN
1083C                 AN INTEGER COMPUTER WORD.
1084C
1085C ALTERNATIVELY, ANY OR ALL OF THE CONSTANTS CAN BE GIVEN
1086C THE VALUE 0 (0.0D0 FOR DZERO). IF A CONSTANT IS ZERO, DXSET TRIES
1087C TO ASSIGN AN APPROPRIATE VALUE BY CALLING I1MACH
1088C (SEE P.A.FOX, A.D.HALL, N.L.SCHRYER, ALGORITHM 528 FRAMEWORK
1089C FOR A PORTABLE LIBRARY, ACM TRANSACTIONS ON MATH SOFTWARE,
1090C V.4, NO.2, JUNE 1978, 177-188).
1091C
1092C   THIS IS THE SETTING-UP SUBROUTINE FOR A PACKAGE OF SUBROUTINES
1093C THAT FACILITATE THE USE OF EXTENDED-RANGE ARITHMETIC. EXTENDED-RANGE
1094C ARITHMETIC ON A PARTICULAR COMPUTER IS DEFINED ON THE SET OF NUMBERS
1095C OF THE FORM
1096C
1097C               (X,IX) = X*RADIX**IX
1098C
1099C WHERE X IS A DOUBLE-PRECISION NUMBER CALLED THE PRINCIPAL PART,
1100C IX IS AN INTEGER CALLED THE AUXILIARY INDEX, AND RADIX IS THE
1101C INTERNAL BASE OF THE DOUBLE-PRECISION ARITHMETIC.  OBVIOUSLY,
1102C EACH REAL NUMBER IS REPRESENTABLE WITHOUT ERROR BY MORE THAN ONE
1103C EXTENDED-RANGE FORM.  CONVERSIONS BETWEEN  DIFFERENT FORMS ARE
1104C ESSENTIAL IN CARRYING OUT ARITHMETIC OPERATIONS.  WITH THE CHOICE
1105C OF RADIX WE HAVE MADE, AND THE SUBROUTINES WE HAVE WRITTEN, THESE
1106C CONVERSIONS ARE PERFORMED WITHOUT ERROR (AT LEAST ON MOST COMPUTERS).
1107C (SEE SMITH, J.M., OLVER, F.W.J., AND LOZIER, D.W., EXTENDED-RANGE
1108C ARITHMETIC AND NORMALIZED LEGENDRE POLYNOMIALS, ACM TRANSACTIONS ON
1109C MATHEMATICAL SOFTWARE, MARCH 1981).
1110C
1111C   AN EXTENDED-RANGE NUMBER  (X,IX)  IS SAID TO BE IN ADJUSTED FORM IF
1112C X AND IX ARE ZERO OR
1113C
1114C           RADIX**(-L) .LE. ABS(X) .LT. RADIX**L
1115C
1116C IS SATISFIED, WHERE L IS A COMPUTER-DEPENDENT INTEGER DEFINED IN THIS
1117C SUBROUTINE. TWO EXTENDED-RANGE NUMBERS IN ADJUSTED FORM CAN BE ADDED,
1118C SUBTRACTED, MULTIPLIED OR DIVIDED (IF THE DIVISOR IS NONZERO) WITHOUT
1119C CAUSING OVERFLOW OR UNDERFLOW IN THE PRINCIPAL PART OF THE RESULT.
1120C WITH PROPER USE OF THE EXTENDED-RANGE SUBROUTINES, THE ONLY OVERFLOW
1121C THAT CAN OCCUR IS INTEGER OVERFLOW IN THE AUXILIARY INDEX. IF THIS
1122C IS DETECTED, THE SOFTWARE CALLS XERROR (A GENERAL ERROR-HANDLING
1123C FORTRAN SUBROUTINE PACKAGE).
1124C
1125C   MULTIPLICATION AND DIVISION IS PERFORMED BY SETTING
1126C
1127C                 (X,IX)*(Y,IY) = (X*Y,IX+IY)
1128C OR
1129C                 (X,IX)/(Y,IY) = (X/Y,IX-IY).
1130C
1131C PRE-ADJUSTMENT OF THE OPERANDS IS ESSENTIAL TO AVOID
1132C OVERFLOW OR  UNDERFLOW OF THE PRINCIPAL PART. SUBROUTINE
1133C DXADJ (SEE BELOW) MAY BE CALLED TO TRANSFORM ANY EXTENDED-
1134C RANGE NUMBER INTO ADJUSTED FORM.
1135C
1136C   ADDITION AND SUBTRACTION REQUIRE THE USE OF SUBROUTINE DXADD
1137C (SEE BELOW).  THE INPUT OPERANDS NEED NOT BE IN ADJUSTED FORM.
1138C HOWEVER, THE RESULT OF ADDITION OR SUBTRACTION IS RETURNED
1139C IN ADJUSTED FORM.  THUS, FOR EXAMPLE, IF (X,IX),(Y,IY),
1140C (U,IU),  AND (V,IV) ARE IN ADJUSTED FORM, THEN
1141C
1142C                 (X,IX)*(Y,IY) + (U,IU)*(V,IV)
1143C
1144C CAN BE COMPUTED AND STORED IN ADJUSTED FORM WITH NO EXPLICIT
1145C CALLS TO DXADJ.
1146C
1147C   WHEN AN EXTENDED-RANGE NUMBER IS TO BE PRINTED, IT MUST BE
1148C CONVERTED TO AN EXTENDED-RANGE FORM WITH DECIMAL RADIX.  SUBROUTINE
1149C DXCON IS PROVIDED FOR THIS PURPOSE.
1150C
1151C   THE SUBROUTINES CONTAINED IN THIS PACKAGE ARE
1152C
1153C     SUBROUTINE DXADD
1154C USAGE
1155C                  CALL DXADD(X,IX,Y,IY,Z,IZ,IERROR)
1156C                  IF (IERROR.NE.0) RETURN
1157C DESCRIPTION
1158C                  FORMS THE EXTENDED-RANGE SUM  (Z,IZ) =
1159C                  (X,IX) + (Y,IY).  (Z,IZ) IS ADJUSTED
1160C                  BEFORE RETURNING. THE INPUT OPERANDS
1161C                  NEED NOT BE IN ADJUSTED FORM, BUT THEIR
1162C                  PRINCIPAL PARTS MUST SATISFY
1163C                  RADIX**(-2L).LE.ABS(X).LE.RADIX**(2L),
1164C                  RADIX**(-2L).LE.ABS(Y).LE.RADIX**(2L).
1165C
1166C     SUBROUTINE DXADJ
1167C USAGE
1168C                  CALL DXADJ(X,IX,IERROR)
1169C                  IF (IERROR.NE.0) RETURN
1170C DESCRIPTION
1171C                  TRANSFORMS (X,IX) SO THAT
1172C                  RADIX**(-L) .LE. ABS(X) .LT. RADIX**L.
1173C                  ON MOST COMPUTERS THIS TRANSFORMATION DOES
1174C                  NOT CHANGE THE MANTISSA OF X PROVIDED RADIX IS
1175C                  THE NUMBER BASE OF DOUBLE-PRECISION ARITHMETIC.
1176C
1177C     SUBROUTINE DXC210
1178C USAGE
1179C                  CALL DXC210(K,Z,J,IERROR)
1180C                  IF (IERROR.NE.0) RETURN
1181C DESCRIPTION
1182C                  GIVEN K THIS SUBROUTINE COMPUTES J AND Z
1183C                  SUCH THAT  RADIX**K = Z*10**J, WHERE Z IS IN
1184C                  THE RANGE 1/10 .LE. Z .LT. 1.
1185C                  THE VALUE OF Z WILL BE ACCURATE TO FULL
1186C                  DOUBLE-PRECISION PROVIDED THE NUMBER
1187C                  OF DECIMAL PLACES IN THE LARGEST
1188C                  INTEGER PLUS THE NUMBER OF DECIMAL
1189C                  PLACES CARRIED IN DOUBLE-PRECISION DOES NOT
1190C                  EXCEED 60. DXC210 IS CALLED BY SUBROUTINE
1191C                  DXCON WHEN NECESSARY. THE USER SHOULD
1192C                  NEVER NEED TO CALL DXC210 DIRECTLY.
1193C
1194C     SUBROUTINE DXCON
1195C USAGE
1196C                  CALL DXCON(X,IX,IERROR)
1197C                  IF (IERROR.NE.0) RETURN
1198C DESCRIPTION
1199C                  CONVERTS (X,IX) = X*RADIX**IX
1200C                  TO DECIMAL FORM IN PREPARATION FOR
1201C                  PRINTING, SO THAT (X,IX) = X*10**IX
1202C                  WHERE 1/10 .LE. ABS(X) .LT. 1
1203C                  IS RETURNED, EXCEPT THAT IF
1204C                  (ABS(X),IX) IS BETWEEN RADIX**(-2L)
1205C                  AND RADIX**(2L) THEN THE REDUCED
1206C                  FORM WITH IX = 0 IS RETURNED.
1207C
1208C     SUBROUTINE DXRED
1209C USAGE
1210C                  CALL DXRED(X,IX,IERROR)
1211C                  IF (IERROR.NE.0) RETURN
1212C DESCRIPTION
1213C                  IF
1214C                  RADIX**(-2L) .LE. (ABS(X),IX) .LE. RADIX**(2L)
1215C                  THEN DXRED TRANSFORMS (X,IX) SO THAT IX=0.
1216C                  IF (X,IX) IS OUTSIDE THE ABOVE RANGE,
1217C                  THEN DXRED TAKES NO ACTION.
1218C                  THIS SUBROUTINE IS USEFUL IF THE
1219C                  RESULTS OF EXTENDED-RANGE CALCULATIONS
1220C                  ARE TO BE USED IN SUBSEQUENT ORDINARY
1221C                  DOUBLE-PRECISION CALCULATIONS.
1222C
1223C***REFERENCES  Smith, Olver and Lozier, Extended-Range Arithmetic and
1224C                 Normalized Legendre Polynomials, ACM Trans on Math
1225C                 Softw, v 7, n 1, March 1981, pp 93--105.
1226C***ROUTINES CALLED  I1MACH, XERMSG
1227C***COMMON BLOCKS    DXBLK1, DXBLK2, DXBLK3
1228C***REVISION HISTORY  (YYMMDD)
1229C   820712  DATE WRITTEN
1230C   881020  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
1231C   901019  Revisions to prologue.  (DWL and WRB)
1232C   901106  Changed all specific intrinsics to generic.  (WRB)
1233C           Corrected order of sections in prologue and added TYPE
1234C           section.  (WRB)
1235C           CALLs to XERROR changed to CALLs to XERMSG.  (WRB)
1236C   920127  Revised PURPOSE section of prologue.  (DWL)
1237C***END PROLOGUE  DXSET
1238      INTEGER IRAD, NRADPL, NBITS
1239      DOUBLE PRECISION DZERO, DZEROX
1240      COMMON /DXBLK1/ NBITSF
1241      SAVE /DXBLK1/
1242      DOUBLE PRECISION RADIX, RADIXL, RAD2L, DLG10R
1243      INTEGER L, L2, KMAX
1244      COMMON /DXBLK2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX
1245      SAVE /DXBLK2/
1246      INTEGER NLG102, MLG102, LG102
1247      COMMON /DXBLK3/ NLG102, MLG102, LG102(21)
1248      SAVE /DXBLK3/
1249      INTEGER IFLAG
1250      SAVE IFLAG
1251
1252*     dlamch is used in place of i1mach :
1253      double precision dlamch
1254      external         dlamch
1255C
1256      DIMENSION LOG102(20), LGTEMP(20)
1257      SAVE LOG102
1258C
1259C   LOG102 CONTAINS THE FIRST 60 DIGITS OF LOG10(2) FOR USE IN
1260C CONVERSION OF EXTENDED-RANGE NUMBERS TO BASE 10 .
1261      DATA LOG102 /301,029,995,663,981,195,213,738,894,724,493,026,768,
1262     * 189,881,462,108,541,310,428/
1263C
1264C FOLLOWING CODING PREVENTS DXSET FROM BEING EXECUTED MORE THAN ONCE.
1265C THIS IS IMPORTANT BECAUSE SOME SUBROUTINES (SUCH AS DXNRMP AND
1266C DXLEGF) CALL DXSET TO MAKE SURE EXTENDED-RANGE ARITHMETIC HAS
1267C BEEN INITIALIZED. THE USER MAY WANT TO PRE-EMPT THIS CALL, FOR
1268C EXAMPLE WHEN I1MACH IS NOT AVAILABLE. SEE CODING BELOW.
1269      DATA IFLAG /0/
1270C***FIRST EXECUTABLE STATEMENT  DXSET
1271      IERROR=0
1272      IF (IFLAG .NE. 0) RETURN
1273      IRADX = IRAD
1274      NRDPLC = NRADPL
1275      DZEROX = DZERO
1276      IMINEX = 0
1277      IMAXEX = 0
1278      NBITSX = NBITS
1279C FOLLOWING 5 STATEMENTS SHOULD BE DELETED IF I1MACH IS
1280C NOT AVAILABLE OR NOT CONFIGURED TO RETURN THE CORRECT
1281C MACHINE-DEPENDENT VALUES.
1282*
1283* modif : use a call to dlamch in place of I1MACH
1284      IF (IRADX .EQ. 0) IRADX = int(dlamch('b'))       ! I1MACH (10)
1285      IF (NRDPLC .EQ. 0) NRDPLC = int(dlamch('n'))     ! I1MACH (14)
1286      IF (DZEROX .EQ. 0.0D0) IMINEX = int(dlamch('m')) ! I1MACH (15)
1287      IF (DZEROX .EQ. 0.0D0) IMAXEX = int(dlamch('l')) ! I1MACH (16)
1288      IF (NBITSX .EQ. 0) NBITSX = 31                   ! I1MACH (8)
1289      IF (IRADX.EQ.2) GO TO 10
1290      IF (IRADX.EQ.4) GO TO 10
1291      IF (IRADX.EQ.8) GO TO 10
1292      IF (IRADX.EQ.16) GO TO 10
1293*      CALL XERMSG ('SLATEC', 'DXSET', 'IMPROPER VALUE OF IRAD', 201, 1)
1294      IERROR=201
1295      RETURN
1296   10 CONTINUE
1297      LOG2R=0
1298      IF (IRADX.EQ.2) LOG2R = 1
1299      IF (IRADX.EQ.4) LOG2R = 2
1300      IF (IRADX.EQ.8) LOG2R = 3
1301      IF (IRADX.EQ.16) LOG2R = 4
1302      NBITSF=LOG2R*NRDPLC
1303      RADIX = IRADX
1304      DLG10R = LOG10(RADIX)
1305      IF (DZEROX .NE. 0.0D0) GO TO 14
1306      LX = MIN ((1-IMINEX)/2, (IMAXEX-1)/2)
1307      GO TO 16
1308   14 LX = 0.5D0*LOG10(DZEROX)/DLG10R
1309C RADIX**(2*L) SHOULD NOT OVERFLOW, BUT REDUCE L BY 1 FOR FURTHER
1310C PROTECTION.
1311      LX=LX-1
1312   16 L2 = 2*LX
1313      IF (LX.GE.4) GO TO 20
1314*      CALL XERMSG ('SLATEC', 'DXSET', 'IMPROPER VALUE OF DZERO', 202, 1)
1315      IERROR=202
1316      RETURN
1317   20 L = LX
1318      RADIXL = RADIX**L
1319      RAD2L = RADIXL**2
1320C    IT IS NECESSARY TO RESTRICT NBITS (OR NBITSX) TO BE LESS THAN SOME
1321C UPPER LIMIT BECAUSE OF BINARY-TO-DECIMAL CONVERSION. SUCH CONVERSION
1322C IS DONE BY DXC210 AND REQUIRES A CONSTANT THAT IS STORED TO SOME FIXED
1323C PRECISION. THE STORED CONSTANT (LOG102 IN THIS ROUTINE) PROVIDES
1324C FOR CONVERSIONS ACCURATE TO THE LAST DECIMAL DIGIT WHEN THE INTEGER
1325C WORD LENGTH DOES NOT EXCEED 63. A LOWER LIMIT OF 15 BITS IS IMPOSED
1326C BECAUSE THE SOFTWARE IS DESIGNED TO RUN ON COMPUTERS WITH INTEGER WORD
1327C LENGTH OF AT LEAST 16 BITS.
1328      IF (15.LE.NBITSX .AND. NBITSX.LE.63) GO TO 30
1329*      CALL XERMSG ('SLATEC', 'DXSET', 'IMPROPER VALUE OF NBITS', 203, 1)
1330      IERROR=203
1331      RETURN
1332   30 CONTINUE
1333      KMAX = 2**(NBITSX-1) - L2
1334      NB = (NBITSX-1)/2
1335      MLG102 = 2**NB
1336      IF (1.LE.NRDPLC*LOG2R .AND. NRDPLC*LOG2R.LE.120) GO TO 40
1337*      CALL XERMSG ('SLATEC', 'DXSET', 'IMPROPER VALUE OF NRADPL', 204,
1338*     +             1)
1339      IERROR=204
1340      RETURN
1341   40 CONTINUE
1342      NLG102 = NRDPLC*LOG2R/NB + 3
1343      NP1 = NLG102 + 1
1344C
1345C   AFTER COMPLETION OF THE FOLLOWING LOOP, IC CONTAINS
1346C THE INTEGER PART AND LGTEMP CONTAINS THE FRACTIONAL PART
1347C OF LOG10(IRADX) IN RADIX 1000.
1348      IC = 0
1349      DO 50 II=1,20
1350        I = 21 - II
1351        IT = LOG2R*LOG102(I) + IC
1352        IC = IT/1000
1353        LGTEMP(I) = MOD(IT,1000)
1354   50 CONTINUE
1355C
1356C   AFTER COMPLETION OF THE FOLLOWING LOOP, LG102 CONTAINS
1357C LOG10(IRADX) IN RADIX MLG102. THE RADIX POINT IS
1358C BETWEEN LG102(1) AND LG102(2).
1359      LG102(1) = IC
1360      DO 80 I=2,NP1
1361        LG102X = 0
1362        DO 70 J=1,NB
1363          IC = 0
1364          DO 60 KK=1,20
1365            K = 21 - KK
1366            IT = 2*LGTEMP(K) + IC
1367            IC = IT/1000
1368            LGTEMP(K) = MOD(IT,1000)
1369   60     CONTINUE
1370          LG102X = 2*LG102X + IC
1371   70   CONTINUE
1372        LG102(I) = LG102X
1373   80 CONTINUE
1374C
1375C CHECK SPECIAL CONDITIONS REQUIRED BY SUBROUTINES...
1376      IF (NRDPLC.LT.L) GO TO 90
1377*      CALL XERMSG ('SLATEC', 'DXSET', 'NRADPL .GE. L', 205, 1)
1378      IERROR=205
1379      RETURN
1380   90 IF (6*L.LE.KMAX) GO TO 100
1381*      CALL XERMSG ('SLATEC', 'DXSET', '6*L .GT. KMAX', 206, 1)
1382      IERROR=206
1383      RETURN
1384  100 CONTINUE
1385      IFLAG = 1
1386      RETURN
1387      END
1388
1389      SUBROUTINE DXADD (X, IX, Y, IY, Z, IZ, IERROR)
1390C***BEGIN PROLOGUE  DXADD
1391C***PURPOSE  To provide double-precision floating-point arithmetic
1392C            with an extended exponent range.
1393C***LIBRARY   SLATEC
1394C***CATEGORY  A3D
1395C***TYPE      DOUBLE PRECISION (XADD-S, DXADD-D)
1396C***KEYWORDS  EXTENDED-RANGE DOUBLE-PRECISION ARITHMETIC
1397C***AUTHOR  Lozier, Daniel W., (National Bureau of Standards)
1398C           Smith, John M., (NBS and George Mason University)
1399C***DESCRIPTION
1400C     DOUBLE PRECISION X, Y, Z
1401C     INTEGER IX, IY, IZ
1402C
1403C                  FORMS THE EXTENDED-RANGE SUM  (Z,IZ) =
1404C                  (X,IX) + (Y,IY).  (Z,IZ) IS ADJUSTED
1405C                  BEFORE RETURNING. THE INPUT OPERANDS
1406C                  NEED NOT BE IN ADJUSTED FORM, BUT THEIR
1407C                  PRINCIPAL PARTS MUST SATISFY
1408C                  RADIX**(-2L).LE.ABS(X).LE.RADIX**(2L),
1409C                  RADIX**(-2L).LE.ABS(Y).LE.RADIX**(2L).
1410C
1411C***SEE ALSO  DXSET
1412C***REFERENCES  (NONE)
1413C***ROUTINES CALLED  DXADJ
1414C***COMMON BLOCKS    DXBLK2
1415C***REVISION HISTORY  (YYMMDD)
1416C   820712  DATE WRITTEN
1417C   881020  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
1418C   901019  Revisions to prologue.  (DWL and WRB)
1419C   901106  Changed all specific intrinsics to generic.  (WRB)
1420C           Corrected order of sections in prologue and added TYPE
1421C           section.  (WRB)
1422C   920127  Revised PURPOSE section of prologue.  (DWL)
1423C***END PROLOGUE  DXADD
1424      DOUBLE PRECISION X, Y, Z
1425      INTEGER IX, IY, IZ
1426      DOUBLE PRECISION RADIX, RADIXL, RAD2L, DLG10R
1427      INTEGER L, L2, KMAX
1428      COMMON /DXBLK2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX
1429      SAVE /DXBLK2/
1430      DOUBLE PRECISION S, T
1431C
1432C   THE CONDITIONS IMPOSED ON L AND KMAX BY THIS SUBROUTINE
1433C ARE
1434C     (1) 1 .LT. L .LE. 0.5D0*LOGR(0.5D0*DZERO)
1435C
1436C     (2) NRADPL .LT. L .LE. KMAX/6
1437C
1438C     (3) KMAX .LE. (2**NBITS - 4*L - 1)/2
1439C
1440C THESE CONDITIONS MUST BE MET BY APPROPRIATE CODING
1441C IN SUBROUTINE DXSET.
1442C
1443C***FIRST EXECUTABLE STATEMENT  DXADD
1444      IERROR=0
1445      IF (X.NE.0.0D0) GO TO 10
1446      Z = Y
1447      IZ = IY
1448      GO TO 220
1449   10 IF (Y.NE.0.0D0) GO TO 20
1450      Z = X
1451      IZ = IX
1452      GO TO 220
1453   20 CONTINUE
1454      IF (IX.GE.0 .AND. IY.GE.0) GO TO 40
1455      IF (IX.LT.0 .AND. IY.LT.0) GO TO 40
1456      IF (ABS(IX).LE.6*L .AND. ABS(IY).LE.6*L) GO TO 40
1457      IF (IX.GE.0) GO TO 30
1458      Z = Y
1459      IZ = IY
1460      GO TO 220
1461   30 CONTINUE
1462      Z = X
1463      IZ = IX
1464      GO TO 220
1465   40 I = IX - IY
1466      if (I .lt. 0) then
1467         goto 80
1468      elseif (I .eq. 0) then
1469         goto 50
1470      else
1471         goto 90
1472      endif
1473   50 IF (ABS(X).GT.1.0D0 .AND. ABS(Y).GT.1.0D0) GO TO 60
1474      IF (ABS(X).LT.1.0D0 .AND. ABS(Y).LT.1.0D0) GO TO 70
1475      Z = X + Y
1476      IZ = IX
1477      GO TO 220
1478   60 S = X/RADIXL
1479      T = Y/RADIXL
1480      Z = S + T
1481      IZ = IX + L
1482      GO TO 220
1483   70 S = X*RADIXL
1484      T = Y*RADIXL
1485      Z = S + T
1486      IZ = IX - L
1487      GO TO 220
1488   80 S = Y
1489      IS = IY
1490      T = X
1491      GO TO 100
1492   90 S = X
1493      IS = IX
1494      T = Y
1495  100 CONTINUE
1496C
1497C  AT THIS POINT, THE ONE OF (X,IX) OR (Y,IY) THAT HAS THE
1498C LARGER AUXILIARY INDEX IS STORED IN (S,IS). THE PRINCIPAL
1499C PART OF THE OTHER INPUT IS STORED IN T.
1500C
1501      I1 = ABS(I)/L
1502      I2 = MOD(ABS(I),L)
1503      IF (ABS(T).GE.RADIXL) GO TO 130
1504      IF (ABS(T).GE.1.0D0) GO TO 120
1505      IF (RADIXL*ABS(T).GE.1.0D0) GO TO 110
1506      J = I1 + 1
1507      T = T*RADIX**(L-I2)
1508      GO TO 140
1509  110 J = I1
1510      T = T*RADIX**(-I2)
1511      GO TO 140
1512  120 J = I1 - 1
1513      IF (J.LT.0) GO TO 110
1514      T = T*RADIX**(-I2)/RADIXL
1515      GO TO 140
1516  130 J = I1 - 2
1517      IF (J.LT.0) GO TO 120
1518      T = T*RADIX**(-I2)/RAD2L
1519  140 CONTINUE
1520C
1521C  AT THIS POINT, SOME OR ALL OF THE DIFFERENCE IN THE
1522C AUXILIARY INDICES HAS BEEN USED TO EFFECT A LEFT SHIFT
1523C OF T.  THE SHIFTED VALUE OF T SATISFIES
1524C
1525C       RADIX**(-2*L) .LE. ABS(T) .LE. 1.0D0
1526C
1527C AND, IF J=0, NO FURTHER SHIFTING REMAINS TO BE DONE.
1528C
1529      IF (J.EQ.0) GO TO 190
1530      IF (ABS(S).GE.RADIXL .OR. J.GT.3) GO TO 150
1531      IF (ABS(S).GE.1.0D0) GO TO (180, 150, 150), J
1532      IF (RADIXL*ABS(S).GE.1.0D0) GO TO (180, 170, 150), J
1533      GO TO (180, 170, 160), J
1534  150 Z = S
1535      IZ = IS
1536      GO TO 220
1537  160 S = S*RADIXL
1538  170 S = S*RADIXL
1539  180 S = S*RADIXL
1540  190 CONTINUE
1541C
1542C   AT THIS POINT, THE REMAINING DIFFERENCE IN THE
1543C AUXILIARY INDICES HAS BEEN USED TO EFFECT A RIGHT SHIFT
1544C OF S.  IF THE SHIFTED VALUE OF S WOULD HAVE EXCEEDED
1545C RADIX**L, THEN (S,IS) IS RETURNED AS THE VALUE OF THE
1546C SUM.
1547C
1548      IF (ABS(S).GT.1.0D0 .AND. ABS(T).GT.1.0D0) GO TO 200
1549      IF (ABS(S).LT.1.0D0 .AND. ABS(T).LT.1.0D0) GO TO 210
1550      Z = S + T
1551      IZ = IS - J*L
1552      GO TO 220
1553  200 S = S/RADIXL
1554      T = T/RADIXL
1555      Z = S + T
1556      IZ = IS - J*L + L
1557      GO TO 220
1558  210 S = S*RADIXL
1559      T = T*RADIXL
1560      Z = S + T
1561      IZ = IS - J*L - L
1562  220 CALL DXADJ(Z, IZ,IERROR)
1563      RETURN
1564      END
1565*DECK DXADJ
1566      SUBROUTINE DXADJ (X, IX, IERROR)
1567C***BEGIN PROLOGUE  DXADJ
1568C***PURPOSE  To provide double-precision floating-point arithmetic
1569C            with an extended exponent range.
1570C***LIBRARY   SLATEC
1571C***CATEGORY  A3D
1572C***TYPE      DOUBLE PRECISION (XADJ-S, DXADJ-D)
1573C***KEYWORDS  EXTENDED-RANGE DOUBLE-PRECISION ARITHMETIC
1574C***AUTHOR  Lozier, Daniel W., (National Bureau of Standards)
1575C           Smith, John M., (NBS and George Mason University)
1576C***DESCRIPTION
1577C     DOUBLE PRECISION X
1578C     INTEGER IX
1579C
1580C                  TRANSFORMS (X,IX) SO THAT
1581C                  RADIX**(-L) .LE. ABS(X) .LT. RADIX**L.
1582C                  ON MOST COMPUTERS THIS TRANSFORMATION DOES
1583C                  NOT CHANGE THE MANTISSA OF X PROVIDED RADIX IS
1584C                  THE NUMBER BASE OF DOUBLE-PRECISION ARITHMETIC.
1585C
1586C***SEE ALSO  DXSET
1587C***REFERENCES  (NONE)
1588C***ROUTINES CALLED  XERMSG
1589C***COMMON BLOCKS    DXBLK2
1590C***REVISION HISTORY  (YYMMDD)
1591C   820712  DATE WRITTEN
1592C   881020  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
1593C   901019  Revisions to prologue.  (DWL and WRB)
1594C   901106  Changed all specific intrinsics to generic.  (WRB)
1595C           Corrected order of sections in prologue and added TYPE
1596C           section.  (WRB)
1597C           CALLs to XERROR changed to CALLs to XERMSG.  (WRB)
1598C   920127  Revised PURPOSE section of prologue.  (DWL)
1599C***END PROLOGUE  DXADJ
1600      DOUBLE PRECISION X
1601      INTEGER IX
1602      DOUBLE PRECISION RADIX, RADIXL, RAD2L, DLG10R
1603      INTEGER L, L2, KMAX
1604      COMMON /DXBLK2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX
1605      SAVE /DXBLK2/
1606C
1607C   THE CONDITION IMPOSED ON L AND KMAX BY THIS SUBROUTINE
1608C IS
1609C     2*L .LE. KMAX
1610C
1611C THIS CONDITION MUST BE MET BY APPROPRIATE CODING
1612C IN SUBROUTINE DXSET.
1613C
1614C***FIRST EXECUTABLE STATEMENT  DXADJ
1615      IERROR=0
1616      IF (X.EQ.0.0D0) GO TO 50
1617      IF (ABS(X).GE.1.0D0) GO TO 20
1618      IF (RADIXL*ABS(X).GE.1.0D0) GO TO 60
1619      X = X*RAD2L
1620      IF (IX.LT.0) GO TO 10
1621      IX = IX - L2
1622      GO TO 70
1623   10 IF (IX.LT.-KMAX+L2) GO TO 40
1624      IX = IX - L2
1625      GO TO 70
1626   20 IF (ABS(X).LT.RADIXL) GO TO 60
1627      X = X/RAD2L
1628      IF (IX.GT.0) GO TO 30
1629      IX = IX + L2
1630      GO TO 70
1631   30 IF (IX.GT.KMAX-L2) GO TO 40
1632      IX = IX + L2
1633      GO TO 70
1634 40   continue
1635*   40 CALL XERMSG ('SLATEC', 'DXADJ', 'overflow in auxiliary index',
1636*     +             207, 1)
1637      IERROR=207
1638      RETURN
1639   50 IX = 0
1640   60 IF (ABS(IX).GT.KMAX) GO TO 40
1641   70 RETURN
1642      END
1643