fix and refactor derivative bounds for legendre polynomials

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);

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_hypgeom_legendre_p_ui_deriv_bound(err1, err2, n, x, x2sub1);
-            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);
 
             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);

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 <http://www.gnu.org/licenses/>.
+*/
+
+#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;
+}
+

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)