Partial implementation of Bessel phase function

["管理员" <c2976501934@qq.com>]

GEN besselphase(GEN n, GEN z, long prec) {
	pari_sp av = get_avma(), av2;
	GEN om, lim, est, phi, J, Y;
	long l = precision(z), s;

	if (l) prec = l;

	n = gtofp(n, prec);
	if (typ(n) != t_REAL) pari_err_IMPL("besselphase for complex order");
	z = gtofp(z, prec);
	if (typ(z) != t_REAL || signe(z) < 0) pari_err_IMPL("besselphase for nonpositive argument");
	om = Pi2n(1, prec);

	s = signe(n);
	if (s < 0) n = mpabs(n);

	if (!s || !signe(z) || expo(n) < -1) {
		est = subrr(z, shiftr(om, -2));
	} else {
		int cond;
		lim = n;
		est = gen_0;
		cond = cmprr(z, lim) > 0;
		av2 = get_avma();
		while (cond) {
			GEN dz, dphi, dphim1, lim2;
			J = real_i(jbessel(n, lim, prec));
			Y = real_i(ybessel(n, lim, prec));
			phi = garg(mkcomplex(J, Y), prec);
			phi = addrr(phi, mulri(om, roundr(divrr(gsub(est, phi), om))));
			dphi = divsr(4, mulrr(mulrr(om, z), addrr(sqrr(J), sqrr(Y))));
			dphim1 = subsr(1, dphi);
			dz = subrr(z, lim);
			cond = 0;
			if (signe(dphim1) > 0) {
				lim2 = shiftr(divrr(om, dphim1), -1);
				cond = cmprr(dz, lim2) > 0;
				if (cond) dz = lim2;
				dz = mulrr(dz, dphi);
				lim = addrr(lim, lim2);
			}
			est = addrr(phi, dz);
			if (gc_needed(av, 2)) {
				if (DEBUGMEM>1) pari_warn(warnmem, "besselphase");
				gerepileall(av2, 2, &est, &lim);
			}
		}
	}

	if (signe(z)) {
		J = real_i(jbessel(n, z, prec));
		Y = real_i(ybessel(n, z, prec));
		phi = garg(mkcomplex(J, Y), prec);
		phi = addrr(phi, mulri(om, roundr(divrr(gsub(est, phi), om))));
	} else {
		phi = est;
	}

	if (s < 0) {
		shiftr_inplace(om, -1);
		phi = addrr(phi, mulrr(om, n));
	}

	return gerepilecopy(av, phi);
}