github.com/gopherd/gonum@v0.0.4/mathext/internal/amos/amoslib/zuni1.f (about) 1 SUBROUTINE ZUNI1(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, 2 * TOL, ELIM, ALIM) 3 C***BEGIN PROLOGUE ZUNI1 4 C***REFER TO ZBESI,ZBESK 5 C 6 C ZUNI1 COMPUTES I(FNU,Z) BY MEANS OF THE UNIFORM ASYMPTOTIC 7 C EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3. 8 C 9 C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC 10 C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. 11 C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER 12 C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. 13 C Y(I)=CZERO FOR I=NLAST+1,N 14 C 15 C***ROUTINES CALLED ZUCHK,ZUNIK,ZUOIK,D1MACH,ZABS 16 C***END PROLOGUE ZUNI1 17 C COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1, 18 C *S2,Y,Z,ZETA1,ZETA2 19 DOUBLE PRECISION ALIM, APHI, ASCLE, BRY, CONER, CRSC, 20 * CSCL, CSRR, CSSR, CWRKI, CWRKR, C1R, C2I, C2M, C2R, ELIM, FN, 21 * FNU, FNUL, PHII, PHIR, RAST, RS1, RZI, RZR, STI, STR, SUMI, 22 * SUMR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I, 23 * ZETA1R, ZETA2I, ZETA2R, ZI, ZR, CYR, CYI, D1MACH, ZABS 24 INTEGER I, IFLAG, INIT, K, KODE, M, N, ND, NLAST, NN, NUF, NW, NZ 25 DIMENSION BRY(3), YR(N), YI(N), CWRKR(16), CWRKI(16), CSSR(3), 26 * CSRR(3), CYR(2), CYI(2) 27 DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 / 28 C 29 NZ = 0 30 ND = N 31 NLAST = 0 32 C----------------------------------------------------------------------- 33 C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG- 34 C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE, 35 C EXP(ALIM)=EXP(ELIM)*TOL 36 C----------------------------------------------------------------------- 37 CSCL = 1.0D0/TOL 38 CRSC = TOL 39 CSSR(1) = CSCL 40 CSSR(2) = CONER 41 CSSR(3) = CRSC 42 CSRR(1) = CRSC 43 CSRR(2) = CONER 44 CSRR(3) = CSCL 45 BRY(1) = 1.0D+3*D1MACH(1)/TOL 46 C----------------------------------------------------------------------- 47 C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER 48 C----------------------------------------------------------------------- 49 FN = DMAX1(FNU,1.0D0) 50 INIT = 0 51 CALL ZUNIK(ZR, ZI, FN, 1, 1, TOL, INIT, PHIR, PHII, ZETA1R, 52 * ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI) 53 IF (KODE.EQ.1) GO TO 10 54 STR = ZR + ZETA2R 55 STI = ZI + ZETA2I 56 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 57 STR = STR*RAST*RAST 58 STI = -STI*RAST*RAST 59 S1R = -ZETA1R + STR 60 S1I = -ZETA1I + STI 61 GO TO 20 62 10 CONTINUE 63 S1R = -ZETA1R + ZETA2R 64 S1I = -ZETA1I + ZETA2I 65 20 CONTINUE 66 RS1 = S1R 67 IF (DABS(RS1).GT.ELIM) GO TO 130 68 30 CONTINUE 69 NN = MIN0(2,ND) 70 DO 80 I=1,NN 71 FN = FNU + DBLE(FLOAT(ND-I)) 72 INIT = 0 73 CALL ZUNIK(ZR, ZI, FN, 1, 0, TOL, INIT, PHIR, PHII, ZETA1R, 74 * ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI) 75 IF (KODE.EQ.1) GO TO 40 76 STR = ZR + ZETA2R 77 STI = ZI + ZETA2I 78 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 79 STR = STR*RAST*RAST 80 STI = -STI*RAST*RAST 81 S1R = -ZETA1R + STR 82 S1I = -ZETA1I + STI + ZI 83 GO TO 50 84 40 CONTINUE 85 S1R = -ZETA1R + ZETA2R 86 S1I = -ZETA1I + ZETA2I 87 50 CONTINUE 88 C----------------------------------------------------------------------- 89 C TEST FOR UNDERFLOW AND OVERFLOW 90 C----------------------------------------------------------------------- 91 RS1 = S1R 92 IF (DABS(RS1).GT.ELIM) GO TO 110 93 IF (I.EQ.1) IFLAG = 2 94 IF (DABS(RS1).LT.ALIM) GO TO 60 95 C----------------------------------------------------------------------- 96 C REFINE TEST AND SCALE 97 C----------------------------------------------------------------------- 98 APHI = ZABS(CMPLX(PHIR,PHII,kind=KIND(1.0D0))) 99 RS1 = RS1 + DLOG(APHI) 100 IF (DABS(RS1).GT.ELIM) GO TO 110 101 IF (I.EQ.1) IFLAG = 1 102 IF (RS1.LT.0.0D0) GO TO 60 103 IF (I.EQ.1) IFLAG = 3 104 60 CONTINUE 105 C----------------------------------------------------------------------- 106 C SCALE S1 IF CABS(S1).LT.ASCLE 107 C----------------------------------------------------------------------- 108 S2R = PHIR*SUMR - PHII*SUMI 109 S2I = PHIR*SUMI + PHII*SUMR 110 STR = DEXP(S1R)*CSSR(IFLAG) 111 S1R = STR*DCOS(S1I) 112 S1I = STR*DSIN(S1I) 113 STR = S2R*S1R - S2I*S1I 114 S2I = S2R*S1I + S2I*S1R 115 S2R = STR 116 IF (IFLAG.NE.1) GO TO 70 117 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) 118 IF (NW.NE.0) GO TO 110 119 70 CONTINUE 120 CYR(I) = S2R 121 CYI(I) = S2I 122 M = ND - I + 1 123 YR(M) = S2R*CSRR(IFLAG) 124 YI(M) = S2I*CSRR(IFLAG) 125 80 CONTINUE 126 IF (ND.LE.2) GO TO 100 127 RAST = 1.0D0/ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 128 STR = ZR*RAST 129 STI = -ZI*RAST 130 RZR = (STR+STR)*RAST 131 RZI = (STI+STI)*RAST 132 BRY(2) = 1.0D0/BRY(1) 133 BRY(3) = D1MACH(2) 134 S1R = CYR(1) 135 S1I = CYI(1) 136 S2R = CYR(2) 137 S2I = CYI(2) 138 C1R = CSRR(IFLAG) 139 ASCLE = BRY(IFLAG) 140 K = ND - 2 141 FN = DBLE(FLOAT(K)) 142 DO 90 I=3,ND 143 C2R = S2R 144 C2I = S2I 145 S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I) 146 S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R) 147 S1R = C2R 148 S1I = C2I 149 C2R = S2R*C1R 150 C2I = S2I*C1R 151 YR(K) = C2R 152 YI(K) = C2I 153 K = K - 1 154 FN = FN - 1.0D0 155 IF (IFLAG.GE.3) GO TO 90 156 STR = DABS(C2R) 157 STI = DABS(C2I) 158 C2M = DMAX1(STR,STI) 159 IF (C2M.LE.ASCLE) GO TO 90 160 IFLAG = IFLAG + 1 161 ASCLE = BRY(IFLAG) 162 S1R = S1R*C1R 163 S1I = S1I*C1R 164 S2R = C2R 165 S2I = C2I 166 S1R = S1R*CSSR(IFLAG) 167 S1I = S1I*CSSR(IFLAG) 168 S2R = S2R*CSSR(IFLAG) 169 S2I = S2I*CSSR(IFLAG) 170 C1R = CSRR(IFLAG) 171 90 CONTINUE 172 100 CONTINUE 173 RETURN 174 C----------------------------------------------------------------------- 175 C SET UNDERFLOW AND UPDATE PARAMETERS 176 C----------------------------------------------------------------------- 177 110 CONTINUE 178 IF (RS1.GT.0.0D0) GO TO 120 179 YR(ND) = ZEROR 180 YI(ND) = ZEROI 181 NZ = NZ + 1 182 ND = ND - 1 183 IF (ND.EQ.0) GO TO 100 184 CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM) 185 IF (NUF.LT.0) GO TO 120 186 ND = ND - NUF 187 NZ = NZ + NUF 188 IF (ND.EQ.0) GO TO 100 189 FN = FNU + DBLE(FLOAT(ND-1)) 190 IF (FN.GE.FNUL) GO TO 30 191 NLAST = ND 192 RETURN 193 120 CONTINUE 194 NZ = -1 195 RETURN 196 130 CONTINUE 197 IF (RS1.GT.0.0D0) GO TO 120 198 NZ = N 199 DO 140 I=1,N 200 YR(I) = ZEROR 201 YI(I) = ZEROI 202 140 CONTINUE 203 RETURN 204 END