mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
complex Lambert W function (acb_lambertw)
This commit is contained in:
Fredrik Johansson
parent
b30933ace9
commit
d77b6c5c73
7 changed files with 1121 additions and 0 deletions
12
acb.h
12
acb.h
|
|
@ -391,6 +391,13 @@ acb_get_rad_ubound_arf(arf_t u, const acb_t z, slong prec)
|
||||||
arf_mul_2exp_si(u, u, 1);
|
arf_mul_2exp_si(u, u, 1);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
ACB_INLINE void
|
||||||
|
acb_union(acb_t res, const acb_t x, const acb_t y, slong prec)
|
||||||
|
{
|
||||||
|
arb_union(acb_realref(res), acb_realref(x), acb_realref(y), prec);
|
||||||
|
arb_union(acb_imagref(res), acb_imagref(x), acb_imagref(y), prec);
|
||||||
|
}
|
||||||
|
|
||||||
void acb_arg(arb_t r, const acb_t z, slong prec);
|
void acb_arg(arb_t r, const acb_t z, slong prec);
|
||||||
|
|
||||||
void acb_sgn(acb_t res, const acb_t z, slong prec);
|
void acb_sgn(acb_t res, const acb_t z, slong prec);
|
||||||
|
|
@ -782,6 +789,11 @@ void acb_polylog_si(acb_t w, slong s, const acb_t z, slong prec);
|
||||||
void acb_agm1(acb_t m, const acb_t z, slong prec);
|
void acb_agm1(acb_t m, const acb_t z, slong prec);
|
||||||
void acb_agm1_cpx(acb_ptr m, const acb_t z, slong len, slong prec);
|
void acb_agm1_cpx(acb_ptr m, const acb_t z, slong len, slong prec);
|
||||||
|
|
||||||
|
void acb_lambertw_asymp(acb_t res, const acb_t z, const fmpz_t k, slong L, slong M, slong prec);
|
||||||
|
int acb_lambertw_check_branch(const acb_t w, const fmpz_t k, slong prec);
|
||||||
|
void acb_lambertw_bound_deriv(mag_t res, const acb_t z, const acb_t ez1, const fmpz_t k);
|
||||||
|
void acb_lambertw(acb_t res, const acb_t z, const fmpz_t k, int flags, slong prec);
|
||||||
|
|
||||||
ACB_INLINE void
|
ACB_INLINE void
|
||||||
acb_sqr(acb_t res, const acb_t val, slong prec)
|
acb_sqr(acb_t res, const acb_t val, slong prec)
|
||||||
{
|
{
|
||||||
|
|
|
||||||
593
acb/lambertw.c
Normal file
593
acb/lambertw.c
Normal file
|
|
@ -0,0 +1,593 @@
|
||||||
|
/*
|
||||||
|
Copyright (C) 2017 Fredrik Johansson
|
||||||
|
|
||||||
|
This file is part of Arb.
|
||||||
|
|
||||||
|
Arb is free software: you can redistribute it and/or modify it under
|
||||||
|
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||||
|
by the Free Software Foundation; either version 2.1 of the License, or
|
||||||
|
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#include "acb.h"
|
||||||
|
|
||||||
|
/* Check if z crosses a branch cut. */
|
||||||
|
int
|
||||||
|
acb_lambertw_branch_crossing(const acb_t z, const acb_t ez1, const fmpz_t k)
|
||||||
|
{
|
||||||
|
if (arb_contains_zero(acb_imagref(z)) && !arb_is_nonnegative(acb_imagref(z)))
|
||||||
|
{
|
||||||
|
if (fmpz_is_zero(k))
|
||||||
|
{
|
||||||
|
if (!arb_is_positive(acb_realref(ez1)))
|
||||||
|
{
|
||||||
|
return 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
else if (!arb_is_positive(acb_realref(z)))
|
||||||
|
{
|
||||||
|
return 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* todo: remove radii */
|
||||||
|
void
|
||||||
|
acb_lambertw_halley_step(acb_t res, const acb_t z, const acb_t w, slong prec)
|
||||||
|
{
|
||||||
|
acb_t ew, t, u, v;
|
||||||
|
|
||||||
|
acb_init(ew);
|
||||||
|
acb_init(t);
|
||||||
|
acb_init(u);
|
||||||
|
acb_init(v);
|
||||||
|
|
||||||
|
acb_exp(ew, w, prec);
|
||||||
|
acb_add_ui(u, w, 2, prec);
|
||||||
|
acb_add_ui(v, w, 1, prec);
|
||||||
|
acb_mul_2exp_si(v, v, 1);
|
||||||
|
acb_div(v, u, v, prec);
|
||||||
|
acb_mul(t, ew, w, prec);
|
||||||
|
acb_sub(u, t, z, prec);
|
||||||
|
acb_mul(v, v, u, prec);
|
||||||
|
acb_neg(v, v);
|
||||||
|
acb_add(v, v, t, prec);
|
||||||
|
acb_add(v, v, ew, prec);
|
||||||
|
acb_div(t, u, v, prec);
|
||||||
|
|
||||||
|
acb_sub(t, w, t, prec);
|
||||||
|
|
||||||
|
acb_swap(res, t);
|
||||||
|
|
||||||
|
acb_clear(ew);
|
||||||
|
acb_clear(t);
|
||||||
|
acb_clear(u);
|
||||||
|
acb_clear(v);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* assumes no aliasing of w and p */
|
||||||
|
void
|
||||||
|
acb_lambertw_branchpoint_series(acb_t w, const acb_t t, int bound, slong prec)
|
||||||
|
{
|
||||||
|
slong i;
|
||||||
|
static const int coeffs[] = {-130636800,130636800,-43545600,19958400,
|
||||||
|
-10402560,5813640,-3394560,2042589,-1256320};
|
||||||
|
|
||||||
|
acb_zero(w);
|
||||||
|
|
||||||
|
for (i = 8; i >= 0; i--)
|
||||||
|
{
|
||||||
|
acb_mul(w, w, t, prec);
|
||||||
|
acb_add_si(w, w, coeffs[i], prec);
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_div_si(w, w, -coeffs[0], prec);
|
||||||
|
|
||||||
|
if (bound)
|
||||||
|
{
|
||||||
|
mag_t err;
|
||||||
|
mag_init(err);
|
||||||
|
acb_get_mag(err, t);
|
||||||
|
mag_geom_series(err, err, 9);
|
||||||
|
|
||||||
|
if (acb_is_real(t))
|
||||||
|
arb_add_error_mag(acb_realref(w), err);
|
||||||
|
else
|
||||||
|
acb_add_error_mag(w, err);
|
||||||
|
mag_clear(err);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void
|
||||||
|
acb_lambertw_principal_d(acb_t res, const acb_t z)
|
||||||
|
{
|
||||||
|
double za, zb, wa, wb, ewa, ewb, t, u, q, r;
|
||||||
|
int k, maxk = 15;
|
||||||
|
|
||||||
|
za = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
|
||||||
|
zb = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
|
||||||
|
|
||||||
|
/* make sure we end up on the right branch */
|
||||||
|
if (za < -0.367 && zb > -1e-20 && zb <= 0.0
|
||||||
|
&& arf_sgn(arb_midref(acb_imagref(z))) < 0)
|
||||||
|
zb = -1e-20;
|
||||||
|
|
||||||
|
wa = za;
|
||||||
|
wb = zb;
|
||||||
|
|
||||||
|
if (fabs(wa) > 2.0 || fabs(wb) > 2.0)
|
||||||
|
{
|
||||||
|
t = atan2(wb, wa);
|
||||||
|
wa = 0.5 * log(wa * wa + wb * wb);
|
||||||
|
wb = t;
|
||||||
|
}
|
||||||
|
else if (fabs(wa) > 0.25 || fabs(wb) > 0.25)
|
||||||
|
{
|
||||||
|
/* We have W(z) ~= -1 + (2(ez+1))^(1/2) near the branch point.
|
||||||
|
Changing the exponent to 1/4 gives a much worse local guess
|
||||||
|
which however does the job on a larger domain. */
|
||||||
|
wa *= 5.43656365691809;
|
||||||
|
wb *= 5.43656365691809;
|
||||||
|
wa += 2.0;
|
||||||
|
t = atan2(wb, wa);
|
||||||
|
r = pow(wa * wa + wb * wb, 0.125);
|
||||||
|
wa = r * cos(0.25 * t);
|
||||||
|
wb = r * sin(0.25 * t);
|
||||||
|
wa -= 1.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
for (k = 0; k < maxk; k++)
|
||||||
|
{
|
||||||
|
t = exp(wa);
|
||||||
|
ewa = t * cos(wb);
|
||||||
|
ewb = t * sin(wb);
|
||||||
|
t = (ewa * wa - ewb * wb); q = t + ewa; t -= za;
|
||||||
|
u = (ewb * wa + ewa * wb); r = u + ewb; u -= zb;
|
||||||
|
ewa = q * t + r * u; ewb = q * u - r * t;
|
||||||
|
r = 1.0 / (q * q + r * r);
|
||||||
|
ewa *= r; ewb *= r;
|
||||||
|
if ((ewa*ewa + ewb*ewb) < (wa*wa + wb*wb) * 1e-12)
|
||||||
|
maxk = FLINT_MIN(maxk, k + 2);
|
||||||
|
wa -= ewa; wb -= ewb;
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_set_d_d(res, wa, wb);
|
||||||
|
}
|
||||||
|
|
||||||
|
void
|
||||||
|
acb_lambertw_initial_asymp(acb_t w, const acb_t z, const fmpz_t k, slong prec)
|
||||||
|
{
|
||||||
|
acb_t L1, L2, t;
|
||||||
|
|
||||||
|
acb_init(L1);
|
||||||
|
acb_init(L2);
|
||||||
|
acb_init(t);
|
||||||
|
|
||||||
|
acb_const_pi(L2, prec);
|
||||||
|
acb_mul_2exp_si(L2, L2, 1);
|
||||||
|
acb_mul_fmpz(L2, L2, k, prec);
|
||||||
|
acb_mul_onei(L2, L2);
|
||||||
|
acb_log(L1, z, prec);
|
||||||
|
acb_add(L1, L1, L2, prec);
|
||||||
|
acb_log(L2, L1, prec);
|
||||||
|
|
||||||
|
/* L1 - L2 + L2/L1 + L2(L2-2)/(2 L1^2) */
|
||||||
|
acb_inv(t, L1, prec);
|
||||||
|
acb_mul_2exp_si(w, L2, 1);
|
||||||
|
acb_submul(w, L2, L2, prec);
|
||||||
|
acb_neg(w, w);
|
||||||
|
acb_mul(w, w, t, prec);
|
||||||
|
acb_mul_2exp_si(w, w, -1);
|
||||||
|
acb_add(w, w, L2, prec);
|
||||||
|
acb_mul(w, w, t, prec);
|
||||||
|
acb_sub(w, w, L2, prec);
|
||||||
|
acb_add(w, w, L1, prec);
|
||||||
|
|
||||||
|
acb_clear(L1);
|
||||||
|
acb_clear(L2);
|
||||||
|
acb_clear(t);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* assumes no aliasing */
|
||||||
|
slong
|
||||||
|
acb_lambertw_initial(acb_t res, const acb_t z, const acb_t ez1, const fmpz_t k, slong prec)
|
||||||
|
{
|
||||||
|
/* Handle z very close to 0 on the principal branch. */
|
||||||
|
if (fmpz_is_zero(k) &&
|
||||||
|
(arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), -20) <= 0 &&
|
||||||
|
arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), -20) <= 0))
|
||||||
|
{
|
||||||
|
acb_set(res, z);
|
||||||
|
acb_submul(res, res, res, prec);
|
||||||
|
return 40; /* could be tightened... */
|
||||||
|
}
|
||||||
|
|
||||||
|
/* For moderate input not close to the branch point, compute a double
|
||||||
|
approximation as the initial value. */
|
||||||
|
if (fmpz_is_zero(k) &&
|
||||||
|
arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 400) < 0 &&
|
||||||
|
arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 400) < 0 &&
|
||||||
|
(arf_cmp_d(arb_midref(acb_realref(z)), -0.37) < 0 ||
|
||||||
|
arf_cmp_d(arb_midref(acb_realref(z)), -0.36) > 0 ||
|
||||||
|
arf_cmpabs_d(arb_midref(acb_imagref(z)), 0.01) > 0))
|
||||||
|
{
|
||||||
|
acb_lambertw_principal_d(res, z);
|
||||||
|
return 48;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Check if we are close to the branch point at -1/e. */
|
||||||
|
if ((fmpz_is_zero(k) || (fmpz_is_one(k) && arb_is_negative(acb_imagref(z)))
|
||||||
|
|| (fmpz_equal_si(k, -1) && arb_is_nonnegative(acb_imagref(z))))
|
||||||
|
&& ((arf_cmpabs_2exp_si(arb_midref(acb_realref(ez1)), -2) <= 0 &&
|
||||||
|
arf_cmpabs_2exp_si(arb_midref(acb_imagref(ez1)), -2) <= 0)))
|
||||||
|
{
|
||||||
|
acb_t t;
|
||||||
|
acb_init(t);
|
||||||
|
acb_mul_2exp_si(t, ez1, 1);
|
||||||
|
mag_zero(arb_radref(acb_realref(t)));
|
||||||
|
mag_zero(arb_radref(acb_imagref(t)));
|
||||||
|
acb_mul_ui(t, t, 3, prec);
|
||||||
|
acb_sqrt(t, t, prec);
|
||||||
|
if (!fmpz_is_zero(k))
|
||||||
|
acb_neg(t, t);
|
||||||
|
acb_lambertw_branchpoint_series(res, t, 0, prec);
|
||||||
|
acb_clear(t);
|
||||||
|
return 1; /* todo: estimate */
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_lambertw_initial_asymp(res, z, k, prec);
|
||||||
|
return 1; /* todo: estimate */
|
||||||
|
}
|
||||||
|
|
||||||
|
/* note: z should be exact here */
|
||||||
|
void acb_lambertw_main(acb_t res, const acb_t z,
|
||||||
|
const acb_t ez1, const fmpz_t k, int flags, slong prec)
|
||||||
|
{
|
||||||
|
acb_t w, t;
|
||||||
|
mag_t err;
|
||||||
|
slong i, wp, accuracy, ebits, kbits, mbits, wp_initial, extraprec;
|
||||||
|
|
||||||
|
acb_init(t);
|
||||||
|
acb_init(w);
|
||||||
|
mag_init(err);
|
||||||
|
|
||||||
|
/* We need higher precision for large k, large exponents, or very close
|
||||||
|
to the branch point at -1/e. todo: we should be recomputing
|
||||||
|
ez1 to higher precision when close... */
|
||||||
|
acb_get_mag(err, z);
|
||||||
|
if (fmpz_is_zero(k) && mag_cmp_2exp_si(err, 0) < 0)
|
||||||
|
ebits = 0;
|
||||||
|
else
|
||||||
|
ebits = fmpz_bits(MAG_EXPREF(err));
|
||||||
|
|
||||||
|
acb_get_mag(err, ez1);
|
||||||
|
mbits = fmpz_bits(MAG_EXPREF(err));
|
||||||
|
kbits = fmpz_bits(k);
|
||||||
|
|
||||||
|
extraprec = FLINT_MAX(ebits, kbits);
|
||||||
|
extraprec = FLINT_MAX(extraprec, mbits);
|
||||||
|
|
||||||
|
wp = wp_initial = 40 + extraprec;
|
||||||
|
|
||||||
|
accuracy = acb_lambertw_initial(w, z, ez1, k, wp_initial);
|
||||||
|
mag_zero(arb_radref(acb_realref(w)));
|
||||||
|
mag_zero(arb_radref(acb_imagref(w)));
|
||||||
|
|
||||||
|
for (i = 0; i < 5 + FLINT_BIT_COUNT(prec + extraprec); i++)
|
||||||
|
{
|
||||||
|
/* todo: should we restart? */
|
||||||
|
if (!acb_is_finite(w))
|
||||||
|
break;
|
||||||
|
|
||||||
|
wp = FLINT_MIN(3 * accuracy, 1.1 * prec + 10);
|
||||||
|
wp = FLINT_MAX(wp, 40);
|
||||||
|
wp += extraprec;
|
||||||
|
|
||||||
|
acb_lambertw_halley_step(t, z, w, wp);
|
||||||
|
|
||||||
|
/* estimate the error (conservatively) */
|
||||||
|
acb_sub(w, w, t, wp);
|
||||||
|
acb_get_mag(err, w);
|
||||||
|
acb_set(w, t);
|
||||||
|
acb_add_error_mag(t, err);
|
||||||
|
accuracy = acb_rel_accuracy_bits(t);
|
||||||
|
|
||||||
|
if (accuracy > 2 * extraprec)
|
||||||
|
accuracy *= 2.9; /* less conservatively */
|
||||||
|
|
||||||
|
accuracy = FLINT_MIN(accuracy, wp);
|
||||||
|
accuracy = FLINT_MAX(accuracy, 0);
|
||||||
|
|
||||||
|
if (accuracy > prec + extraprec)
|
||||||
|
break;
|
||||||
|
|
||||||
|
mag_zero(arb_radref(acb_realref(w)));
|
||||||
|
mag_zero(arb_radref(acb_imagref(w)));
|
||||||
|
}
|
||||||
|
|
||||||
|
wp = FLINT_MIN(3 * accuracy, 1.1 * prec + 10);
|
||||||
|
wp = FLINT_MAX(wp, 40);
|
||||||
|
wp += extraprec;
|
||||||
|
|
||||||
|
if (acb_lambertw_check_branch(w, k, wp))
|
||||||
|
{
|
||||||
|
acb_t u, r, eu1;
|
||||||
|
mag_t err, rad;
|
||||||
|
|
||||||
|
acb_init(u);
|
||||||
|
acb_init(r);
|
||||||
|
acb_init(eu1);
|
||||||
|
|
||||||
|
mag_init(err);
|
||||||
|
mag_init(rad);
|
||||||
|
|
||||||
|
/* todo: avoid recomputing e^w */
|
||||||
|
acb_exp(t, w, wp);
|
||||||
|
acb_mul(t, t, w, wp);
|
||||||
|
|
||||||
|
acb_sub(r, t, z, wp);
|
||||||
|
|
||||||
|
/* Bound W' on the straight line path between t and z */
|
||||||
|
acb_union(u, t, z, wp);
|
||||||
|
|
||||||
|
arb_const_e(acb_realref(eu1), wp);
|
||||||
|
arb_zero(acb_imagref(eu1));
|
||||||
|
acb_mul(eu1, eu1, u, wp);
|
||||||
|
acb_add_ui(eu1, eu1, 1, wp);
|
||||||
|
|
||||||
|
if (acb_lambertw_branch_crossing(u, eu1, k))
|
||||||
|
{
|
||||||
|
mag_inf(err);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
acb_lambertw_bound_deriv(err, u, eu1, k);
|
||||||
|
acb_get_mag(rad, r);
|
||||||
|
mag_mul(err, err, rad);
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_add_error_mag(w, err);
|
||||||
|
|
||||||
|
acb_set(res, w);
|
||||||
|
|
||||||
|
acb_clear(u);
|
||||||
|
acb_clear(r);
|
||||||
|
acb_clear(eu1);
|
||||||
|
mag_clear(err);
|
||||||
|
mag_clear(rad);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
acb_indeterminate(res);
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_clear(t);
|
||||||
|
acb_clear(w);
|
||||||
|
mag_clear(err);
|
||||||
|
}
|
||||||
|
|
||||||
|
void
|
||||||
|
acb_lambertw_cleared_cut(acb_t res, const acb_t z, const fmpz_t k, int flags, slong prec)
|
||||||
|
{
|
||||||
|
acb_t ez1;
|
||||||
|
acb_init(ez1);
|
||||||
|
|
||||||
|
/* compute z*e + 1 */
|
||||||
|
arb_const_e(acb_realref(ez1), prec);
|
||||||
|
acb_mul(ez1, ez1, z, prec);
|
||||||
|
acb_add_ui(ez1, ez1, 1, prec);
|
||||||
|
|
||||||
|
if (acb_is_exact(z))
|
||||||
|
{
|
||||||
|
acb_lambertw_main(res, z, ez1, k, flags, prec);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
acb_t zz;
|
||||||
|
mag_t err, rad;
|
||||||
|
|
||||||
|
mag_init(err);
|
||||||
|
mag_init(rad);
|
||||||
|
acb_init(zz);
|
||||||
|
|
||||||
|
acb_lambertw_bound_deriv(err, z, ez1, k);
|
||||||
|
mag_hypot(rad, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
|
||||||
|
mag_mul(err, err, rad);
|
||||||
|
|
||||||
|
acb_set(zz, z);
|
||||||
|
mag_zero(arb_radref(acb_realref(zz)));
|
||||||
|
mag_zero(arb_radref(acb_imagref(zz))); /* todo: recompute ez1? */
|
||||||
|
|
||||||
|
acb_lambertw_main(res, zz, ez1, k, flags, prec);
|
||||||
|
acb_add_error_mag(res, err);
|
||||||
|
|
||||||
|
mag_clear(err);
|
||||||
|
mag_clear(rad);
|
||||||
|
acb_clear(zz);
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_clear(ez1);
|
||||||
|
}
|
||||||
|
|
||||||
|
void
|
||||||
|
_acb_lambertw(acb_t res, const acb_t z, const acb_t ez1, const fmpz_t k, int flags, slong prec)
|
||||||
|
{
|
||||||
|
slong goal, ebits, ebits2, ls, lt;
|
||||||
|
const fmpz * expo;
|
||||||
|
|
||||||
|
/* Estimated accuracy goal. */
|
||||||
|
/* todo: account for exponent bits and bits in k. */
|
||||||
|
goal = acb_rel_accuracy_bits(z);
|
||||||
|
goal = FLINT_MAX(goal, 10);
|
||||||
|
goal = FLINT_MIN(goal, prec);
|
||||||
|
|
||||||
|
/* Handle tiny z directly. For k >= 2, |c_k| <= 4^k / 16. */
|
||||||
|
if (fmpz_is_zero(k)
|
||||||
|
&& arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), -goal / 2) < 0
|
||||||
|
&& arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), -goal / 2) < 0)
|
||||||
|
{
|
||||||
|
mag_t err;
|
||||||
|
mag_init(err);
|
||||||
|
acb_get_mag(err, z);
|
||||||
|
mag_mul_2exp_si(err, err, 2);
|
||||||
|
acb_set(res, z);
|
||||||
|
acb_submul(res, res, res, prec);
|
||||||
|
mag_geom_series(err, err, 3);
|
||||||
|
mag_mul_2exp_si(err, err, -4);
|
||||||
|
acb_add_error_mag(res, err);
|
||||||
|
mag_clear(err);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (arf_cmpabs(arb_midref(acb_realref(z)), arb_midref(acb_imagref(z))) >= 0)
|
||||||
|
expo = ARF_EXPREF(arb_midref(acb_realref(z)));
|
||||||
|
else
|
||||||
|
expo = ARF_EXPREF(arb_midref(acb_imagref(z)));
|
||||||
|
|
||||||
|
ebits = fmpz_bits(expo);
|
||||||
|
|
||||||
|
/* ebits ~= log2(|log(z) + 2 pi i k|) */
|
||||||
|
/* ebits2 ~= log2(log(log(z))) */
|
||||||
|
ebits = FLINT_MAX(ebits, fmpz_bits(k));
|
||||||
|
ebits = FLINT_MAX(ebits, 1) - 1;
|
||||||
|
ebits2 = FLINT_BIT_COUNT(ebits);
|
||||||
|
ebits2 = FLINT_MAX(ebits2, 1) - 1;
|
||||||
|
|
||||||
|
/* We gain accuracy from the exponent when W ~ log - log log */
|
||||||
|
if (fmpz_sgn(expo) > 0 || (fmpz_sgn(expo) < 0 && !fmpz_is_zero(k)))
|
||||||
|
{
|
||||||
|
goal += ebits - ebits2;
|
||||||
|
goal = FLINT_MAX(goal, 10);
|
||||||
|
goal = FLINT_MIN(goal, prec);
|
||||||
|
|
||||||
|
/* The asymptotic series with truncation L, M gives us about
|
||||||
|
t - max(2+lt+L*(2+ls), M*(2+lt)) bits of accuracy where
|
||||||
|
ls = -ebits, lt = ebits2 - ebits. */
|
||||||
|
ls = 2 - ebits;
|
||||||
|
lt = 2 + ebits2 - ebits;
|
||||||
|
|
||||||
|
if (ebits - FLINT_MAX(lt + 1*ls, 1*lt) > goal)
|
||||||
|
{
|
||||||
|
acb_lambertw_asymp(res, z, k, 1, 1, goal);
|
||||||
|
acb_set_round(res, res, prec);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
else if (ebits - FLINT_MAX(lt + 3*ls, 5*lt) > goal)
|
||||||
|
{
|
||||||
|
acb_lambertw_asymp(res, z, k, 3, 5, goal);
|
||||||
|
acb_set_round(res, res, prec);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Extremely close to the branch point at -1/e, use the series expansion directly. */
|
||||||
|
if (fmpz_is_zero(k) || (fmpz_is_one(k) && arb_is_negative(acb_imagref(z)))
|
||||||
|
|| (fmpz_equal_si(k, -1) && arb_is_nonnegative(acb_imagref(z))))
|
||||||
|
{
|
||||||
|
if (acb_contains_zero(ez1) ||
|
||||||
|
(arf_cmpabs_2exp_si(arb_midref(acb_realref(ez1)), -goal / 4.5) < 0 &&
|
||||||
|
arf_cmpabs_2exp_si(arb_midref(acb_imagref(ez1)), -goal / 4.5) < 0))
|
||||||
|
{
|
||||||
|
acb_t t;
|
||||||
|
acb_init(t);
|
||||||
|
acb_mul_2exp_si(t, ez1, 1);
|
||||||
|
acb_sqrt(t, t, goal);
|
||||||
|
if (!fmpz_is_zero(k))
|
||||||
|
acb_neg(t, t);
|
||||||
|
acb_lambertw_branchpoint_series(res, t, 1, goal);
|
||||||
|
acb_set_round(res, res, prec);
|
||||||
|
acb_clear(t);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* todo: compute union of two results */
|
||||||
|
if (acb_lambertw_branch_crossing(z, ez1, k))
|
||||||
|
{
|
||||||
|
acb_indeterminate(res);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
acb_t zz, zmid, zmide1;
|
||||||
|
arf_t eps;
|
||||||
|
|
||||||
|
acb_init(zz);
|
||||||
|
acb_init(zmid);
|
||||||
|
acb_init(zmide1);
|
||||||
|
arf_init(eps);
|
||||||
|
|
||||||
|
arf_mul_2exp_si(eps, arb_midref(acb_realref(z)), -goal);
|
||||||
|
acb_set(zz, z);
|
||||||
|
|
||||||
|
if (arf_sgn(arb_midref(acb_realref(zz))) < 0 &&
|
||||||
|
(!fmpz_is_zero(k) || arf_sgn(arb_midref(acb_realref(ez1))) < 0) &&
|
||||||
|
arf_cmpabs(arb_midref(acb_imagref(zz)), eps) < 0)
|
||||||
|
{
|
||||||
|
/* now the value must be in [0,2eps] */
|
||||||
|
arf_get_mag(arb_radref(acb_imagref(zz)), eps);
|
||||||
|
arf_set_mag(arb_midref(acb_imagref(zz)), arb_radref(acb_imagref(zz)));
|
||||||
|
|
||||||
|
if (arf_sgn(arb_midref(acb_imagref(z))) >= 0)
|
||||||
|
{
|
||||||
|
acb_lambertw_cleared_cut(res, zz, k, flags, goal);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
fmpz_t kk;
|
||||||
|
fmpz_init(kk);
|
||||||
|
fmpz_neg(kk, k);
|
||||||
|
acb_lambertw_cleared_cut(res, zz, kk, flags, goal);
|
||||||
|
acb_conj(res, res);
|
||||||
|
fmpz_clear(kk);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
acb_lambertw_cleared_cut(res, zz, k, flags, goal);
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_set_round(res, res, prec);
|
||||||
|
|
||||||
|
acb_clear(zz);
|
||||||
|
acb_clear(zmid);
|
||||||
|
acb_clear(zmide1);
|
||||||
|
arf_clear(eps);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void
|
||||||
|
acb_lambertw(acb_t res, const acb_t z, const fmpz_t k, int flags, slong prec)
|
||||||
|
{
|
||||||
|
acb_t ez1;
|
||||||
|
|
||||||
|
if (!acb_is_finite(z) || (!fmpz_is_zero(k) && acb_contains_zero(z)))
|
||||||
|
{
|
||||||
|
acb_indeterminate(res);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_init(ez1);
|
||||||
|
|
||||||
|
/* precompute z*e + 1 */
|
||||||
|
arb_const_e(acb_realref(ez1), prec);
|
||||||
|
acb_mul(ez1, ez1, z, prec);
|
||||||
|
acb_add_ui(ez1, ez1, 1, prec);
|
||||||
|
|
||||||
|
/* use real code when possible */
|
||||||
|
if (acb_is_real(z) && arb_is_positive(acb_realref(ez1)) &&
|
||||||
|
(fmpz_is_zero(k) ||
|
||||||
|
(fmpz_equal_si(k, -1) && arb_is_negative(acb_realref(z)))))
|
||||||
|
{
|
||||||
|
arb_lambertw(acb_realref(res), acb_realref(z), !fmpz_is_zero(k), prec);
|
||||||
|
arb_zero(acb_imagref(res));
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
_acb_lambertw(res, z, ez1, k, flags, prec);
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_clear(ez1);
|
||||||
|
}
|
||||||
|
|
||||||
144
acb/lambertw_asymp.c
Normal file
144
acb/lambertw_asymp.c
Normal file
|
|
@ -0,0 +1,144 @@
|
||||||
|
/*
|
||||||
|
Copyright (C) 2017 Fredrik Johansson
|
||||||
|
|
||||||
|
This file is part of Arb.
|
||||||
|
|
||||||
|
Arb is free software: you can redistribute it and/or modify it under
|
||||||
|
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||||
|
by the Free Software Foundation; either version 2.1 of the License, or
|
||||||
|
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#include "acb.h"
|
||||||
|
|
||||||
|
void
|
||||||
|
acb_lambertw_asymp(acb_t res, const acb_t z, const fmpz_t k, slong L, slong M, slong prec)
|
||||||
|
{
|
||||||
|
acb_t L1, L2, sigma, tau, s, c, u;
|
||||||
|
slong l, m;
|
||||||
|
fmpz_t t;
|
||||||
|
fmpz * sc;
|
||||||
|
|
||||||
|
/* For k = 0, the asymptotic expansion is not valid near 0. */
|
||||||
|
/* (It is sufficient to look at the midpoint as a test here.) */
|
||||||
|
if (fmpz_is_zero(k) && arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 0) < 0
|
||||||
|
&& arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 0) < 0)
|
||||||
|
{
|
||||||
|
acb_indeterminate(res);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_init(L1);
|
||||||
|
acb_init(L2);
|
||||||
|
acb_init(sigma);
|
||||||
|
acb_init(tau);
|
||||||
|
acb_init(s);
|
||||||
|
acb_init(c);
|
||||||
|
acb_init(u);
|
||||||
|
fmpz_init(t);
|
||||||
|
|
||||||
|
acb_const_pi(L2, prec);
|
||||||
|
acb_mul_2exp_si(L2, L2, 1);
|
||||||
|
acb_mul_fmpz(L2, L2, k, prec);
|
||||||
|
acb_mul_onei(L2, L2);
|
||||||
|
acb_log(L1, z, prec);
|
||||||
|
acb_add(L1, L1, L2, prec);
|
||||||
|
acb_log(L2, L1, prec);
|
||||||
|
|
||||||
|
acb_inv(sigma, L1, prec);
|
||||||
|
acb_mul(tau, L2, sigma, prec);
|
||||||
|
|
||||||
|
acb_zero(s);
|
||||||
|
|
||||||
|
/* Stirling numbers */
|
||||||
|
sc = _fmpz_vec_init(L);
|
||||||
|
|
||||||
|
acb_one(u);
|
||||||
|
|
||||||
|
for (m = 1; m < M; m++)
|
||||||
|
{
|
||||||
|
if (m == 1)
|
||||||
|
{
|
||||||
|
for (l = 0; l < L; l++)
|
||||||
|
fmpz_one(sc + l);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
for (l = 0; l < L; l++)
|
||||||
|
{
|
||||||
|
fmpz_mul_ui(sc + l, sc + l, m + l - 1);
|
||||||
|
if (l > 0)
|
||||||
|
fmpz_add(sc + l, sc + l, sc + l - 1);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_zero(c);
|
||||||
|
|
||||||
|
/* todo: precompute powers instead of horner */
|
||||||
|
for (l = L - 1; l >= 0; l--)
|
||||||
|
{
|
||||||
|
acb_mul(c, c, sigma, prec);
|
||||||
|
if (l % 2)
|
||||||
|
acb_sub_fmpz(c, c, sc + l, prec);
|
||||||
|
else
|
||||||
|
acb_add_fmpz(c, c, sc + l, prec);
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_mul(u, u, tau, prec);
|
||||||
|
acb_div_ui(u, u, m, prec);
|
||||||
|
acb_addmul(s, c, u, prec);
|
||||||
|
}
|
||||||
|
|
||||||
|
_fmpz_vec_clear(sc, L);
|
||||||
|
|
||||||
|
acb_sub(s, s, L2, prec);
|
||||||
|
acb_add(s, s, L1, prec);
|
||||||
|
|
||||||
|
{
|
||||||
|
mag_t m4s, m4t, one, q, r;
|
||||||
|
|
||||||
|
mag_init(m4s);
|
||||||
|
mag_init(m4t);
|
||||||
|
mag_init(one);
|
||||||
|
mag_init(q);
|
||||||
|
mag_init(r);
|
||||||
|
|
||||||
|
acb_get_mag(m4s, sigma);
|
||||||
|
mag_mul_2exp_si(m4s, m4s, 2);
|
||||||
|
acb_get_mag(m4t, tau);
|
||||||
|
mag_mul_2exp_si(m4t, m4t, 2);
|
||||||
|
|
||||||
|
mag_one(one);
|
||||||
|
|
||||||
|
mag_sub_lower(q, one, m4s);
|
||||||
|
mag_sub_lower(r, one, m4t);
|
||||||
|
mag_mul(q, q, r);
|
||||||
|
|
||||||
|
mag_pow_ui(r, m4s, L);
|
||||||
|
mag_mul(r, r, m4t);
|
||||||
|
mag_pow_ui(m4t, m4t, M);
|
||||||
|
mag_add(r, r, m4t);
|
||||||
|
|
||||||
|
mag_div(q, r, q);
|
||||||
|
|
||||||
|
acb_add_error_mag(s, q);
|
||||||
|
|
||||||
|
mag_clear(m4s);
|
||||||
|
mag_clear(m4t);
|
||||||
|
mag_clear(one);
|
||||||
|
mag_clear(q);
|
||||||
|
mag_clear(r);
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_set(res, s);
|
||||||
|
|
||||||
|
acb_clear(sigma);
|
||||||
|
acb_clear(tau);
|
||||||
|
acb_clear(s);
|
||||||
|
acb_clear(c);
|
||||||
|
acb_clear(L1);
|
||||||
|
acb_clear(L2);
|
||||||
|
acb_clear(u);
|
||||||
|
fmpz_clear(t);
|
||||||
|
}
|
||||||
|
|
||||||
90
acb/lambertw_bound_deriv.c
Normal file
90
acb/lambertw_bound_deriv.c
Normal file
|
|
@ -0,0 +1,90 @@
|
||||||
|
/*
|
||||||
|
Copyright (C) 2017 Fredrik Johansson
|
||||||
|
|
||||||
|
This file is part of Arb.
|
||||||
|
|
||||||
|
Arb is free software: you can redistribute it and/or modify it under
|
||||||
|
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||||
|
by the Free Software Foundation; either version 2.1 of the License, or
|
||||||
|
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#include "acb.h"
|
||||||
|
|
||||||
|
void
|
||||||
|
acb_lambertw_bound_deriv(mag_t res, const acb_t z, const acb_t ez1, const fmpz_t k)
|
||||||
|
{
|
||||||
|
mag_t t, u;
|
||||||
|
|
||||||
|
mag_init(t);
|
||||||
|
mag_init(u);
|
||||||
|
|
||||||
|
/* Main approximation: W'(z) ~ 1/|z| */
|
||||||
|
if (fmpz_is_zero(k))
|
||||||
|
{
|
||||||
|
/* 1/(1+|z|) */
|
||||||
|
acb_get_mag_lower(res, z);
|
||||||
|
mag_one(t);
|
||||||
|
mag_add_lower(res, res, t);
|
||||||
|
mag_div(res, t, res);
|
||||||
|
}
|
||||||
|
else if (fmpz_is_pm1(k))
|
||||||
|
{
|
||||||
|
/* 2/|z| */
|
||||||
|
mag_set_ui(t, 2);
|
||||||
|
acb_get_mag_lower(res, z);
|
||||||
|
mag_div(res, t, res);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
/* |W'(z)| = |1/z| |W(z)/(1+W(z))|
|
||||||
|
<= |1/z| (2pi) / (2pi-1) */
|
||||||
|
mag_set_ui_2exp_si(t, 77, -6);
|
||||||
|
acb_get_mag_lower(res, z);
|
||||||
|
mag_div(res, t, res);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Compute correction near the branch point */
|
||||||
|
if (fmpz_is_zero(k) || fmpz_equal_si(k,-1) ||
|
||||||
|
(fmpz_is_one(k) && !arb_is_nonnegative(acb_imagref(z))))
|
||||||
|
{
|
||||||
|
acb_t b;
|
||||||
|
acb_init(b);
|
||||||
|
|
||||||
|
if (fmpz_is_zero(k))
|
||||||
|
{
|
||||||
|
/* [-4,1] + [-2,2]i */
|
||||||
|
arf_set_si_2exp_si(arb_midref(acb_realref(b)), -3, -1);
|
||||||
|
mag_set_ui_2exp_si(arb_radref(acb_realref(b)), 5, -1);
|
||||||
|
arf_zero(arb_midref(acb_imagref(b)));
|
||||||
|
mag_set_ui(arb_radref(acb_imagref(b)), 2);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
/* k = 1 [-1/2,-1/4] + [-1/8,0) i */
|
||||||
|
/* k = -1 [-1/2,-1/4] + [0,1/8] i */
|
||||||
|
arf_set_si_2exp_si(arb_midref(acb_realref(b)), -3, -3);
|
||||||
|
mag_set_ui_2exp_si(arb_radref(acb_realref(b)), 1, -3);
|
||||||
|
if (fmpz_is_one(k))
|
||||||
|
arf_set_si_2exp_si(arb_midref(acb_imagref(b)), -1, -4);
|
||||||
|
else
|
||||||
|
arf_set_ui_2exp_si(arb_midref(acb_imagref(b)), 1, -4);
|
||||||
|
mag_set_ui_2exp_si(arb_radref(acb_imagref(b)), 1, -4);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (acb_overlaps(z, b))
|
||||||
|
{
|
||||||
|
/* add 2/sqrt(|ez+1|) */
|
||||||
|
acb_get_mag_lower(u, ez1);
|
||||||
|
mag_rsqrt(u, u);
|
||||||
|
mag_mul_2exp_si(u, u, 1);
|
||||||
|
mag_add(res, res, u);
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_clear(b);
|
||||||
|
}
|
||||||
|
|
||||||
|
mag_clear(t);
|
||||||
|
mag_clear(u);
|
||||||
|
}
|
||||||
|
|
||||||
88
acb/lambertw_check_branch.c
Normal file
88
acb/lambertw_check_branch.c
Normal file
|
|
@ -0,0 +1,88 @@
|
||||||
|
/*
|
||||||
|
Copyright (C) 2017 Fredrik Johansson
|
||||||
|
|
||||||
|
This file is part of Arb.
|
||||||
|
|
||||||
|
Arb is free software: you can redistribute it and/or modify it under
|
||||||
|
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||||
|
by the Free Software Foundation; either version 2.1 of the License, or
|
||||||
|
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#include "acb.h"
|
||||||
|
|
||||||
|
/* todo: use reduced precision test first */
|
||||||
|
int
|
||||||
|
acb_lambertw_check_branch(const acb_t w, const fmpz_t k, slong prec)
|
||||||
|
{
|
||||||
|
arb_t t, u, v, a, b;
|
||||||
|
int res = 0;
|
||||||
|
|
||||||
|
arb_init(t);
|
||||||
|
arb_init(u);
|
||||||
|
arb_init(v);
|
||||||
|
arb_init(a);
|
||||||
|
arb_init(b);
|
||||||
|
|
||||||
|
/* t = x sinc(y), v = -cos(y) */
|
||||||
|
arb_sinc(t, acb_imagref(w), prec);
|
||||||
|
arb_mul(t, t, acb_realref(w), prec);
|
||||||
|
arb_cos(v, acb_imagref(w), prec);
|
||||||
|
arb_neg(v, v);
|
||||||
|
|
||||||
|
/* u = y / pi, with conjugate relation for k */
|
||||||
|
arb_const_pi(u, prec);
|
||||||
|
arb_div(u, acb_imagref(w), u, prec);
|
||||||
|
if (fmpz_sgn(k) < 0)
|
||||||
|
arb_neg(u, u);
|
||||||
|
|
||||||
|
if (fmpz_is_zero(k))
|
||||||
|
{
|
||||||
|
/* -1 < u < 1 and t > v */
|
||||||
|
arb_set_si(a, -1);
|
||||||
|
arb_set_si(b, 1);
|
||||||
|
|
||||||
|
if (arb_gt(u, a) && arb_lt(u, b) && arb_gt(t, v))
|
||||||
|
{
|
||||||
|
res = 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
arb_set_fmpz(a, k);
|
||||||
|
arb_abs(a, a);
|
||||||
|
arb_mul_2exp_si(a, a, 1);
|
||||||
|
arb_add_ui(b, a, 1, prec);
|
||||||
|
arb_sub_ui(a, a, 2, prec);
|
||||||
|
|
||||||
|
/* if u > 2|k|-2 and u < 2|k|+1 */
|
||||||
|
if (arb_gt(u, a) && arb_lt(u, b))
|
||||||
|
{
|
||||||
|
arb_add_ui(a, a, 1, prec);
|
||||||
|
arb_sub_ui(b, b, 1, prec);
|
||||||
|
|
||||||
|
/* if u > 2|k|-1 and u < 2|k| */
|
||||||
|
if (arb_gt(u, a) && arb_lt(u, b))
|
||||||
|
{
|
||||||
|
res = 1;
|
||||||
|
}
|
||||||
|
else if (arb_lt(u, b) && arb_lt(t, v)) /* u < 2|k| and t < v */
|
||||||
|
{
|
||||||
|
res = 1;
|
||||||
|
}
|
||||||
|
else if (arb_gt(u, a) && arb_gt(t, v)) /* u > 2|k|-1 and t > v */
|
||||||
|
{
|
||||||
|
res = 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
arb_clear(t);
|
||||||
|
arb_clear(u);
|
||||||
|
arb_clear(v);
|
||||||
|
arb_clear(a);
|
||||||
|
arb_clear(b);
|
||||||
|
|
||||||
|
return res;
|
||||||
|
}
|
||||||
|
|
||||||
152
acb/test/t-lambertw.c
Normal file
152
acb/test/t-lambertw.c
Normal file
|
|
@ -0,0 +1,152 @@
|
||||||
|
/*
|
||||||
|
Copyright (C) 2017 Fredrik Johansson
|
||||||
|
|
||||||
|
This file is part of Arb.
|
||||||
|
|
||||||
|
Arb is free software: you can redistribute it and/or modify it under
|
||||||
|
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||||
|
by the Free Software Foundation; either version 2.1 of the License, or
|
||||||
|
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#include "acb.h"
|
||||||
|
|
||||||
|
int main()
|
||||||
|
{
|
||||||
|
slong iter;
|
||||||
|
flint_rand_t state;
|
||||||
|
|
||||||
|
flint_printf("lambertw....");
|
||||||
|
fflush(stdout);
|
||||||
|
|
||||||
|
flint_randinit(state);
|
||||||
|
|
||||||
|
for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
|
||||||
|
{
|
||||||
|
acb_t x1, x2, t, w1, w2;
|
||||||
|
slong prec1, prec2, ebits;
|
||||||
|
fmpz_t k;
|
||||||
|
|
||||||
|
acb_init(x1);
|
||||||
|
acb_init(x2);
|
||||||
|
acb_init(t);
|
||||||
|
acb_init(w1);
|
||||||
|
acb_init(w2);
|
||||||
|
fmpz_init(k);
|
||||||
|
|
||||||
|
if (n_randint(state, 4) == 0)
|
||||||
|
{
|
||||||
|
prec1 = 2 + n_randint(state, 3000);
|
||||||
|
prec2 = 2 + n_randint(state, 3000);
|
||||||
|
ebits = 1 + n_randint(state, 1000);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
prec1 = 2 + n_randint(state, 300);
|
||||||
|
prec2 = 2 + n_randint(state, 300);
|
||||||
|
ebits = 1 + n_randint(state, 50);
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_randtest(x1, state, 1 + n_randint(state, 1000), ebits);
|
||||||
|
acb_randtest(x2, state, 1 + n_randint(state, 1000), ebits);
|
||||||
|
acb_randtest(t, state, 1 + n_randint(state, 1000), ebits);
|
||||||
|
acb_randtest(w1, state, 1 + n_randint(state, 1000), ebits);
|
||||||
|
acb_randtest(w2, state, 1 + n_randint(state, 1000), ebits);
|
||||||
|
fmpz_randtest(k, state, ebits);
|
||||||
|
|
||||||
|
if (n_randint(state, 4) == 0)
|
||||||
|
{
|
||||||
|
arb_const_e(acb_realref(t), 2 * prec1);
|
||||||
|
arb_inv(acb_realref(t), acb_realref(t), 2 * prec1);
|
||||||
|
arb_sub(acb_realref(x1), acb_realref(x1), acb_realref(t), 2 * prec1);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (n_randint(state, 2))
|
||||||
|
{
|
||||||
|
acb_set(x2, x1);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
acb_add(x2, x1, t, 2 * prec1);
|
||||||
|
acb_sub(x2, x2, t, 2 * prec1);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (n_randint(state, 4) == 0)
|
||||||
|
acb_lambertw_asymp(w1, x1, k,
|
||||||
|
1 + n_randint(state, 10), 1 + n_randint(state, 10), prec1);
|
||||||
|
else
|
||||||
|
acb_lambertw(w1, x1, k, 0, prec1);
|
||||||
|
|
||||||
|
if (n_randint(state, 2))
|
||||||
|
{
|
||||||
|
acb_set(w2, x2);
|
||||||
|
acb_lambertw(w2, w2, k, 0, prec1);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
acb_lambertw(w2, x2, k, 0, prec1);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (!acb_overlaps(w1, w2))
|
||||||
|
{
|
||||||
|
flint_printf("FAIL: overlap\n\n");
|
||||||
|
flint_printf("iter %wd, branch = ", iter); fmpz_print(k);
|
||||||
|
flint_printf("prec1 = %wd, prec2 = %wd\n\n", prec1, prec2);
|
||||||
|
flint_printf("x1 = "); acb_printd(x1, 50); flint_printf("\n\n");
|
||||||
|
flint_printf("x2 = "); acb_printd(x2, 50); flint_printf("\n\n");
|
||||||
|
flint_printf("w1 = "); acb_printd(w1, 50); flint_printf("\n\n");
|
||||||
|
flint_printf("w2 = "); acb_printd(w2, 50); flint_printf("\n\n");
|
||||||
|
flint_abort();
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_exp(t, w1, prec1);
|
||||||
|
acb_mul(t, t, w1, prec1);
|
||||||
|
|
||||||
|
if (!acb_contains(t, x1))
|
||||||
|
{
|
||||||
|
flint_printf("FAIL: functional equation\n\n");
|
||||||
|
flint_printf("iter %wd, branch = ", iter); fmpz_print(k);
|
||||||
|
flint_printf("prec1 = %wd, prec2 = %wd\n\n", prec1, prec2);
|
||||||
|
flint_printf("x1 = "); acb_printd(x1, 50); flint_printf("\n\n");
|
||||||
|
flint_printf("x2 = "); acb_printd(x2, 50); flint_printf("\n\n");
|
||||||
|
flint_printf("w1 = "); acb_printd(w1, 50); flint_printf("\n\n");
|
||||||
|
flint_printf("w2 = "); acb_printd(w2, 50); flint_printf("\n\n");
|
||||||
|
flint_printf("t = "); acb_printd(t, 50); flint_printf("\n\n");
|
||||||
|
flint_abort();
|
||||||
|
}
|
||||||
|
|
||||||
|
if (!arb_contains_zero(acb_imagref(x1)))
|
||||||
|
{
|
||||||
|
acb_conj(x2, x1);
|
||||||
|
fmpz_neg(k, k);
|
||||||
|
acb_lambertw(w2, x2, k, 0, prec2);
|
||||||
|
fmpz_neg(k, k);
|
||||||
|
acb_conj(w2, w2);
|
||||||
|
|
||||||
|
if (!acb_overlaps(w1, w2))
|
||||||
|
{
|
||||||
|
flint_printf("FAIL: conjugation\n\n");
|
||||||
|
flint_printf("iter %wd, branch = ", iter); fmpz_print(k);
|
||||||
|
flint_printf("prec1 = %wd, prec2 = %wd\n\n", prec1, prec2);
|
||||||
|
flint_printf("x1 = "); acb_printd(x1, 50); flint_printf("\n\n");
|
||||||
|
flint_printf("x2 = "); acb_printd(x2, 50); flint_printf("\n\n");
|
||||||
|
flint_printf("w1 = "); acb_printd(w1, 50); flint_printf("\n\n");
|
||||||
|
flint_printf("w2 = "); acb_printd(w2, 50); flint_printf("\n\n");
|
||||||
|
flint_abort();
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
acb_clear(x1);
|
||||||
|
acb_clear(x2);
|
||||||
|
acb_clear(t);
|
||||||
|
acb_clear(w1);
|
||||||
|
acb_clear(w2);
|
||||||
|
fmpz_clear(k);
|
||||||
|
}
|
||||||
|
|
||||||
|
flint_randclear(state);
|
||||||
|
flint_cleanup();
|
||||||
|
flint_printf("PASS\n");
|
||||||
|
return EXIT_SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
@ -260,6 +260,10 @@ Precision and comparisons
|
||||||
|
|
||||||
Returns nonzero iff *x* and *y* have some point in common.
|
Returns nonzero iff *x* and *y* have some point in common.
|
||||||
|
|
||||||
|
.. function:: void acb_union(acb_t z, const acb_t x, const acb_t y, slong prec)
|
||||||
|
|
||||||
|
Sets *z* to a complex interval containing both *x* and *y*.
|
||||||
|
|
||||||
.. function:: void acb_get_abs_ubound_arf(arf_t u, const acb_t z, slong prec)
|
.. function:: void acb_get_abs_ubound_arf(arf_t u, const acb_t z, slong prec)
|
||||||
|
|
||||||
Sets *u* to an upper bound for the absolute value of *z*, computed
|
Sets *u* to an upper bound for the absolute value of *z*, computed
|
||||||
|
|
@ -672,6 +676,44 @@ Inverse hyperbolic functions
|
||||||
|
|
||||||
Sets *res* to `\operatorname{atanh}(z) = -i \operatorname{atan}(iz)`.
|
Sets *res* to `\operatorname{atanh}(z) = -i \operatorname{atan}(iz)`.
|
||||||
|
|
||||||
|
Lambert W function
|
||||||
|
-------------------------------------------------------------------------------
|
||||||
|
|
||||||
|
.. function:: void acb_lambertw_asymp(acb_t res, const acb_t z, const fmpz_t k, slong L, slong M, slong prec)
|
||||||
|
|
||||||
|
Sets *res* to the Lambert W function `W_k(z)` computed using *L* and *M*
|
||||||
|
terms in the bivariate series giving the asymptotic expansion at
|
||||||
|
zero or infinity. This algorithm is valid
|
||||||
|
everywhere, but the error bound is only finite when `|\log(z)|` is
|
||||||
|
sufficiently large.
|
||||||
|
|
||||||
|
.. function:: int acb_lambertw_check_branch(const acb_t w, const fmpz_t k, slong prec)
|
||||||
|
|
||||||
|
Tests if *w* definitely lies in the image of the branch `W_k(z)`.
|
||||||
|
This function is used internally to verify that a computed approximation
|
||||||
|
of the Lambert W function lies on the intended branch. Note that this will
|
||||||
|
necessarily evaluate to false for points exactly on (or overlapping) the
|
||||||
|
branch cuts, where a different algorithm has to be used.
|
||||||
|
|
||||||
|
.. function:: void acb_lambertw_bound_deriv(mag_t res, const acb_t z, const acb_t ez1, const fmpz_t k)
|
||||||
|
|
||||||
|
Sets *res* to an upper bound for `|W_k'(z)|`. The input *ez1* should
|
||||||
|
contain the precomputed value of `ez+1`.
|
||||||
|
|
||||||
|
Along the real line, the directional derivative of `W_k(z)` is understood
|
||||||
|
to be taken. As a result, the user must handle the branch cut
|
||||||
|
discontinuity separately when using this function to bound perturbations
|
||||||
|
in the value of `W_k(z)`.
|
||||||
|
|
||||||
|
.. function:: void acb_lambertw(acb_t res, const acb_t z, const fmpz_t k, int flags, slong prec)
|
||||||
|
|
||||||
|
Sets *res* to the Lambert W function `W_k(z)` where the index *k* selects
|
||||||
|
an arbitrary branch (with `k = 0` giving the principal branch).
|
||||||
|
The placement of branch cuts follows [CGHJK1996]_.
|
||||||
|
|
||||||
|
The *flags* argument is currently unused. In a future version, it might
|
||||||
|
allow further control over branch cuts.
|
||||||
|
|
||||||
Rising factorials
|
Rising factorials
|
||||||
-------------------------------------------------------------------------------
|
-------------------------------------------------------------------------------
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Add table
Reference in a new issue