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