From 8947049b77a4d8f01c45025b7bd45d90972d0881 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Tue, 10 Oct 2017 22:34:09 +0200 Subject: [PATCH] fix and refactor derivative bounds for legendre polynomials --- arb_hypgeom.h | 1 + arb_hypgeom/legendre_p_ui.c | 77 ++++++++++----- .../test/t-legendre_p_ui_deriv_bound.c | 98 +++++++++++++++++++ doc/source/arb_hypgeom.rst | 8 ++ 4 files changed, 158 insertions(+), 26 deletions(-) create mode 100644 arb_hypgeom/test/t-legendre_p_ui_deriv_bound.c diff --git a/arb_hypgeom.h b/arb_hypgeom.h index 73597c89..5500994f 100644 --- a/arb_hypgeom.h +++ b/arb_hypgeom.h @@ -104,6 +104,7 @@ void arb_hypgeom_hermite_h(arb_t res, const arb_t nu, const arb_t z, slong prec) void arb_hypgeom_legendre_p(arb_t res, const arb_t n, const arb_t m, const arb_t z, int type, slong prec); void arb_hypgeom_legendre_q(arb_t res, const arb_t n, const arb_t m, const arb_t z, int type, slong prec); +void arb_hypgeom_legendre_p_ui_deriv_bound(mag_t dp, mag_t dp2, ulong n, const arb_t x, const arb_t x2sub1); void arb_hypgeom_legendre_p_ui_asymp(arb_t res, arb_t res2, ulong n, const arb_t x, slong K, slong prec); void arb_hypgeom_legendre_p_ui_one(arb_t res, arb_t res2, ulong n, const arb_t x, slong K, slong prec); void arb_hypgeom_legendre_p_ui_zero(arb_t res, arb_t res2, ulong n, const arb_t x, slong K, slong prec); diff --git a/arb_hypgeom/legendre_p_ui.c b/arb_hypgeom/legendre_p_ui.c index f0721b12..dc39a4c8 100644 --- a/arb_hypgeom/legendre_p_ui.c +++ b/arb_hypgeom/legendre_p_ui.c @@ -26,6 +26,54 @@ static double log2_bin_uiui_fast(ulong n, ulong k) return n * htab[FLINT_MIN(k, 15)]; } +void +arb_hypgeom_legendre_p_ui_deriv_bound(mag_t dp, mag_t dp2, ulong n, const arb_t x, const arb_t x2sub1) +{ + mag_t t, u, xm, x2sub1m; + + mag_init(t); + mag_init(u); + mag_init(xm); + mag_init(x2sub1m); + + arb_get_mag(xm, x); + arb_get_mag_lower(x2sub1m, x2sub1); + + /* |P'(x)| <= min(n(n+1)/2, n/sqrt(1-x^2)) */ + mag_set_ui(t, n); + mag_add_ui(t, t, 1); + mag_mul_2exp_si(t, t, -1); + mag_rsqrt(u, x2sub1m); + mag_max(t, t, u); + mag_mul_ui(dp, t, n); + + /* |P''(x)| <= min((n+2)(n+1)n(n-1)/8, (2x|P'(x)| + n(n+1))/(1-x^2)) */ + mag_mul(dp2, dp, xm); + mag_mul_2exp_si(dp2, dp2, 1); + + mag_set_ui(t, n); + mag_add_ui(t, t, 1); + mag_mul_ui(t, t, n); + mag_add(dp2, dp2, t); + mag_div(dp2, dp2, x2sub1m); + + mag_set_ui(t, n); + mag_add_ui(t, t, 2); + mag_set_ui(u, n); + mag_add_ui(u, u, 1); + mag_mul(t, t, u); + mag_mul_ui(t, t, n); + mag_mul_ui(t, t, n - 1); + mag_mul_2exp_si(t, t, -3); + mag_min(dp2, dp2, t); + + + mag_clear(t); + mag_clear(u); + mag_clear(xm); + mag_clear(x2sub1m); +} + void arb_hypgeom_legendre_p_ui(arb_t res, arb_t res_prime, ulong n, const arb_t x, slong prec) { @@ -274,11 +322,9 @@ arb_hypgeom_legendre_p_ui(arb_t res, arb_t res_prime, ulong n, const arb_t x, sl } else if (FLINT_MIN(cost_zero, cost_one) < (1e6 * prec) * prec && n < UWORD_MAX / 4) { - mag_t t, u, err1, err2, xrad; + mag_t err1, err2, xrad; arb_t xmid; - mag_init(t); - mag_init(u); mag_init(err1); mag_init(err2); mag_init(xrad); @@ -288,14 +334,7 @@ arb_hypgeom_legendre_p_ui(arb_t res, arb_t res_prime, ulong n, const arb_t x, sl mag_zero(arb_radref(xmid)); mag_set(xrad, arb_radref(x)); - /* |P'(x)| <= min(n(n+1)/2, n/sqrt(1-x^2)) */ - arb_get_mag_lower(u, x2sub1); - mag_rsqrt(u, u); - mag_set_ui(t, n + 1); - mag_mul_2exp_si(t, t, -1); - mag_max(t, t, u); - mag_mul_ui(t, t, n); - mag_mul(err1, t, xrad); + arb_hypgeom_legendre_p_ui_deriv_bound(err1, err2, n, x, x2sub1); if (cost_zero < cost_one) arb_hypgeom_legendre_p_ui_zero(res, res_prime, n, xmid, K_zero, wp + cancellation_zero); @@ -304,32 +343,18 @@ arb_hypgeom_legendre_p_ui(arb_t res, arb_t res_prime, ulong n, const arb_t x, sl if (res != NULL) { + mag_mul(err1, err1, xrad); arb_add_error_mag(res, err1); arb_set_round(res, res, prec); } if (res_prime != NULL) { - /* |P''(x)| <= min((n+2)(n+1)n(n-1)/8, (2x|P'(x)| + n(n+1))/(1-x^2) */ - arb_get_mag_lower(u, x2sub1); - mag_set_ui(t, n); - mag_mul_ui(t, t, n + 1); - mag_div(t, t, u); - mag_add(err2, err2, t); - mag_set_ui(t, n + 2); - mag_mul_ui(t, t, n + 1); - mag_mul_ui(t, t, n); - mag_mul_ui(t, t, n - 1); - mag_mul_2exp_si(t, t, -3); - mag_min(err2, err2, t); mag_mul(err2, err2, xrad); - arb_add_error_mag(res_prime, err2); arb_set_round(res_prime, res_prime, prec); } - mag_clear(t); - mag_clear(u); mag_clear(err1); mag_clear(err2); mag_clear(xrad); diff --git a/arb_hypgeom/test/t-legendre_p_ui_deriv_bound.c b/arb_hypgeom/test/t-legendre_p_ui_deriv_bound.c new file mode 100644 index 00000000..1c15dfc0 --- /dev/null +++ b/arb_hypgeom/test/t-legendre_p_ui_deriv_bound.c @@ -0,0 +1,98 @@ +/* + 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 . +*/ + +#include "arb_hypgeom.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("legendre_p_ui_deriv_bound...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++) + { + arb_t x, x2sub1, y, z; + mag_t m1, m2, exact; + ulong n; + slong prec; + + arb_init(x); + arb_init(x2sub1); + arb_init(y); + arb_init(z); + mag_init(m1); + mag_init(m2); + mag_init(exact); + + n = n_randtest(state) % 10000; + prec = 64; + + arb_randtest(x, state, 2 * prec, 0); + mag_zero(arb_radref(x)); + while (arf_cmpabs_2exp_si(arb_midref(x), 0) >= 0) + arb_mul_2exp_si(x, x, -1); + arb_mul_2exp_si(x, x, -n_randint(state, 10)); + + arb_mul(x2sub1, x, x, 2 * prec); + arb_sub_ui(x2sub1, x2sub1, 1, 2 * prec); + arb_neg(x2sub1, x2sub1); + + arb_hypgeom_legendre_p_ui_deriv_bound(m1, m2, n, x, x2sub1); + arb_hypgeom_legendre_p_ui(y, z, n, x, prec); + arb_get_mag_lower(exact, z); + + if (mag_cmp(m1, exact) < 0) + { + flint_printf("FAIL: first derivative\n\n"); + flint_printf("n = %wu\n\n", n); + flint_printf("x = "); arb_printn(x, 50, 0); flint_printf("\n\n"); + flint_printf("z = "); arb_printn(z, 50, 0); flint_printf("\n\n"); + flint_printf("m1 = "); mag_printd(m1, 15); flint_printf("\n\n"); + flint_abort(); + } + + arb_mul(z, z, x, prec); + arb_mul_2exp_si(z, z, 1); + arb_mul_ui(y, y, n, prec); + arb_mul_ui(y, y, n + 1, prec); + arb_sub(z, z, y, prec); + arb_div(z, z, x2sub1, prec); + arb_get_mag_lower(exact, z); + + if (mag_cmp(m2, exact) < 0) + { + flint_printf("FAIL: second derivative\n\n"); + flint_printf("n = %wu\n\n", n); + flint_printf("x = "); arb_printn(x, 50, 0); flint_printf("\n\n"); + flint_printf("z = "); arb_printn(z, 50, 0); flint_printf("\n\n"); + flint_printf("m2 = "); mag_printd(m2, 15); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(x); + arb_clear(x2sub1); + arb_clear(y); + arb_clear(z); + mag_clear(m1); + mag_clear(m2); + mag_clear(exact); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/doc/source/arb_hypgeom.rst b/doc/source/arb_hypgeom.rst index 64fbb9c3..74a18695 100644 --- a/doc/source/arb_hypgeom.rst +++ b/doc/source/arb_hypgeom.rst @@ -336,6 +336,14 @@ Orthogonal polynomials and functions See :func:`acb_hypgeom_legendre_p` and :func:`acb_hypgeom_legendre_q` for definitions. +.. function:: void arb_hypgeom_legendre_p_ui_deriv_bound(mag_t dp, mag_t dp2, ulong n, const arb_t x, const arb_t x2sub1) + + Sets *dp* to an upper bound for `P'_n(x)` and *dp2* to an upper + bound for `P''_n(x)` on given *x* assumed to represent a real + number with `|x| \le 1`. The variable *x2sub1* must contain + the precomputed value `1-x^2` (or `x^2-1`). This method is used + internally to bound the propagated error for Legendre polynomials. + .. function:: void arb_hypgeom_legendre_p_ui_zero(arb_t res, arb_t res_prime, ulong n, const arb_t x, slong K, slong prec) .. function:: void arb_hypgeom_legendre_p_ui_one(arb_t res, arb_t res_prime, ulong n, const arb_t x, slong K, slong prec)