zseri.go

func Zseri(z complex128, fnu float64, kode, n int, y []complex128, tol, elim, alim float64) (nz int) {
	var AA, AK, ATOL,
		DFNU, FNUP, RS, S, SS float64
	var I, IB, IDUM, IFLAG, IL, L, M, NN, NW int
	var w [2]complex128
	var sin, cos float64
	var first bool

	var scale float64
	var ak1, coef, rz, s1, s2, st complex128

	// TOOD(btracey): The original fortran line is "ARM = 1.0D+3*D1MACH(1)". Evidently, in Fortran
	// this is interpreted as one to the power of +3*D1MACH(1). While it is possible
	// this was intentional, it seems unlikely.
	arm := 1000 * dmach[1]
	az := cmplx.Abs(z)
	if az < arm {
		for i := 0; i < n; i++ {
			y[i] = 0
		}
		if fnu == 0 {
			y[0] = 1
			n--
		}
		if az == 0 {
			return 0
		}
		return n
	}
	crscr := 1.0
	IFLAG = 0
	hz := 0.5 * z
	var cz complex128
	var ACZ float64
	if az > math.Sqrt(arm) {
		cz = hz * hz
		ACZ = cmplx.Abs(cz)
	}
	NN = n
	ck := cmplx.Log(hz)
	rotate := false
retry:
	for {
		if !rotate {
			DFNU = fnu + float64(NN-1)
			FNUP = DFNU + 1.0E0

			// UNDERFLOW TEST.
			ak1 = ck * complex(DFNU, 0)
			AK = dgamln(FNUP, IDUM)
			ak1 -= complex(AK, 0)
			if kode == 2 {
				ak1 -= complex(real(z), 0)
			}
			if real(ak1) > -elim {
				break
			}
		}
		rotate = false
		nz++
		y[NN-1] = 0
		if ACZ > DFNU {
			// Return with nz < 0 if abs(Z*Z/4)>fnu+u-nz-1 complete the calculation
			// in cbinu with n = n - abs(nz).
			nz *= -1
			return nz
		}
		NN--
		if NN == 0 {
			return nz
		}
	}
	if real(ak1) <= -alim {
		IFLAG = 1
		SS = 1 / tol
		crscr = tol
		scale = arm * SS
	}
	AA = math.Exp(real(ak1))
	if IFLAG == 1 {
		AA *= SS
	}
	sin, cos = math.Sincos(imag(ak1))
	coef = complex(AA*cos, AA*sin)
	ATOL = tol * ACZ / FNUP
	IL = min(2, NN)
	for I = 1; I <= IL; I++ {
		DFNU = fnu + float64(NN-I)
		FNUP = DFNU + 1.0E0
		s1 = 1
		if ACZ >= tol*FNUP {
			ak1 = 1
			AK = FNUP + 2.0E0
			S = FNUP
			AA = 2.0E0
			first = true
			for first || AA > ATOL {
				RS = 1.0E0 / S
				st = ak1 * cz
				ak1 = st * complex(RS, 0)
				s1 += ak1
				S += AK
				AK += 2
				AA *= ACZ * RS
				first = false
			}
		}
		s2 = s1 * coef
		w[I-1] = s2
		if IFLAG != 0 {
			NW = Zuchk(s2, scale, tol)
			if NW != 0 {
				rotate = true
				goto retry
			}
		}
		M = NN - I + 1
		y[M-1] = s2 * complex(crscr, 0)
		if I == IL {
			continue
		}
		st = coef / hz
		coef = st * complex(DFNU, 0)
	}
	if NN <= 2 {
		return nz
	}
	rz = complex(2*real(z)/(az*az), -2*imag(z)/(az*az))
	if IFLAG != 1 {
		for i := NN - 3; i >= 0; i-- {
			y[i] = complex(float64(i+1)+fnu, 0)*rz*y[i+1] + y[i+2]
		}
		return nz
	}
	// EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
	// UNDERFLOW LIMIT = ASCLE = dmach[1)*SS*1.0D+3.
	K := NN - 2
	AK = float64(K)
	s1 = w[0]
	s2 = w[1]
	for L = 3; L <= NN; L++ {
		s1, s2 = s2, s1+complex(AK+fnu, 0)*(rz*s2)
		ck := s2 * complex(crscr, 0)
		y[K-1] = ck
		AK--
		K--
		if cmplx.Abs(ck) > scale {
			IB = L + 1
			for I = IB; I <= NN; I++ {
				y[K-1] = complex(AK+fnu, 0)*rz*y[K] + y[K+1]
				AK--
				K--
			}
			break
		}
	}
	return nz
}