1// Copyright ©2016 The Gonum Authors. All rights reserved. 2// Use of this source code is governed by a BSD-style 3// license that can be found in the LICENSE file. 4 5package amos 6 7import ( 8 "math" 9 "math/cmplx" 10) 11 12// These routines are the versions directly modified from the Fortran code. 13// They are used to ensure that code style improvements do not change the 14// code output. 15 16func iabs(a int) int { 17 if a >= 0 { 18 return a 19 } 20 return -a 21} 22 23func min0(a, b int) int { 24 if a < b { 25 return a 26 } 27 return b 28} 29 30func max0(a, b int) int { 31 if a > b { 32 return a 33 } 34 return b 35} 36 37func zairyOrig(ZR, ZI float64, ID, KODE int) (AIR, AII float64, NZ int) { 38 // zairy is adapted from the original Netlib code by Donald Amos. 39 // http://www.netlib.no/netlib/amos/zairy.f 40 41 // Original comment: 42 /* 43 C***BEGIN PROLOGUE ZAIRY 44 C***DATE WRITTEN 830501 (YYMMDD) 45 C***REVISION DATE 890801 (YYMMDD) 46 C***CATEGORY NO. B5K 47 C***KEYWORDS AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD 48 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 49 C***PURPOSE TO COMPUTE AIRY FUNCTIONS AI(Z) AND DAI(Z) FOR COMPLEX Z 50 C***DESCRIPTION 51 C 52 C ***A DOUBLE PRECISION ROUTINE*** 53 C ON KODE=1, ZAIRY COMPUTES THE COMPLEX AIRY FUNCTION AI(Z) OR 54 C ITS DERIVATIVE DAI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON 55 C KODE=2, A SCALING OPTION CEXP(ZTA)*AI(Z) OR CEXP(ZTA)* 56 C DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN 57 C -PI/3<ARG(Z)<PI/3 AND THE EXPONENTIAL GROWTH IN 58 C PI/3<ABS(ARG(Z))<PI WHERE ZTA=(2/3)*Z*CSQRT(Z). 59 C 60 C WHILE THE AIRY FUNCTIONS AI(Z) AND DAI(Z)/DZ ARE ANALYTIC IN 61 C THE WHOLE Z PLANE, THE CORRESPONDING SCALED FUNCTIONS DEFINED 62 C FOR KODE=2 HAVE A CUT ALONG THE NEGATIVE REAL AXIS. 63 C DEFINTIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF 64 C MATHEMATICAL FUNCTIONS (REF. 1). 65 C 66 C INPUT ZR,ZI ARE DOUBLE PRECISION 67 C ZR,ZI - Z=CMPLX(ZR,ZI) 68 C ID - ORDER OF DERIVATIVE, ID=0 OR ID=1 69 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 70 C KODE= 1 returnS 71 C AI=AI(Z) ON ID=0 OR 72 C AI=DAI(Z)/DZ ON ID=1 73 C = 2 returnS 74 C AI=CEXP(ZTA)*AI(Z) ON ID=0 OR 75 C AI=CEXP(ZTA)*DAI(Z)/DZ ON ID=1 WHERE 76 C ZTA=(2/3)*Z*CSQRT(Z) 77 C 78 C OUTPUT AIR,AII ARE DOUBLE PRECISION 79 C AIR,AII- COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND 80 C KODE 81 C NZ - UNDERFLOW INDICATOR 82 C NZ= 0 , NORMAL return 83 C NZ= 1 , AI=CMPLX(0.0E0,0.0E0) DUE TO UNDERFLOW IN 84 C -PI/3<ARG(Z)<PI/3 ON KODE=1 85 C IERR - ERROR FLAG 86 C IERR=0, NORMAL return - COMPUTATION COMPLETED 87 C IERR=1, INPUT ERROR - NO COMPUTATION 88 C IERR=2, OVERFLOW - NO COMPUTATION, REAL(ZTA) 89 C TOO LARGE ON KODE=1 90 C IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED 91 C LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION 92 C PRODUCE LESS THAN HALF OF MACHINE ACCURACY 93 C IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION 94 C COMPLETE LOSS OF ACCURACY BY ARGUMENT 95 C REDUCTION 96 C IERR=5, ERROR - NO COMPUTATION, 97 C ALGORITHM TERMINATION CONDITION NOT MET 98 C 99 C***LONG DESCRIPTION 100 C 101 C AI AND DAI ARE COMPUTED FOR CABS(Z)>1.0 FROM THE K BESSEL 102 C FUNCTIONS BY 103 C 104 C AI(Z)=C*SQRT(Z)*K(1/3,ZTA) , DAI(Z)=-C*Z*K(2/3,ZTA) 105 C C=1.0/(PI*SQRT(3.0)) 106 C ZTA=(2/3)*Z**(3/2) 107 C 108 C WITH THE POWER SERIES FOR CABS(Z)<=1.0. 109 C 110 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 111 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES 112 C OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF 113 C THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR), 114 C THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR 115 C FLAG IERR=3 IS TRIGGERED WHERE UR=dmax(dmach[4),1.0D-18) IS 116 C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. 117 C ALSO, if THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN 118 C ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT 119 C FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE 120 C LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA 121 C MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, 122 C AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE 123 C PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE 124 C PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT- 125 C ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG- 126 C NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN 127 C DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN 128 C EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, 129 C NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE 130 C PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER 131 C MACHINES. 132 C 133 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 134 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 135 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 136 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 137 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 138 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 139 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 140 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 141 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 142 C SEVERAL ORDERS OF MAGNITUDE. if ONE COMPONENT IS 10**K LARGER 143 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 144 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 145 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 146 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 147 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 148 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 149 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 150 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 151 C OR -PI/2+P. 152 C 153 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 154 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 155 C COMMERCE, 1955. 156 C 157 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 158 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983 159 C 160 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 161 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 162 C 1018, MAY, 1985 163 C 164 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 165 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 166 C MATH. SOFTWARE, 1986 167 */ 168 var AI, CONE, CSQ, CY, S1, S2, TRM1, TRM2, Z, ZTA, Z3 complex128 169 var AA, AD, AK, ALIM, ATRM, AZ, AZ3, BK, 170 CC, CK, COEF, CONEI, CONER, CSQI, CSQR, C1, C2, DIG, 171 DK, D1, D2, ELIM, FID, FNU, PTR, RL, R1M5, SFAC, STI, STR, 172 S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I, TRM2R, TTH, ZEROI, 173 ZEROR, ZTAI, ZTAR, Z3I, Z3R, ALAZ, BB float64 174 var IERR, IFLAG, K, K1, K2, MR, NN int 175 var tmp complex128 176 177 // Extra element for padding. 178 CYR := []float64{math.NaN(), 0} 179 CYI := []float64{math.NaN(), 0} 180 181 _ = AI 182 _ = CONE 183 _ = CSQ 184 _ = CY 185 _ = S1 186 _ = S2 187 _ = TRM1 188 _ = TRM2 189 _ = Z 190 _ = ZTA 191 _ = Z3 192 193 TTH = 6.66666666666666667e-01 194 C1 = 3.55028053887817240e-01 195 C2 = 2.58819403792806799e-01 196 COEF = 1.83776298473930683e-01 197 ZEROR = 0 198 ZEROI = 0 199 CONER = 1 200 CONEI = 0 201 202 NZ = 0 203 if ID < 0 || ID > 1 { 204 IERR = 1 205 } 206 if KODE < 1 || KODE > 2 { 207 IERR = 1 208 } 209 if IERR != 0 { 210 return 211 } 212 AZ = zabs(complex(ZR, ZI)) 213 TOL = dmax(dmach[4], 1.0e-18) 214 FID = float64(ID) 215 if AZ > 1.0e0 { 216 goto Seventy 217 } 218 219 // POWER SERIES FOR CABS(Z)<=1. 220 S1R = CONER 221 S1I = CONEI 222 S2R = CONER 223 S2I = CONEI 224 if AZ < TOL { 225 goto OneSeventy 226 } 227 AA = AZ * AZ 228 if AA < TOL/AZ { 229 goto Forty 230 } 231 TRM1R = CONER 232 TRM1I = CONEI 233 TRM2R = CONER 234 TRM2I = CONEI 235 ATRM = 1.0e0 236 STR = ZR*ZR - ZI*ZI 237 STI = ZR*ZI + ZI*ZR 238 Z3R = STR*ZR - STI*ZI 239 Z3I = STR*ZI + STI*ZR 240 AZ3 = AZ * AA 241 AK = 2.0e0 + FID 242 BK = 3.0e0 - FID - FID 243 CK = 4.0e0 - FID 244 DK = 3.0e0 + FID + FID 245 D1 = AK * DK 246 D2 = BK * CK 247 AD = dmin(D1, D2) 248 AK = 24.0e0 + 9.0e0*FID 249 BK = 30.0e0 - 9.0e0*FID 250 for K = 1; K <= 25; K++ { 251 STR = (TRM1R*Z3R - TRM1I*Z3I) / D1 252 TRM1I = (TRM1R*Z3I + TRM1I*Z3R) / D1 253 TRM1R = STR 254 S1R = S1R + TRM1R 255 S1I = S1I + TRM1I 256 STR = (TRM2R*Z3R - TRM2I*Z3I) / D2 257 TRM2I = (TRM2R*Z3I + TRM2I*Z3R) / D2 258 TRM2R = STR 259 S2R = S2R + TRM2R 260 S2I = S2I + TRM2I 261 ATRM = ATRM * AZ3 / AD 262 D1 = D1 + AK 263 D2 = D2 + BK 264 AD = dmin(D1, D2) 265 if ATRM < TOL*AD { 266 goto Forty 267 } 268 AK = AK + 18.0e0 269 BK = BK + 18.0e0 270 } 271Forty: 272 if ID == 1 { 273 goto Fifty 274 } 275 AIR = S1R*C1 - C2*(ZR*S2R-ZI*S2I) 276 AII = S1I*C1 - C2*(ZR*S2I+ZI*S2R) 277 if KODE == 1 { 278 return 279 } 280 tmp = zsqrt(complex(ZR, ZI)) 281 STR = real(tmp) 282 STI = imag(tmp) 283 ZTAR = TTH * (ZR*STR - ZI*STI) 284 ZTAI = TTH * (ZR*STI + ZI*STR) 285 tmp = zexp(complex(ZTAR, ZTAI)) 286 STR = real(tmp) 287 STI = imag(tmp) 288 PTR = AIR*STR - AII*STI 289 AII = AIR*STI + AII*STR 290 AIR = PTR 291 return 292 293Fifty: 294 AIR = -S2R * C2 295 AII = -S2I * C2 296 if AZ <= TOL { 297 goto Sixty 298 } 299 STR = ZR*S1R - ZI*S1I 300 STI = ZR*S1I + ZI*S1R 301 CC = C1 / (1.0e0 + FID) 302 AIR = AIR + CC*(STR*ZR-STI*ZI) 303 AII = AII + CC*(STR*ZI+STI*ZR) 304 305Sixty: 306 if KODE == 1 { 307 return 308 } 309 tmp = zsqrt(complex(ZR, ZI)) 310 STR = real(tmp) 311 STI = imag(tmp) 312 ZTAR = TTH * (ZR*STR - ZI*STI) 313 ZTAI = TTH * (ZR*STI + ZI*STR) 314 tmp = zexp(complex(ZTAR, ZTAI)) 315 STR = real(tmp) 316 STI = imag(tmp) 317 PTR = STR*AIR - STI*AII 318 AII = STR*AII + STI*AIR 319 AIR = PTR 320 return 321 322 // CASE FOR CABS(Z)>1.0. 323Seventy: 324 FNU = (1.0e0 + FID) / 3.0e0 325 326 /* 327 SET PARAMETERS RELATED TO MACHINE CONSTANTS. 328 TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0D-18. 329 ELIM IS THE APPROXIMATE EXPONENTIAL OVER-&&UNDERFLOW LIMIT. 330 EXP(-ELIM)<EXP(-ALIM)=EXP(-ELIM)/TOL AND 331 EXP(ELIM)>EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 332 UNDERFLOW&&OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 333 RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LA>=Z. 334 DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 335 */ 336 K1 = imach[15] 337 K2 = imach[16] 338 R1M5 = dmach[5] 339 340 K = min0(iabs(K1), iabs(K2)) 341 ELIM = 2.303e0 * (float64(K)*R1M5 - 3.0e0) 342 K1 = imach[14] - 1 343 AA = R1M5 * float64(K1) 344 DIG = dmin(AA, 18.0e0) 345 AA = AA * 2.303e0 346 ALIM = ELIM + dmax(-AA, -41.45e0) 347 RL = 1.2e0*DIG + 3.0e0 348 ALAZ = dlog(AZ) 349 350 // TEST FOR PROPER RANGE. 351 AA = 0.5e0 / TOL 352 BB = float64(float32(imach[9])) * 0.5e0 353 AA = dmin(AA, BB) 354 AA = math.Pow(AA, TTH) 355 if AZ > AA { 356 goto TwoSixty 357 } 358 AA = dsqrt(AA) 359 if AZ > AA { 360 IERR = 3 361 } 362 tmp = zsqrt(complex(ZR, ZI)) 363 CSQR = real(tmp) 364 CSQI = imag(tmp) 365 ZTAR = TTH * (ZR*CSQR - ZI*CSQI) 366 ZTAI = TTH * (ZR*CSQI + ZI*CSQR) 367 368 // RE(ZTA)<=0 WHEN RE(Z)<0, ESPECIALLY WHEN IM(Z) IS SMALL. 369 IFLAG = 0 370 SFAC = 1.0e0 371 AK = ZTAI 372 if ZR >= 0.0e0 { 373 goto Eighty 374 } 375 BK = ZTAR 376 CK = -dabs(BK) 377 ZTAR = CK 378 ZTAI = AK 379 380Eighty: 381 if ZI != 0.0e0 { 382 goto Ninety 383 } 384 if ZR > 0.0e0 { 385 goto Ninety 386 } 387 ZTAR = 0.0e0 388 ZTAI = AK 389Ninety: 390 AA = ZTAR 391 if AA >= 0.0e0 && ZR > 0.0e0 { 392 goto OneTen 393 } 394 if KODE == 2 { 395 goto OneHundred 396 } 397 398 // OVERFLOW TEST. 399 if AA > (-ALIM) { 400 goto OneHundred 401 } 402 AA = -AA + 0.25e0*ALAZ 403 IFLAG = 1 404 SFAC = TOL 405 if AA > ELIM { 406 goto TwoSeventy 407 } 408 409OneHundred: 410 // CBKNU AND CACON return EXP(ZTA)*K(FNU,ZTA) ON KODE=2. 411 MR = 1 412 if ZI < 0.0e0 { 413 MR = -1 414 } 415 ZTAR, ZTAI, FNU, KODE, MR, _, CYR, CYI, NN, RL, TOL, ELIM, ALIM = zacaiOrig(ZTAR, ZTAI, FNU, KODE, MR, 1, CYR, CYI, NN, RL, TOL, ELIM, ALIM) 416 if NN < 0 { 417 goto TwoEighty 418 } 419 NZ = NZ + NN 420 goto OneThirty 421 422OneTen: 423 if KODE == 2 { 424 goto OneTwenty 425 } 426 427 // UNDERFLOW TEST. 428 if AA < ALIM { 429 goto OneTwenty 430 } 431 AA = -AA - 0.25e0*ALAZ 432 IFLAG = 2 433 SFAC = 1.0e0 / TOL 434 if AA < (-ELIM) { 435 goto TwoTen 436 } 437OneTwenty: 438 ZTAR, ZTAI, FNU, KODE, _, CYR, CYI, NZ, TOL, ELIM, ALIM = zbknuOrig(ZTAR, ZTAI, FNU, KODE, 1, CYR, CYI, NZ, TOL, ELIM, ALIM) 439 440OneThirty: 441 S1R = CYR[1] * COEF 442 S1I = CYI[1] * COEF 443 if IFLAG != 0 { 444 goto OneFifty 445 } 446 if ID == 1 { 447 goto OneFourty 448 } 449 AIR = CSQR*S1R - CSQI*S1I 450 AII = CSQR*S1I + CSQI*S1R 451 return 452OneFourty: 453 AIR = -(ZR*S1R - ZI*S1I) 454 AII = -(ZR*S1I + ZI*S1R) 455 return 456OneFifty: 457 S1R = S1R * SFAC 458 S1I = S1I * SFAC 459 if ID == 1 { 460 goto OneSixty 461 } 462 STR = S1R*CSQR - S1I*CSQI 463 S1I = S1R*CSQI + S1I*CSQR 464 S1R = STR 465 AIR = S1R / SFAC 466 AII = S1I / SFAC 467 return 468OneSixty: 469 STR = -(S1R*ZR - S1I*ZI) 470 S1I = -(S1R*ZI + S1I*ZR) 471 S1R = STR 472 AIR = S1R / SFAC 473 AII = S1I / SFAC 474 return 475OneSeventy: 476 AA = 1.0e+3 * dmach[1] 477 S1R = ZEROR 478 S1I = ZEROI 479 if ID == 1 { 480 goto OneNinety 481 } 482 if AZ <= AA { 483 goto OneEighty 484 } 485 S1R = C2 * ZR 486 S1I = C2 * ZI 487OneEighty: 488 AIR = C1 - S1R 489 AII = -S1I 490 return 491OneNinety: 492 AIR = -C2 493 AII = 0.0e0 494 AA = dsqrt(AA) 495 if AZ <= AA { 496 goto TwoHundred 497 } 498 S1R = 0.5e0 * (ZR*ZR - ZI*ZI) 499 S1I = ZR * ZI 500TwoHundred: 501 AIR = AIR + C1*S1R 502 AII = AII + C1*S1I 503 return 504TwoTen: 505 NZ = 1 506 AIR = ZEROR 507 AII = ZEROI 508 return 509TwoSeventy: 510 NZ = 0 511 IERR = 2 512 return 513TwoEighty: 514 if NN == (-1) { 515 goto TwoSeventy 516 } 517 NZ = 0 518 IERR = 5 519 return 520TwoSixty: 521 IERR = 4 522 NZ = 0 523 return 524} 525 526// sbknu computes the k bessel function in the right half z plane. 527func zbknuOrig(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL, ELIM, ALIM float64) (ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZout int, TOLout, ELIMout, ALIMout float64) { 528 /* Old dimension comment. 529 DIMENSION YR(N), YI(N), CC(8), CSSR(3), CSRR(3), BRY(3), CYR(2), 530 * CYI(2) 531 */ 532 533 // TODO(btracey): Find which of these are inputs/outputs/both and clean up 534 // the function call. 535 // YR and YI have length n (but n+1 with better indexing) 536 var AA, AK, ASCLE, A1, A2, BB, BK, CAZ, 537 CBI, CBR, CCHI, CCHR, CKI, CKR, COEFI, COEFR, CONEI, CONER, 538 CRSCR, CSCLR, CSHI, CSHR, CSI, CSR, CTWOR, 539 CZEROI, CZEROR, CZI, CZR, DNU, DNU2, DPI, ETEST, FC, FHS, 540 FI, FK, FKS, FMUI, FMUR, FPI, FR, G1, G2, HPI, PI, PR, PTI, 541 PTR, P1I, P1R, P2I, P2M, P2R, QI, QR, RAK, RCAZ, RTHPI, RZI, 542 RZR, R1, S, SMUI, SMUR, SPI, STI, STR, S1I, S1R, S2I, S2R, TM, 543 TTH, T1, T2, ELM, CELMR, ZDR, ZDI, AS, ALAS, HELIM float64 544 545 var I, IFLAG, INU, K, KFLAG, KK, KMAX, KODED, IDUM, J, IC, INUB, NW int 546 547 var tmp complex128 548 var CSSR, CSRR, BRY [4]float64 549 var CYR, CYI [3]float64 550 551 KMAX = 30 552 CZEROR = 0 553 CZEROI = 0 554 CONER = 1 555 CONEI = 0 556 CTWOR = 2 557 R1 = 2 558 559 DPI = 3.14159265358979324e0 560 RTHPI = 1.25331413731550025e0 561 SPI = 1.90985931710274403e0 562 HPI = 1.57079632679489662e0 563 FPI = 1.89769999331517738e0 564 TTH = 6.66666666666666666e-01 565 566 CC := [9]float64{math.NaN(), 5.77215664901532861e-01, -4.20026350340952355e-02, 567 -4.21977345555443367e-02, 7.21894324666309954e-03, 568 -2.15241674114950973e-04, -2.01348547807882387e-05, 569 1.13302723198169588e-06, 6.11609510448141582e-09} 570 571 CAZ = zabs(complex(ZR, ZI)) 572 CSCLR = 1.0e0 / TOL 573 CRSCR = TOL 574 CSSR[1] = CSCLR 575 CSSR[2] = 1.0e0 576 CSSR[3] = CRSCR 577 CSRR[1] = CRSCR 578 CSRR[2] = 1.0e0 579 CSRR[3] = CSCLR 580 BRY[1] = 1.0e+3 * dmach[1] / TOL 581 BRY[2] = 1.0e0 / BRY[1] 582 BRY[3] = dmach[2] 583 NZ = 0 584 IFLAG = 0 585 KODED = KODE 586 RCAZ = 1.0e0 / CAZ 587 STR = ZR * RCAZ 588 STI = -ZI * RCAZ 589 RZR = (STR + STR) * RCAZ 590 RZI = (STI + STI) * RCAZ 591 INU = int(float32(FNU + 0.5)) 592 DNU = FNU - float64(INU) 593 if dabs(DNU) == 0.5e0 { 594 goto OneTen 595 } 596 DNU2 = 0.0e0 597 if dabs(DNU) > TOL { 598 DNU2 = DNU * DNU 599 } 600 if CAZ > R1 { 601 goto OneTen 602 } 603 604 // SERIES FOR CABS(Z)<=R1. 605 FC = 1.0e0 606 tmp = zlog(complex(RZR, RZI)) 607 SMUR = real(tmp) 608 SMUI = imag(tmp) 609 FMUR = SMUR * DNU 610 FMUI = SMUI * DNU 611 FMUR, FMUI, CSHR, CSHI, CCHR, CCHI = zshchOrig(FMUR, FMUI, CSHR, CSHI, CCHR, CCHI) 612 if DNU == 0.0e0 { 613 goto Ten 614 } 615 FC = DNU * DPI 616 FC = FC / dsin(FC) 617 SMUR = CSHR / DNU 618 SMUI = CSHI / DNU 619Ten: 620 A2 = 1.0e0 + DNU 621 622 // GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU). 623 T2 = dexp(-dgamln(A2, IDUM)) 624 T1 = 1.0e0 / (T2 * FC) 625 if dabs(DNU) > 0.1e0 { 626 goto Forty 627 } 628 629 // SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU). 630 AK = 1.0e0 631 S = CC[1] 632 for K = 2; K <= 8; K++ { 633 AK = AK * DNU2 634 TM = CC[K] * AK 635 S = S + TM 636 if dabs(TM) < TOL { 637 goto Thirty 638 } 639 } 640Thirty: 641 G1 = -S 642 goto Fifty 643Forty: 644 G1 = (T1 - T2) / (DNU + DNU) 645Fifty: 646 G2 = (T1 + T2) * 0.5e0 647 FR = FC * (CCHR*G1 + SMUR*G2) 648 FI = FC * (CCHI*G1 + SMUI*G2) 649 tmp = zexp(complex(FMUR, FMUI)) 650 STR = real(tmp) 651 STI = imag(tmp) 652 PR = 0.5e0 * STR / T2 653 PI = 0.5e0 * STI / T2 654 tmp = zdiv(complex(0.5, 0), complex(STR, STI)) 655 PTR = real(tmp) 656 PTI = imag(tmp) 657 QR = PTR / T1 658 QI = PTI / T1 659 S1R = FR 660 S1I = FI 661 S2R = PR 662 S2I = PI 663 AK = 1.0e0 664 A1 = 1.0e0 665 CKR = CONER 666 CKI = CONEI 667 BK = 1.0e0 - DNU2 668 if INU > 0 || N > 1 { 669 goto Eighty 670 } 671 672 // GENERATE K(FNU,Z), 0.0E0 <= FNU < 0.5E0 AND N=1. 673 if CAZ < TOL { 674 goto Seventy 675 } 676 tmp = zmlt(complex(ZR, ZI), complex(ZR, ZI)) 677 CZR = real(tmp) 678 CZI = imag(tmp) 679 CZR = 0.25e0 * CZR 680 CZI = 0.25e0 * CZI 681 T1 = 0.25e0 * CAZ * CAZ 682Sixty: 683 FR = (FR*AK + PR + QR) / BK 684 FI = (FI*AK + PI + QI) / BK 685 STR = 1.0e0 / (AK - DNU) 686 PR = PR * STR 687 PI = PI * STR 688 STR = 1.0e0 / (AK + DNU) 689 QR = QR * STR 690 QI = QI * STR 691 STR = CKR*CZR - CKI*CZI 692 RAK = 1.0e0 / AK 693 CKI = (CKR*CZI + CKI*CZR) * RAK 694 CKR = STR * RAK 695 S1R = CKR*FR - CKI*FI + S1R 696 S1I = CKR*FI + CKI*FR + S1I 697 A1 = A1 * T1 * RAK 698 BK = BK + AK + AK + 1.0e0 699 AK = AK + 1.0e0 700 if A1 > TOL { 701 goto Sixty 702 } 703Seventy: 704 YR[1] = S1R 705 YI[1] = S1I 706 if KODED == 1 { 707 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 708 } 709 tmp = zexp(complex(ZR, ZI)) 710 STR = real(tmp) 711 STI = imag(tmp) 712 tmp = zmlt(complex(S1R, S1I), complex(STR, STI)) 713 YR[1] = real(tmp) 714 YI[1] = imag(tmp) 715 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 716 717 // GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE. 718Eighty: 719 if CAZ < TOL { 720 goto OneHundred 721 } 722 tmp = zmlt(complex(ZR, ZI), complex(ZR, ZI)) 723 CZR = real(tmp) 724 CZI = imag(tmp) 725 CZR = 0.25e0 * CZR 726 CZI = 0.25e0 * CZI 727 T1 = 0.25e0 * CAZ * CAZ 728Ninety: 729 FR = (FR*AK + PR + QR) / BK 730 FI = (FI*AK + PI + QI) / BK 731 STR = 1.0e0 / (AK - DNU) 732 PR = PR * STR 733 PI = PI * STR 734 STR = 1.0e0 / (AK + DNU) 735 QR = QR * STR 736 QI = QI * STR 737 STR = CKR*CZR - CKI*CZI 738 RAK = 1.0e0 / AK 739 CKI = (CKR*CZI + CKI*CZR) * RAK 740 CKR = STR * RAK 741 S1R = CKR*FR - CKI*FI + S1R 742 S1I = CKR*FI + CKI*FR + S1I 743 STR = PR - FR*AK 744 STI = PI - FI*AK 745 S2R = CKR*STR - CKI*STI + S2R 746 S2I = CKR*STI + CKI*STR + S2I 747 A1 = A1 * T1 * RAK 748 BK = BK + AK + AK + 1.0e0 749 AK = AK + 1.0e0 750 if A1 > TOL { 751 goto Ninety 752 } 753OneHundred: 754 KFLAG = 2 755 A1 = FNU + 1.0e0 756 AK = A1 * dabs(SMUR) 757 if AK > ALIM { 758 KFLAG = 3 759 } 760 STR = CSSR[KFLAG] 761 P2R = S2R * STR 762 P2I = S2I * STR 763 tmp = zmlt(complex(P2R, P2I), complex(RZR, RZI)) 764 S2R = real(tmp) 765 S2I = imag(tmp) 766 S1R = S1R * STR 767 S1I = S1I * STR 768 if KODED == 1 { 769 goto TwoTen 770 } 771 tmp = zexp(complex(ZR, ZI)) 772 FR = real(tmp) 773 FI = imag(tmp) 774 tmp = zmlt(complex(S1R, S1I), complex(FR, FI)) 775 S1R = real(tmp) 776 S1I = imag(tmp) 777 tmp = zmlt(complex(S2R, S2I), complex(FR, FI)) 778 S2R = real(tmp) 779 S2I = imag(tmp) 780 goto TwoTen 781 782 // IFLAG=0 MEANS NO UNDERFLOW OCCURRED 783 // IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH 784 // KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD RECURSION 785OneTen: 786 tmp = zsqrt(complex(ZR, ZI)) 787 STR = real(tmp) 788 STI = imag(tmp) 789 tmp = zdiv(complex(RTHPI, CZEROI), complex(STR, STI)) 790 COEFR = real(tmp) 791 COEFI = imag(tmp) 792 KFLAG = 2 793 if KODED == 2 { 794 goto OneTwenty 795 } 796 if ZR > ALIM { 797 goto TwoNinety 798 } 799 800 STR = dexp(-ZR) * CSSR[KFLAG] 801 STI = -STR * dsin(ZI) 802 STR = STR * dcos(ZI) 803 tmp = zmlt(complex(COEFR, COEFI), complex(STR, STI)) 804 COEFR = real(tmp) 805 COEFI = imag(tmp) 806OneTwenty: 807 if dabs(DNU) == 0.5e0 { 808 goto ThreeHundred 809 } 810 // MILLER ALGORITHM FOR CABS(Z)>R1. 811 AK = dcos(DPI * DNU) 812 AK = dabs(AK) 813 if AK == CZEROR { 814 goto ThreeHundred 815 } 816 FHS = dabs(0.25e0 - DNU2) 817 if FHS == CZEROR { 818 goto ThreeHundred 819 } 820 821 // COMPUTE R2=F(E). if CABS(Z)>=R2, USE FORWARD RECURRENCE TO 822 // DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON 823 // 12<=E<=60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(14))= 824 // TOL WHERE B IS THE BASE OF THE ARITHMETIC. 825 T1 = float64(imach[14] - 1) 826 T1 = T1 * dmach[5] * 3.321928094e0 827 T1 = dmax(T1, 12.0e0) 828 T1 = dmin(T1, 60.0e0) 829 T2 = TTH*T1 - 6.0e0 830 if ZR != 0.0e0 { 831 goto OneThirty 832 } 833 T1 = HPI 834 goto OneFourty 835OneThirty: 836 T1 = datan(ZI / ZR) 837 T1 = dabs(T1) 838OneFourty: 839 if T2 > CAZ { 840 goto OneSeventy 841 } 842 // FORWARD RECURRENCE LOOP WHEN CABS(Z)>=R2. 843 ETEST = AK / (DPI * CAZ * TOL) 844 FK = CONER 845 if ETEST < CONER { 846 goto OneEighty 847 } 848 FKS = CTWOR 849 CKR = CAZ + CAZ + CTWOR 850 P1R = CZEROR 851 P2R = CONER 852 for I = 1; I <= KMAX; I++ { 853 AK = FHS / FKS 854 CBR = CKR / (FK + CONER) 855 PTR = P2R 856 P2R = CBR*P2R - P1R*AK 857 P1R = PTR 858 CKR = CKR + CTWOR 859 FKS = FKS + FK + FK + CTWOR 860 FHS = FHS + FK + FK 861 FK = FK + CONER 862 STR = dabs(P2R) * FK 863 if ETEST < STR { 864 goto OneSixty 865 } 866 } 867 goto ThreeTen 868OneSixty: 869 FK = FK + SPI*T1*dsqrt(T2/CAZ) 870 FHS = dabs(0.25 - DNU2) 871 goto OneEighty 872OneSeventy: 873 // COMPUTE BACKWARD INDEX K FOR CABS(Z)<R2. 874 A2 = dsqrt(CAZ) 875 AK = FPI * AK / (TOL * dsqrt(A2)) 876 AA = 3.0e0 * T1 / (1.0e0 + CAZ) 877 BB = 14.7e0 * T1 / (28.0e0 + CAZ) 878 AK = (dlog(AK) + CAZ*dcos(AA)/(1.0e0+0.008e0*CAZ)) / dcos(BB) 879 FK = 0.12125e0*AK*AK/CAZ + 1.5e0 880OneEighty: 881 // BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM. 882 K = int(float32(FK)) 883 FK = float64(K) 884 FKS = FK * FK 885 P1R = CZEROR 886 P1I = CZEROI 887 P2R = TOL 888 P2I = CZEROI 889 CSR = P2R 890 CSI = P2I 891 for I = 1; I <= K; I++ { 892 A1 = FKS - FK 893 AK = (FKS + FK) / (A1 + FHS) 894 RAK = 2.0e0 / (FK + CONER) 895 CBR = (FK + ZR) * RAK 896 CBI = ZI * RAK 897 PTR = P2R 898 PTI = P2I 899 P2R = (PTR*CBR - PTI*CBI - P1R) * AK 900 P2I = (PTI*CBR + PTR*CBI - P1I) * AK 901 P1R = PTR 902 P1I = PTI 903 CSR = CSR + P2R 904 CSI = CSI + P2I 905 FKS = A1 - FK + CONER 906 FK = FK - CONER 907 } 908 // COMPUTE (P2/CS)=(P2/CABS(CS))*(CONJG(CS)/CABS(CS)) FOR BETTER SCALING. 909 TM = zabs(complex(CSR, CSI)) 910 PTR = 1.0e0 / TM 911 S1R = P2R * PTR 912 S1I = P2I * PTR 913 CSR = CSR * PTR 914 CSI = -CSI * PTR 915 tmp = zmlt(complex(COEFR, COEFI), complex(S1R, S1I)) 916 STR = real(tmp) 917 STI = imag(tmp) 918 tmp = zmlt(complex(STR, STI), complex(CSR, CSI)) 919 S1R = real(tmp) 920 S1I = imag(tmp) 921 if INU > 0 || N > 1 { 922 goto TwoHundred 923 } 924 ZDR = ZR 925 ZDI = ZI 926 if IFLAG == 1 { 927 goto TwoSeventy 928 } 929 goto TwoFourty 930TwoHundred: 931 // COMPUTE P1/P2=(P1/CABS(P2)*CONJG(P2)/CABS(P2) FOR SCALING. 932 TM = zabs(complex(P2R, P2I)) 933 PTR = 1.0e0 / TM 934 P1R = P1R * PTR 935 P1I = P1I * PTR 936 P2R = P2R * PTR 937 P2I = -P2I * PTR 938 tmp = zmlt(complex(P1R, P1I), complex(P2R, P2I)) 939 PTR = real(tmp) 940 PTI = imag(tmp) 941 STR = DNU + 0.5e0 - PTR 942 STI = -PTI 943 tmp = zdiv(complex(STR, STI), complex(ZR, ZI)) 944 STR = real(tmp) 945 STI = imag(tmp) 946 STR = STR + 1.0e0 947 tmp = zmlt(complex(STR, STI), complex(S1R, S1I)) 948 S2R = real(tmp) 949 S2I = imag(tmp) 950 951 // FORWARD RECURSION ON THE THREE TERM RECURSION WITH RELATION WITH 952 // SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3 953TwoTen: 954 STR = DNU + 1.0e0 955 CKR = STR * RZR 956 CKI = STR * RZI 957 if N == 1 { 958 INU = INU - 1 959 } 960 if INU > 0 { 961 goto TwoTwenty 962 } 963 if N > 1 { 964 goto TwoFifteen 965 } 966 S1R = S2R 967 S1I = S2I 968TwoFifteen: 969 ZDR = ZR 970 ZDI = ZI 971 if IFLAG == 1 { 972 goto TwoSeventy 973 } 974 goto TwoFourty 975TwoTwenty: 976 INUB = 1 977 if IFLAG == 1 { 978 goto TwoSixtyOne 979 } 980TwoTwentyFive: 981 P1R = CSRR[KFLAG] 982 ASCLE = BRY[KFLAG] 983 for I = INUB; I <= INU; I++ { 984 STR = S2R 985 STI = S2I 986 S2R = CKR*STR - CKI*STI + S1R 987 S2I = CKR*STI + CKI*STR + S1I 988 S1R = STR 989 S1I = STI 990 CKR = CKR + RZR 991 CKI = CKI + RZI 992 if KFLAG >= 3 { 993 continue 994 } 995 P2R = S2R * P1R 996 P2I = S2I * P1R 997 STR = dabs(P2R) 998 STI = dabs(P2I) 999 P2M = dmax(STR, STI) 1000 if P2M <= ASCLE { 1001 continue 1002 } 1003 KFLAG = KFLAG + 1 1004 ASCLE = BRY[KFLAG] 1005 S1R = S1R * P1R 1006 S1I = S1I * P1R 1007 S2R = P2R 1008 S2I = P2I 1009 STR = CSSR[KFLAG] 1010 S1R = S1R * STR 1011 S1I = S1I * STR 1012 S2R = S2R * STR 1013 S2I = S2I * STR 1014 P1R = CSRR[KFLAG] 1015 } 1016 if N != 1 { 1017 goto TwoFourty 1018 } 1019 S1R = S2R 1020 S1I = S2I 1021TwoFourty: 1022 STR = CSRR[KFLAG] 1023 YR[1] = S1R * STR 1024 YI[1] = S1I * STR 1025 if N == 1 { 1026 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 1027 } 1028 YR[2] = S2R * STR 1029 YI[2] = S2I * STR 1030 if N == 2 { 1031 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 1032 } 1033 KK = 2 1034TwoFifty: 1035 KK = KK + 1 1036 if KK > N { 1037 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 1038 } 1039 P1R = CSRR[KFLAG] 1040 ASCLE = BRY[KFLAG] 1041 for I = KK; I <= N; I++ { 1042 P2R = S2R 1043 P2I = S2I 1044 S2R = CKR*P2R - CKI*P2I + S1R 1045 S2I = CKI*P2R + CKR*P2I + S1I 1046 S1R = P2R 1047 S1I = P2I 1048 CKR = CKR + RZR 1049 CKI = CKI + RZI 1050 P2R = S2R * P1R 1051 P2I = S2I * P1R 1052 YR[I] = P2R 1053 YI[I] = P2I 1054 if KFLAG >= 3 { 1055 continue 1056 } 1057 STR = dabs(P2R) 1058 STI = dabs(P2I) 1059 P2M = dmax(STR, STI) 1060 if P2M <= ASCLE { 1061 continue 1062 } 1063 KFLAG = KFLAG + 1 1064 ASCLE = BRY[KFLAG] 1065 S1R = S1R * P1R 1066 S1I = S1I * P1R 1067 S2R = P2R 1068 S2I = P2I 1069 STR = CSSR[KFLAG] 1070 S1R = S1R * STR 1071 S1I = S1I * STR 1072 S2R = S2R * STR 1073 S2I = S2I * STR 1074 P1R = CSRR[KFLAG] 1075 } 1076 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 1077 1078 // IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW. 1079TwoSixtyOne: 1080 HELIM = 0.5e0 * ELIM 1081 ELM = dexp(-ELIM) 1082 CELMR = ELM 1083 ASCLE = BRY[1] 1084 ZDR = ZR 1085 ZDI = ZI 1086 IC = -1 1087 J = 2 1088 for I = 1; I <= INU; I++ { 1089 STR = S2R 1090 STI = S2I 1091 S2R = STR*CKR - STI*CKI + S1R 1092 S2I = STI*CKR + STR*CKI + S1I 1093 S1R = STR 1094 S1I = STI 1095 CKR = CKR + RZR 1096 CKI = CKI + RZI 1097 AS = zabs(complex(S2R, S2I)) 1098 ALAS = dlog(AS) 1099 P2R = -ZDR + ALAS 1100 if P2R < (-ELIM) { 1101 goto TwoSixtyThree 1102 } 1103 tmp = zlog(complex(S2R, S2I)) 1104 STR = real(tmp) 1105 STI = imag(tmp) 1106 P2R = -ZDR + STR 1107 P2I = -ZDI + STI 1108 P2M = dexp(P2R) / TOL 1109 P1R = P2M * dcos(P2I) 1110 P1I = P2M * dsin(P2I) 1111 P1R, P1I, NW, ASCLE, TOL = zuchkOrig(P1R, P1I, NW, ASCLE, TOL) 1112 if NW != 0 { 1113 goto TwoSixtyThree 1114 } 1115 J = 3 - J 1116 CYR[J] = P1R 1117 CYI[J] = P1I 1118 if IC == (I - 1) { 1119 goto TwoSixtyFour 1120 } 1121 IC = I 1122 continue 1123 TwoSixtyThree: 1124 if ALAS < HELIM { 1125 continue 1126 } 1127 ZDR = ZDR - ELIM 1128 S1R = S1R * CELMR 1129 S1I = S1I * CELMR 1130 S2R = S2R * CELMR 1131 S2I = S2I * CELMR 1132 } 1133 if N != 1 { 1134 goto TwoSeventy 1135 } 1136 S1R = S2R 1137 S1I = S2I 1138 goto TwoSeventy 1139TwoSixtyFour: 1140 KFLAG = 1 1141 INUB = I + 1 1142 S2R = CYR[J] 1143 S2I = CYI[J] 1144 J = 3 - J 1145 S1R = CYR[J] 1146 S1I = CYI[J] 1147 if INUB <= INU { 1148 goto TwoTwentyFive 1149 } 1150 if N != 1 { 1151 goto TwoFourty 1152 } 1153 S1R = S2R 1154 S1I = S2I 1155 goto TwoFourty 1156TwoSeventy: 1157 YR[1] = S1R 1158 YI[1] = S1I 1159 if N == 1 { 1160 goto TwoEighty 1161 } 1162 YR[2] = S2R 1163 YI[2] = S2I 1164TwoEighty: 1165 ASCLE = BRY[1] 1166 ZDR, ZDI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM = zksclOrig(ZDR, ZDI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM) 1167 INU = N - NZ 1168 if INU <= 0 { 1169 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 1170 } 1171 KK = NZ + 1 1172 S1R = YR[KK] 1173 S1I = YI[KK] 1174 YR[KK] = S1R * CSRR[1] 1175 YI[KK] = S1I * CSRR[1] 1176 if INU == 1 { 1177 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 1178 } 1179 KK = NZ + 2 1180 S2R = YR[KK] 1181 S2I = YI[KK] 1182 YR[KK] = S2R * CSRR[1] 1183 YI[KK] = S2I * CSRR[1] 1184 if INU == 2 { 1185 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 1186 } 1187 T2 = FNU + float64(float32(KK-1)) 1188 CKR = T2 * RZR 1189 CKI = T2 * RZI 1190 KFLAG = 1 1191 goto TwoFifty 1192TwoNinety: 1193 1194 // SCALE BY dexp(Z), IFLAG = 1 CASES. 1195 1196 KODED = 2 1197 IFLAG = 1 1198 KFLAG = 2 1199 goto OneTwenty 1200 1201 // FNU=HALF ODD INTEGER CASE, DNU=-0.5 1202ThreeHundred: 1203 S1R = COEFR 1204 S1I = COEFI 1205 S2R = COEFR 1206 S2I = COEFI 1207 goto TwoTen 1208 1209ThreeTen: 1210 NZ = -2 1211 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 1212} 1213 1214// SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE 1215// ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN 1216// return WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL. 1217func zksclOrig(ZRR, ZRI, FNU float64, N int, YR, YI []float64, NZ int, RZR, RZI, ASCLE, TOL, ELIM float64) ( 1218 ZRRout, ZRIout, FNUout float64, Nout int, YRout, YIout []float64, NZout int, RZRout, RZIout, ASCLEout, TOLout, ELIMout float64) { 1219 var ACS, AS, CKI, CKR, CSI, CSR, FN, STR, S1I, S1R, S2I, 1220 S2R, ZEROI, ZEROR, ZDR, ZDI, CELMR, ELM, HELIM, ALAS float64 1221 1222 var I, IC, KK, NN, NW int 1223 var tmp complex128 1224 var CYR, CYI [3]float64 1225 // DIMENSION YR(N), YI(N), CYR(2), CYI(2) 1226 ZEROR = 0 1227 ZEROI = 0 1228 NZ = 0 1229 IC = 0 1230 NN = min0(2, N) 1231 for I = 1; I <= NN; I++ { 1232 S1R = YR[I] 1233 S1I = YI[I] 1234 CYR[I] = S1R 1235 CYI[I] = S1I 1236 AS = zabs(complex(S1R, S1I)) 1237 ACS = -ZRR + dlog(AS) 1238 NZ = NZ + 1 1239 YR[I] = ZEROR 1240 YI[I] = ZEROI 1241 if ACS < (-ELIM) { 1242 continue 1243 } 1244 1245 tmp = zlog(complex(S1R, S1I)) 1246 CSR = real(tmp) 1247 CSI = imag(tmp) 1248 CSR = CSR - ZRR 1249 CSI = CSI - ZRI 1250 STR = dexp(CSR) / TOL 1251 CSR = STR * dcos(CSI) 1252 CSI = STR * dsin(CSI) 1253 CSR, CSI, NW, ASCLE, TOL = zuchkOrig(CSR, CSI, NW, ASCLE, TOL) 1254 if NW != 0 { 1255 continue 1256 } 1257 YR[I] = CSR 1258 YI[I] = CSI 1259 IC = I 1260 NZ = NZ - 1 1261 } 1262 if N == 1 { 1263 return ZRR, ZRI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM 1264 } 1265 if IC > 1 { 1266 goto Twenty 1267 } 1268 YR[1] = ZEROR 1269 YI[1] = ZEROI 1270 NZ = 2 1271Twenty: 1272 if N == 2 { 1273 return ZRR, ZRI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM 1274 } 1275 if NZ == 0 { 1276 return ZRR, ZRI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM 1277 } 1278 FN = FNU + 1.0e0 1279 CKR = FN * RZR 1280 CKI = FN * RZI 1281 S1R = CYR[1] 1282 S1I = CYI[1] 1283 S2R = CYR[2] 1284 S2I = CYI[2] 1285 HELIM = 0.5e0 * ELIM 1286 ELM = dexp(-ELIM) 1287 CELMR = ELM 1288 ZDR = ZRR 1289 ZDI = ZRI 1290 1291 // FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF 1292 // S2 GETS LARGER THAN EXP(ELIM/2) 1293 for I = 3; I <= N; I++ { 1294 KK = I 1295 CSR = S2R 1296 CSI = S2I 1297 S2R = CKR*CSR - CKI*CSI + S1R 1298 S2I = CKI*CSR + CKR*CSI + S1I 1299 S1R = CSR 1300 S1I = CSI 1301 CKR = CKR + RZR 1302 CKI = CKI + RZI 1303 AS = zabs(complex(S2R, S2I)) 1304 ALAS = dlog(AS) 1305 ACS = -ZDR + ALAS 1306 NZ = NZ + 1 1307 YR[I] = ZEROR 1308 YI[I] = ZEROI 1309 if ACS < (-ELIM) { 1310 goto TwentyFive 1311 } 1312 tmp = zlog(complex(S2R, S2I)) 1313 CSR = real(tmp) 1314 CSI = imag(tmp) 1315 CSR = CSR - ZDR 1316 CSI = CSI - ZDI 1317 STR = dexp(CSR) / TOL 1318 CSR = STR * dcos(CSI) 1319 CSI = STR * dsin(CSI) 1320 CSR, CSI, NW, ASCLE, TOL = zuchkOrig(CSR, CSI, NW, ASCLE, TOL) 1321 if NW != 0 { 1322 goto TwentyFive 1323 } 1324 YR[I] = CSR 1325 YI[I] = CSI 1326 NZ = NZ - 1 1327 if IC == KK-1 { 1328 goto Forty 1329 } 1330 IC = KK 1331 continue 1332 TwentyFive: 1333 if ALAS < HELIM { 1334 continue 1335 } 1336 ZDR = ZDR - ELIM 1337 S1R = S1R * CELMR 1338 S1I = S1I * CELMR 1339 S2R = S2R * CELMR 1340 S2I = S2I * CELMR 1341 } 1342 NZ = N 1343 if IC == N { 1344 NZ = N - 1 1345 } 1346 goto FourtyFive 1347Forty: 1348 NZ = KK - 2 1349FourtyFive: 1350 for I = 1; I <= NZ; I++ { 1351 YR[I] = ZEROR 1352 YI[I] = ZEROI 1353 } 1354 return ZRR, ZRI, FNU, N, YR, YI, NZ, RZR, RZI, ASCLE, TOL, ELIM 1355} 1356 1357// Y ENTERS AS A SCALED QUANTITY WHOSE MAGNITUDE IS GREATER THAN 1358// EXP(-ALIM)=ASCLE=1.0E+3*dmach[1)/TOL. THE TEST IS MADE TO SEE 1359// if THE MAGNITUDE OF THE REAL OR IMAGINARY PART WOULD UNDERFLOW 1360// WHEN Y IS SCALED (BY TOL) TO ITS PROPER VALUE. Y IS ACCEPTED 1361// if THE UNDERFLOW IS AT LEAST ONE PRECISION BELOW THE MAGNITUDE 1362// OF THE LARGEST COMPONENT; OTHERWISE THE PHASE ANGLE DOES NOT HAVE 1363// ABSOLUTE ACCURACY AND AN UNDERFLOW IS ASSUMED. 1364func zuchkOrig(YR, YI float64, NZ int, ASCLE, TOL float64) (YRout, YIout float64, NZout int, ASCLEout, TOLout float64) { 1365 var SS, ST, WR, WI float64 1366 NZ = 0 1367 WR = dabs(YR) 1368 WI = dabs(YI) 1369 ST = dmin(WR, WI) 1370 if ST > ASCLE { 1371 return YR, YI, NZ, ASCLE, TOL 1372 } 1373 SS = dmax(WR, WI) 1374 ST = ST / TOL 1375 if SS < ST { 1376 NZ = 1 1377 } 1378 return YR, YI, NZ, ASCLE, TOL 1379} 1380 1381// ZACAI APPLIES THE ANALYTIC CONTINUATION FORMULA 1382// 1383// K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN) 1384// MP=PI*MR*CMPLX(0.0,1.0) 1385// 1386// TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT 1387// HALF Z PLANE FOR USE WITH ZAIRY WHERE FNU=1/3 OR 2/3 AND N=1. 1388// ZACAI IS THE SAME AS ZACON WITH THE PARTS FOR LARGER ORDERS AND 1389// RECURRENCE REMOVED. A RECURSIVE CALL TO ZACON CAN RESULT if ZACON 1390// IS CALLED FROM ZAIRY. 1391func zacaiOrig(ZR, ZI, FNU float64, KODE, MR, N int, YR, YI []float64, NZ int, RL, TOL, ELIM, ALIM float64) ( 1392 ZRout, ZIout, FNUout float64, KODEout, MRout, Nout int, YRout, YIout []float64, NZout int, RLout, TOLout, ELIMout, ALIMout float64) { 1393 var ARG, ASCLE, AZ, CSGNR, CSGNI, CSPNR, 1394 CSPNI, C1R, C1I, C2R, C2I, DFNU, FMR, PI, 1395 SGN, YY, ZNR, ZNI float64 1396 var INU, IUF, NN, NW int 1397 CYR := []float64{math.NaN(), 0, 0} 1398 CYI := []float64{math.NaN(), 0, 0} 1399 1400 PI = math.Pi 1401 NZ = 0 1402 ZNR = -ZR 1403 ZNI = -ZI 1404 AZ = zabs(complex(ZR, ZI)) 1405 NN = N 1406 DFNU = FNU + float64(float32(N-1)) 1407 if AZ <= 2.0e0 { 1408 goto Ten 1409 } 1410 if AZ*AZ*0.25 > DFNU+1.0e0 { 1411 goto Twenty 1412 } 1413Ten: 1414 // POWER SERIES FOR THE I FUNCTION. 1415 ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL, ELIM, ALIM = zseriOrig(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL, ELIM, ALIM) 1416 goto Forty 1417Twenty: 1418 if AZ < RL { 1419 goto Thirty 1420 } 1421 // ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION. 1422 ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM, ALIM = zasyiOrig(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, RL, TOL, ELIM, ALIM) 1423 if NW < 0 { 1424 goto Eighty 1425 } 1426 goto Forty 1427Thirty: 1428 // MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION 1429 ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL = zmlriOrig(ZNR, ZNI, FNU, KODE, NN, YR, YI, NW, TOL) 1430 if NW < 0 { 1431 goto Eighty 1432 } 1433Forty: 1434 // ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION. 1435 ZNR, ZNI, FNU, KODE, _, CYR, CYI, NW, TOL, ELIM, ALIM = zbknuOrig(ZNR, ZNI, FNU, KODE, 1, CYR, CYI, NW, TOL, ELIM, ALIM) 1436 if NW != 0 { 1437 goto Eighty 1438 } 1439 FMR = float64(float32(MR)) 1440 SGN = -math.Copysign(PI, FMR) 1441 CSGNR = 0.0e0 1442 CSGNI = SGN 1443 if KODE == 1 { 1444 goto Fifty 1445 } 1446 YY = -ZNI 1447 CSGNR = -CSGNI * dsin(YY) 1448 CSGNI = CSGNI * dcos(YY) 1449Fifty: 1450 // CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE 1451 // WHEN FNU IS LARGE 1452 INU = int(float32(FNU)) 1453 ARG = (FNU - float64(float32(INU))) * SGN 1454 CSPNR = dcos(ARG) 1455 CSPNI = dsin(ARG) 1456 if INU%2 == 0 { 1457 goto Sixty 1458 } 1459 CSPNR = -CSPNR 1460 CSPNI = -CSPNI 1461Sixty: 1462 C1R = CYR[1] 1463 C1I = CYI[1] 1464 C2R = YR[1] 1465 C2I = YI[1] 1466 if KODE == 1 { 1467 goto Seventy 1468 } 1469 IUF = 0 1470 ASCLE = 1.0e+3 * dmach[1] / TOL 1471 ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF = zs1s2Orig(ZNR, ZNI, C1R, C1I, C2R, C2I, NW, ASCLE, ALIM, IUF) 1472 NZ = NZ + NW 1473Seventy: 1474 YR[1] = CSPNR*C1R - CSPNI*C1I + CSGNR*C2R - CSGNI*C2I 1475 YI[1] = CSPNR*C1I + CSPNI*C1R + CSGNR*C2I + CSGNI*C2R 1476 return ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, TOL, ELIM, ALIM 1477Eighty: 1478 NZ = -1 1479 if NW == -2 { 1480 NZ = -2 1481 } 1482 return ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, TOL, ELIM, ALIM 1483} 1484 1485// ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z)>=0.0 BY 1486// MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE 1487// REGION CABS(Z)>MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL return. 1488// NZ<0 INDICATES AN OVERFLOW ON KODE=1. 1489func zasyiOrig(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, RL, TOL, ELIM, ALIM float64) ( 1490 ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZout int, RLout, TOLout, ELIMout, ALIMout float64) { 1491 var AA, AEZ, AK, AK1I, AK1R, ARG, ARM, ATOL, 1492 AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI, 1493 CZR, DFNU, DKI, DKR, DNU2, EZI, EZR, FDN, PI, P1I, 1494 P1R, RAZ, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I, 1495 S2R, TZI, TZR, ZEROI, ZEROR float64 1496 1497 var I, IB, IL, INU, J, JL, K, KODED, M, NN int 1498 var tmp complex128 1499 1500 PI = math.Pi 1501 RTPI = 0.159154943091895336e0 1502 ZEROR = 0 1503 ZEROI = 0 1504 CONER = 1 1505 CONEI = 0 1506 1507 NZ = 0 1508 AZ = zabs(complex(ZR, ZI)) 1509 ARM = 1.0e3 * dmach[1] 1510 RTR1 = dsqrt(ARM) 1511 IL = min0(2, N) 1512 DFNU = FNU + float64(float32(N-IL)) 1513 1514 // OVERFLOW TEST 1515 RAZ = 1.0e0 / AZ 1516 STR = ZR * RAZ 1517 STI = -ZI * RAZ 1518 AK1R = RTPI * STR * RAZ 1519 AK1I = RTPI * STI * RAZ 1520 tmp = zsqrt(complex(AK1R, AK1I)) 1521 AK1R = real(tmp) 1522 AK1I = imag(tmp) 1523 CZR = ZR 1524 CZI = ZI 1525 if KODE != 2 { 1526 goto Ten 1527 } 1528 CZR = ZEROR 1529 CZI = ZI 1530Ten: 1531 if dabs(CZR) > ELIM { 1532 goto OneHundred 1533 } 1534 DNU2 = DFNU + DFNU 1535 KODED = 1 1536 if (dabs(CZR) > ALIM) && (N > 2) { 1537 goto Twenty 1538 } 1539 KODED = 0 1540 tmp = zexp(complex(CZR, CZI)) 1541 STR = real(tmp) 1542 STI = imag(tmp) 1543 tmp = zmlt(complex(AK1R, AK1I), complex(STR, STI)) 1544 AK1R = real(tmp) 1545 AK1I = imag(tmp) 1546Twenty: 1547 FDN = 0.0e0 1548 if DNU2 > RTR1 { 1549 FDN = DNU2 * DNU2 1550 } 1551 EZR = ZR * 8.0e0 1552 EZI = ZI * 8.0e0 1553 1554 // WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE 1555 // FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE 1556 // EXPANSION FOR THE IMAGINARY PART. 1557 AEZ = 8.0e0 * AZ 1558 S = TOL / AEZ 1559 JL = int(float32(RL+RL)) + 2 1560 P1R = ZEROR 1561 P1I = ZEROI 1562 if ZI == 0.0e0 { 1563 goto Thirty 1564 } 1565 1566 // CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF 1567 // SIGNIFICANCE WHEN FNU OR N IS LARGE 1568 INU = int(float32(FNU)) 1569 ARG = (FNU - float64(float32(INU))) * PI 1570 INU = INU + N - IL 1571 AK = -dsin(ARG) 1572 BK = dcos(ARG) 1573 if ZI < 0.0e0 { 1574 BK = -BK 1575 } 1576 P1R = AK 1577 P1I = BK 1578 if INU%2 == 0 { 1579 goto Thirty 1580 } 1581 P1R = -P1R 1582 P1I = -P1I 1583Thirty: 1584 for K = 1; K <= IL; K++ { 1585 SQK = FDN - 1.0e0 1586 ATOL = S * dabs(SQK) 1587 SGN = 1.0e0 1588 CS1R = CONER 1589 CS1I = CONEI 1590 CS2R = CONER 1591 CS2I = CONEI 1592 CKR = CONER 1593 CKI = CONEI 1594 AK = 0.0e0 1595 AA = 1.0e0 1596 BB = AEZ 1597 DKR = EZR 1598 DKI = EZI 1599 // TODO(btracey): This loop is executed tens of thousands of times. Why? 1600 // is that really necessary? 1601 for J = 1; J <= JL; J++ { 1602 tmp = zdiv(complex(CKR, CKI), complex(DKR, DKI)) 1603 STR = real(tmp) 1604 STI = imag(tmp) 1605 CKR = STR * SQK 1606 CKI = STI * SQK 1607 CS2R = CS2R + CKR 1608 CS2I = CS2I + CKI 1609 SGN = -SGN 1610 CS1R = CS1R + CKR*SGN 1611 CS1I = CS1I + CKI*SGN 1612 DKR = DKR + EZR 1613 DKI = DKI + EZI 1614 AA = AA * dabs(SQK) / BB 1615 BB = BB + AEZ 1616 AK = AK + 8.0e0 1617 SQK = SQK - AK 1618 if AA <= ATOL { 1619 goto Fifty 1620 } 1621 } 1622 goto OneTen 1623 Fifty: 1624 S2R = CS1R 1625 S2I = CS1I 1626 if ZR+ZR >= ELIM { 1627 goto Sixty 1628 } 1629 TZR = ZR + ZR 1630 TZI = ZI + ZI 1631 tmp = zexp(complex(-TZR, -TZI)) 1632 STR = real(tmp) 1633 STI = imag(tmp) 1634 tmp = zmlt(complex(STR, STI), complex(P1R, P1I)) 1635 STR = real(tmp) 1636 STI = imag(tmp) 1637 tmp = zmlt(complex(STR, STI), complex(CS2R, CS2I)) 1638 STR = real(tmp) 1639 STI = imag(tmp) 1640 S2R = S2R + STR 1641 S2I = S2I + STI 1642 Sixty: 1643 FDN = FDN + 8.0e0*DFNU + 4.0e0 1644 P1R = -P1R 1645 P1I = -P1I 1646 M = N - IL + K 1647 YR[M] = S2R*AK1R - S2I*AK1I 1648 YI[M] = S2R*AK1I + S2I*AK1R 1649 } 1650 if N <= 2 { 1651 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, ALIM 1652 } 1653 NN = N 1654 K = NN - 2 1655 AK = float64(float32(K)) 1656 STR = ZR * RAZ 1657 STI = -ZI * RAZ 1658 RZR = (STR + STR) * RAZ 1659 RZI = (STI + STI) * RAZ 1660 IB = 3 1661 for I = IB; I <= NN; I++ { 1662 YR[K] = (AK+FNU)*(RZR*YR[K+1]-RZI*YI[K+1]) + YR[K+2] 1663 YI[K] = (AK+FNU)*(RZR*YI[K+1]+RZI*YR[K+1]) + YI[K+2] 1664 AK = AK - 1.0e0 1665 K = K - 1 1666 } 1667 if KODED == 0 { 1668 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, ALIM 1669 } 1670 tmp = zexp(complex(CZR, CZI)) 1671 CKR = real(tmp) 1672 CKI = imag(tmp) 1673 for I = 1; I <= NN; I++ { 1674 STR = YR[I]*CKR - YI[I]*CKI 1675 YI[I] = YR[I]*CKI + YI[I]*CKR 1676 YR[I] = STR 1677 } 1678 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, ALIM 1679OneHundred: 1680 NZ = -1 1681 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, ALIM 1682OneTen: 1683 NZ = -2 1684 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, ALIM 1685} 1686 1687// ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z)>=0.0 BY THE 1688// MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES. 1689func zmlriOrig(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL float64) ( 1690 ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZout int, TOLout float64) { 1691 var ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI, 1692 CNORMR, CONEI, CONER, FKAP, FKK, FLAM, FNF, PTI, PTR, P1I, 1693 P1R, P2I, P2R, RAZ, RHO, RHO2, RZI, RZR, SCLE, STI, STR, SUMI, 1694 SUMR, TFNF, TST, ZEROI, ZEROR float64 1695 var I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, M int 1696 var tmp complex128 1697 ZEROR = 0 1698 ZEROI = 0 1699 CONER = 1 1700 CONEI = 0 1701 1702 SCLE = dmach[1] / TOL 1703 NZ = 0 1704 AZ = zabs(complex(ZR, ZI)) 1705 IAZ = int(float32(AZ)) 1706 IFNU = int(float32(FNU)) 1707 INU = IFNU + N - 1 1708 AT = float64(float32(IAZ)) + 1.0e0 1709 RAZ = 1.0e0 / AZ 1710 STR = ZR * RAZ 1711 STI = -ZI * RAZ 1712 CKR = STR * AT * RAZ 1713 CKI = STI * AT * RAZ 1714 RZR = (STR + STR) * RAZ 1715 RZI = (STI + STI) * RAZ 1716 P1R = ZEROR 1717 P1I = ZEROI 1718 P2R = CONER 1719 P2I = CONEI 1720 ACK = (AT + 1.0e0) * RAZ 1721 RHO = ACK + dsqrt(ACK*ACK-1.0e0) 1722 RHO2 = RHO * RHO 1723 TST = (RHO2 + RHO2) / ((RHO2 - 1.0e0) * (RHO - 1.0e0)) 1724 TST = TST / TOL 1725 1726 // COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES. 1727 //fmt.Println("before loop", P2R, P2I, CKR, CKI, RZR, RZI, TST, AK) 1728 AK = AT 1729 for I = 1; I <= 80; I++ { 1730 PTR = P2R 1731 PTI = P2I 1732 P2R = P1R - (CKR*PTR - CKI*PTI) 1733 P2I = P1I - (CKI*PTR + CKR*PTI) 1734 P1R = PTR 1735 P1I = PTI 1736 CKR = CKR + RZR 1737 CKI = CKI + RZI 1738 AP = zabs(complex(P2R, P2I)) 1739 if AP > TST*AK*AK { 1740 goto Twenty 1741 } 1742 AK = AK + 1.0e0 1743 } 1744 goto OneTen 1745Twenty: 1746 I = I + 1 1747 K = 0 1748 if INU < IAZ { 1749 goto Forty 1750 } 1751 // COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS. 1752 P1R = ZEROR 1753 P1I = ZEROI 1754 P2R = CONER 1755 P2I = CONEI 1756 AT = float64(float32(INU)) + 1.0e0 1757 STR = ZR * RAZ 1758 STI = -ZI * RAZ 1759 CKR = STR * AT * RAZ 1760 CKI = STI * AT * RAZ 1761 ACK = AT * RAZ 1762 TST = dsqrt(ACK / TOL) 1763 ITIME = 1 1764 for K = 1; K <= 80; K++ { 1765 PTR = P2R 1766 PTI = P2I 1767 P2R = P1R - (CKR*PTR - CKI*PTI) 1768 P2I = P1I - (CKR*PTI + CKI*PTR) 1769 P1R = PTR 1770 P1I = PTI 1771 CKR = CKR + RZR 1772 CKI = CKI + RZI 1773 AP = zabs(complex(P2R, P2I)) 1774 if AP < TST { 1775 continue 1776 } 1777 if ITIME == 2 { 1778 goto Forty 1779 } 1780 ACK = zabs(complex(CKR, CKI)) 1781 FLAM = ACK + dsqrt(ACK*ACK-1.0e0) 1782 FKAP = AP / zabs(complex(P1R, P1I)) 1783 RHO = dmin(FLAM, FKAP) 1784 TST = TST * dsqrt(RHO/(RHO*RHO-1.0e0)) 1785 ITIME = 2 1786 } 1787 goto OneTen 1788Forty: 1789 // BACKWARD RECURRENCE AND SUM NORMALIZING RELATION. 1790 K = K + 1 1791 KK = max0(I+IAZ, K+INU) 1792 FKK = float64(float32(KK)) 1793 P1R = ZEROR 1794 P1I = ZEROI 1795 1796 // SCALE P2 AND SUM BY SCLE. 1797 P2R = SCLE 1798 P2I = ZEROI 1799 FNF = FNU - float64(float32(IFNU)) 1800 TFNF = FNF + FNF 1801 BK = dgamln(FKK+TFNF+1.0e0, IDUM) - dgamln(FKK+1.0e0, IDUM) - dgamln(TFNF+1.0e0, IDUM) 1802 BK = dexp(BK) 1803 SUMR = ZEROR 1804 SUMI = ZEROI 1805 KM = KK - INU 1806 for I = 1; I <= KM; I++ { 1807 PTR = P2R 1808 PTI = P2I 1809 P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) 1810 P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI) 1811 P1R = PTR 1812 P1I = PTI 1813 AK = 1.0e0 - TFNF/(FKK+TFNF) 1814 ACK = BK * AK 1815 SUMR = SUMR + (ACK+BK)*P1R 1816 SUMI = SUMI + (ACK+BK)*P1I 1817 BK = ACK 1818 FKK = FKK - 1.0e0 1819 } 1820 YR[N] = P2R 1821 YI[N] = P2I 1822 if N == 1 { 1823 goto Seventy 1824 } 1825 for I = 2; I <= N; I++ { 1826 PTR = P2R 1827 PTI = P2I 1828 P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) 1829 P2I = P1I + (FKK+FNF)*(RZI*PTR+RZR*PTI) 1830 P1R = PTR 1831 P1I = PTI 1832 AK = 1.0e0 - TFNF/(FKK+TFNF) 1833 ACK = BK * AK 1834 SUMR = SUMR + (ACK+BK)*P1R 1835 SUMI = SUMI + (ACK+BK)*P1I 1836 BK = ACK 1837 FKK = FKK - 1.0e0 1838 M = N - I + 1 1839 YR[M] = P2R 1840 YI[M] = P2I 1841 } 1842Seventy: 1843 if IFNU <= 0 { 1844 goto Ninety 1845 } 1846 for I = 1; I <= IFNU; I++ { 1847 PTR = P2R 1848 PTI = P2I 1849 P2R = P1R + (FKK+FNF)*(RZR*PTR-RZI*PTI) 1850 P2I = P1I + (FKK+FNF)*(RZR*PTI+RZI*PTR) 1851 P1R = PTR 1852 P1I = PTI 1853 AK = 1.0e0 - TFNF/(FKK+TFNF) 1854 ACK = BK * AK 1855 SUMR = SUMR + (ACK+BK)*P1R 1856 SUMI = SUMI + (ACK+BK)*P1I 1857 BK = ACK 1858 FKK = FKK - 1.0e0 1859 } 1860Ninety: 1861 PTR = ZR 1862 PTI = ZI 1863 if KODE == 2 { 1864 PTR = ZEROR 1865 } 1866 tmp = zlog(complex(RZR, RZI)) 1867 STR = real(tmp) 1868 STI = imag(tmp) 1869 P1R = -FNF*STR + PTR 1870 P1I = -FNF*STI + PTI 1871 AP = dgamln(1.0e0+FNF, IDUM) 1872 PTR = P1R - AP 1873 PTI = P1I 1874 1875 // THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW 1876 // IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES. 1877 P2R = P2R + SUMR 1878 P2I = P2I + SUMI 1879 AP = zabs(complex(P2R, P2I)) 1880 P1R = 1.0e0 / AP 1881 tmp = zexp(complex(PTR, PTI)) 1882 STR = real(tmp) 1883 STI = imag(tmp) 1884 CKR = STR * P1R 1885 CKI = STI * P1R 1886 PTR = P2R * P1R 1887 PTI = -P2I * P1R 1888 tmp = zmlt(complex(CKR, CKI), complex(PTR, PTI)) 1889 CNORMR = real(tmp) 1890 CNORMI = imag(tmp) 1891 for I = 1; I <= N; I++ { 1892 STR = YR[I]*CNORMR - YI[I]*CNORMI 1893 YI[I] = YR[I]*CNORMI + YI[I]*CNORMR 1894 YR[I] = STR 1895 } 1896 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL 1897OneTen: 1898 NZ = -2 1899 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL 1900} 1901 1902// ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z)>=0.0 BY 1903// MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE 1904// REGION CABS(Z)<=2*SQRT(FNU+1). NZ=0 IS A NORMAL return. 1905// NZ>0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO 1906// DUE TO UNDERFLOW. NZ<0 MEANS UNDERFLOW OCCURRED, BUT THE 1907// CONDITION CABS(Z)<=2*SQRT(FNU+1) WAS VIOLATED AND THE 1908// COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ). 1909func zseriOrig(ZR, ZI, FNU float64, KODE, N int, YR, YI []float64, NZ int, TOL, ELIM, ALIM float64) ( 1910 ZRout, ZIout, FNUout float64, KODEout, Nout int, YRout, YIout []float64, NZout int, TOLout, ELIMout, ALIMout float64) { 1911 var AA, ACZ, AK, AK1I, AK1R, ARM, ASCLE, ATOL, 1912 AZ, CKI, CKR, COEFI, COEFR, CONEI, CONER, CRSCR, CZI, CZR, DFNU, 1913 FNUP, HZI, HZR, RAZ, RS, RTR1, RZI, RZR, S, SS, STI, 1914 STR, S1I, S1R, S2I, S2R, ZEROI, ZEROR float64 1915 var I, IB, IDUM, IFLAG, IL, K, L, M, NN, NW int 1916 var WR, WI [3]float64 1917 var tmp complex128 1918 1919 CONER = 1.0 1920 NZ = 0 1921 AZ = zabs(complex(ZR, ZI)) 1922 if AZ == 0.0e0 { 1923 goto OneSixty 1924 } 1925 // TODO(btracey) 1926 // The original fortran line is "ARM = 1.0D+3*D1MACH(1)". Evidently, in Fortran 1927 // this is interpreted as one to the power of +3*D1MACH(1). While it is possible 1928 // this was intentional, it seems unlikely. 1929 //ARM = 1.0E0 + 3*dmach[1] 1930 //math.Pow(1, 3*dmach[1]) 1931 ARM = 1000 * dmach[1] 1932 RTR1 = dsqrt(ARM) 1933 CRSCR = 1.0e0 1934 IFLAG = 0 1935 if AZ < ARM { 1936 goto OneFifty 1937 } 1938 HZR = 0.5e0 * ZR 1939 HZI = 0.5e0 * ZI 1940 CZR = ZEROR 1941 CZI = ZEROI 1942 if AZ <= RTR1 { 1943 goto Ten 1944 } 1945 tmp = zmlt(complex(HZR, HZI), complex(HZR, HZI)) 1946 CZR = real(tmp) 1947 CZI = imag(tmp) 1948Ten: 1949 ACZ = zabs(complex(CZR, CZI)) 1950 NN = N 1951 tmp = zlog(complex(HZR, HZI)) 1952 CKR = real(tmp) 1953 CKI = imag(tmp) 1954Twenty: 1955 DFNU = FNU + float64(float32(NN-1)) 1956 FNUP = DFNU + 1.0e0 1957 1958 // UNDERFLOW TEST. 1959 AK1R = CKR * DFNU 1960 AK1I = CKI * DFNU 1961 AK = dgamln(FNUP, IDUM) 1962 AK1R = AK1R - AK 1963 if KODE == 2 { 1964 AK1R = AK1R - ZR 1965 } 1966 if AK1R > (-ELIM) { 1967 goto Forty 1968 } 1969Thirty: 1970 NZ = NZ + 1 1971 YR[NN] = ZEROR 1972 YI[NN] = ZEROI 1973 if ACZ > DFNU { 1974 goto OneNinety 1975 } 1976 NN = NN - 1 1977 if NN == 0 { 1978 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 1979 } 1980 goto Twenty 1981Forty: 1982 if AK1R > (-ALIM) { 1983 goto Fifty 1984 } 1985 IFLAG = 1 1986 SS = 1.0e0 / TOL 1987 CRSCR = TOL 1988 ASCLE = ARM * SS 1989Fifty: 1990 AA = dexp(AK1R) 1991 if IFLAG == 1 { 1992 AA = AA * SS 1993 } 1994 COEFR = AA * dcos(AK1I) 1995 COEFI = AA * dsin(AK1I) 1996 ATOL = TOL * ACZ / FNUP 1997 IL = min0(2, NN) 1998 for I = 1; I <= IL; I++ { 1999 DFNU = FNU + float64(float32(NN-I)) 2000 FNUP = DFNU + 1.0e0 2001 S1R = CONER 2002 S1I = CONEI 2003 if ACZ < TOL*FNUP { 2004 goto Seventy 2005 } 2006 AK1R = CONER 2007 AK1I = CONEI 2008 AK = FNUP + 2.0e0 2009 S = FNUP 2010 AA = 2.0e0 2011 Sixty: 2012 RS = 1.0e0 / S 2013 STR = AK1R*CZR - AK1I*CZI 2014 STI = AK1R*CZI + AK1I*CZR 2015 AK1R = STR * RS 2016 AK1I = STI * RS 2017 S1R = S1R + AK1R 2018 S1I = S1I + AK1I 2019 S = S + AK 2020 AK = AK + 2.0e0 2021 AA = AA * ACZ * RS 2022 if AA > ATOL { 2023 goto Sixty 2024 } 2025 Seventy: 2026 S2R = S1R*COEFR - S1I*COEFI 2027 S2I = S1R*COEFI + S1I*COEFR 2028 WR[I] = S2R 2029 WI[I] = S2I 2030 if IFLAG == 0 { 2031 goto Eighty 2032 } 2033 S2R, S2I, NW, ASCLE, TOL = zuchkOrig(S2R, S2I, NW, ASCLE, TOL) 2034 if NW != 0 { 2035 goto Thirty 2036 } 2037 Eighty: 2038 M = NN - I + 1 2039 YR[M] = S2R * CRSCR 2040 YI[M] = S2I * CRSCR 2041 if I == IL { 2042 continue 2043 } 2044 tmp = zdiv(complex(COEFR, COEFI), complex(HZR, HZI)) 2045 STR = real(tmp) 2046 STI = imag(tmp) 2047 COEFR = STR * DFNU 2048 COEFI = STI * DFNU 2049 } 2050 if NN <= 2 { 2051 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 2052 } 2053 K = NN - 2 2054 AK = float64(float32(K)) 2055 RAZ = 1.0e0 / AZ 2056 STR = ZR * RAZ 2057 STI = -ZI * RAZ 2058 RZR = (STR + STR) * RAZ 2059 RZI = (STI + STI) * RAZ 2060 if IFLAG == 1 { 2061 goto OneTwenty 2062 } 2063 IB = 3 2064OneHundred: 2065 for I = IB; I <= NN; I++ { 2066 YR[K] = (AK+FNU)*(RZR*YR[K+1]-RZI*YI[K+1]) + YR[K+2] 2067 YI[K] = (AK+FNU)*(RZR*YI[K+1]+RZI*YR[K+1]) + YI[K+2] 2068 AK = AK - 1.0e0 2069 K = K - 1 2070 } 2071 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 2072 2073 // RECUR BACKWARD WITH SCALED VALUES. 2074OneTwenty: 2075 // EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE 2076 // UNDERFLOW LIMIT = ASCLE = dmach[1)*SS*1.0D+3. 2077 S1R = WR[1] 2078 S1I = WI[1] 2079 S2R = WR[2] 2080 S2I = WI[2] 2081 for L = 3; L <= NN; L++ { 2082 CKR = S2R 2083 CKI = S2I 2084 S2R = S1R + (AK+FNU)*(RZR*CKR-RZI*CKI) 2085 S2I = S1I + (AK+FNU)*(RZR*CKI+RZI*CKR) 2086 S1R = CKR 2087 S1I = CKI 2088 CKR = S2R * CRSCR 2089 CKI = S2I * CRSCR 2090 YR[K] = CKR 2091 YI[K] = CKI 2092 AK = AK - 1.0e0 2093 K = K - 1 2094 if zabs(complex(CKR, CKI)) > ASCLE { 2095 goto OneFourty 2096 } 2097 } 2098 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 2099OneFourty: 2100 IB = L + 1 2101 if IB > NN { 2102 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 2103 } 2104 goto OneHundred 2105OneFifty: 2106 NZ = N 2107 if FNU == 0.0e0 { 2108 NZ = NZ - 1 2109 } 2110OneSixty: 2111 YR[1] = ZEROR 2112 YI[1] = ZEROI 2113 if FNU != 0.0e0 { 2114 goto OneSeventy 2115 } 2116 YR[1] = CONER 2117 YI[1] = CONEI 2118OneSeventy: 2119 if N == 1 { 2120 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 2121 } 2122 for I = 2; I <= N; I++ { 2123 YR[I] = ZEROR 2124 YI[I] = ZEROI 2125 } 2126 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 2127 2128 // return WITH NZ<0 if CABS(Z*Z/4)>FNU+N-NZ-1 COMPLETE 2129 // THE CALCULATION IN CBINU WITH N=N-IABS(NZ) 2130 2131OneNinety: 2132 NZ = -NZ 2133 return ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM 2134} 2135 2136// ZS1S2 TESTS FOR A POSSIBLE UNDERFLOW RESULTING FROM THE 2137// ADDITION OF THE I AND K FUNCTIONS IN THE ANALYTIC CON- 2138// TINUATION FORMULA WHERE S1=K FUNCTION AND S2=I FUNCTION. 2139// ON KODE=1 THE I AND K FUNCTIONS ARE DIFFERENT ORDERS OF 2140// MAGNITUDE, BUT FOR KODE=2 THEY CAN BE OF THE SAME ORDER 2141// OF MAGNITUDE AND THE MAXIMUM MUST BE AT LEAST ONE 2142// PRECISION ABOVE THE UNDERFLOW LIMIT. 2143func zs1s2Orig(ZRR, ZRI, S1R, S1I, S2R, S2I float64, NZ int, ASCLE, ALIM float64, IUF int) ( 2144 ZRRout, ZRIout, S1Rout, S1Iout, S2Rout, S2Iout float64, NZout int, ASCLEout, ALIMout float64, IUFout int) { 2145 var AA, ALN, AS1, AS2, C1I, C1R, S1DI, S1DR, ZEROI, ZEROR float64 2146 var tmp complex128 2147 2148 ZEROR = 0 2149 ZEROI = 0 2150 NZ = 0 2151 AS1 = zabs(complex(S1R, S1I)) 2152 AS2 = zabs(complex(S2R, S2I)) 2153 if S1R == 0.0e0 && S1I == 0.0e0 { 2154 goto Ten 2155 } 2156 if AS1 == 0.0e0 { 2157 goto Ten 2158 } 2159 ALN = -ZRR - ZRR + dlog(AS1) 2160 S1DR = S1R 2161 S1DI = S1I 2162 S1R = ZEROR 2163 S1I = ZEROI 2164 AS1 = ZEROR 2165 if ALN < (-ALIM) { 2166 goto Ten 2167 } 2168 tmp = zlog(complex(S1DR, S1DI)) 2169 C1R = real(tmp) 2170 C1I = imag(tmp) 2171 2172 C1R = C1R - ZRR - ZRR 2173 C1I = C1I - ZRI - ZRI 2174 tmp = zexp(complex(C1R, C1I)) 2175 S1R = real(tmp) 2176 S1I = imag(tmp) 2177 AS1 = zabs(complex(S1R, S1I)) 2178 IUF = IUF + 1 2179Ten: 2180 AA = dmax(AS1, AS2) 2181 if AA > ASCLE { 2182 return ZRR, ZRI, S1R, S1I, S2R, S2I, NZ, ASCLE, ALIM, IUF 2183 } 2184 S1R = ZEROR 2185 S1I = ZEROI 2186 S2R = ZEROR 2187 S2I = ZEROI 2188 NZ = 1 2189 IUF = 0 2190 return ZRR, ZRI, S1R, S1I, S2R, S2I, NZ, ASCLE, ALIM, IUF 2191} 2192 2193// ZSHCH COMPUTES THE COMPLEX HYPERBOLIC FUNCTIONS CSH=SINH(X+iY) AND 2194// CCH=COSH(X+I*Y), WHERE I**2=-1. 2195// TODO(btracey): use cmplx.Sinh and cmplx.Cosh. 2196func zshchOrig(ZR, ZI, CSHR, CSHI, CCHR, CCHI float64) (ZRout, ZIout, CSHRout, CSHIout, CCHRout, CCHIout float64) { 2197 var CH, CN, SH, SN float64 2198 SH = math.Sinh(ZR) 2199 CH = math.Cosh(ZR) 2200 SN = dsin(ZI) 2201 CN = dcos(ZI) 2202 CSHR = SH * CN 2203 CSHI = CH * SN 2204 CCHR = CH * CN 2205 CCHI = SH * SN 2206 return ZR, ZI, CSHR, CSHI, CCHR, CCHI 2207} 2208 2209func dmax(a, b float64) float64 { 2210 return math.Max(a, b) 2211} 2212 2213func dmin(a, b float64) float64 { 2214 return math.Min(a, b) 2215} 2216 2217func dabs(a float64) float64 { 2218 return math.Abs(a) 2219} 2220 2221func datan(a float64) float64 { 2222 return math.Atan(a) 2223} 2224 2225func dtan(a float64) float64 { 2226 return math.Tan(a) 2227} 2228 2229func dlog(a float64) float64 { 2230 return math.Log(a) 2231} 2232 2233func dsin(a float64) float64 { 2234 return math.Sin(a) 2235} 2236 2237func dcos(a float64) float64 { 2238 return math.Cos(a) 2239} 2240 2241func dexp(a float64) float64 { 2242 return math.Exp(a) 2243} 2244 2245func dsqrt(a float64) float64 { 2246 return math.Sqrt(a) 2247} 2248 2249func zmlt(a, b complex128) complex128 { 2250 return a * b 2251} 2252 2253func zdiv(a, b complex128) complex128 { 2254 return a / b 2255} 2256 2257func zabs(a complex128) float64 { 2258 return cmplx.Abs(a) 2259} 2260 2261func zsqrt(a complex128) complex128 { 2262 return cmplx.Sqrt(a) 2263} 2264 2265func zexp(a complex128) complex128 { 2266 return cmplx.Exp(a) 2267} 2268 2269func zlog(a complex128) complex128 { 2270 return cmplx.Log(a) 2271} 2272 2273// Zshch computes the hyperbolic sin and cosine of the input z. 2274func Zshch(z complex128) (sinh, cosh complex128) { 2275 return cmplx.Sinh(z), cmplx.Cosh(z) 2276} 2277