Merge pull request #291 from p15-git-acc/csc_pi

add csc_pi for arb_t and acb_t types

acb.h

@@ -729,6 +729,7 @@ void acb_cos_pi(acb_t r, const acb_t z, slong prec);
 void acb_sin_cos_pi(acb_t s, acb_t c, const acb_t z, slong prec);
 void acb_tan_pi(acb_t r, const acb_t z, slong prec);
 void acb_cot_pi(acb_t r, const acb_t z, slong prec);
+void acb_csc_pi(acb_t y, const acb_t x, slong prec);
 
 void acb_sinc(acb_t res, const acb_t z, slong prec);
 void acb_sinc_pi(acb_t res, const acb_t z, slong prec);

acb/csc_pi.c

@@ -0,0 +1,70 @@
+/*
+    Copyright (C) 2015 Fredrik Johansson
+    Copyright (C) 2020 D.H.J. Polymath
+
+    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_csc_pi(acb_t res, const acb_t z, slong prec)
+{
+    if (acb_contains_zero(z) || !acb_is_finite(z))
+    {
+        acb_indeterminate(res);
+    }
+    else if (arb_is_zero(acb_imagref(z)))
+    {
+        arb_csc_pi(acb_realref(res), acb_realref(z), prec);
+        arb_zero(acb_imagref(res));
+    }
+    else if (arb_is_zero(acb_realref(z)))
+    {
+        arb_const_pi(acb_realref(res), prec);
+        arb_mul(acb_imagref(res), acb_imagref(z), acb_realref(res), prec);
+        arb_csch(acb_imagref(res), acb_imagref(res), prec);
+        arb_neg(acb_imagref(res), acb_imagref(res));
+        arb_zero(acb_realref(res));
+    }
+    else
+    {
+        if (arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 0) > 0)
+        {
+            acb_t t;
+            acb_init(t);
+
+            if (arf_sgn(arb_midref(acb_imagref(z))) < 0)
+            {
+                acb_neg(t, z);
+                acb_exp_pi_i(t, t, prec + 4);
+                acb_mul(res, t, t, prec + 4);
+                acb_sub_ui(res, res, 1, prec + 4);
+                acb_div(res, t, res, prec);
+                acb_neg(res, res);
+            }
+            else
+            {
+                acb_exp_pi_i(t, z, prec + 4);
+                acb_mul(res, t, t, prec + 4);
+                acb_sub_ui(res, res, 1, prec + 4);
+                acb_div(res, t, res, prec);
+            }
+
+            acb_mul_2exp_si(res, res, 1);
+            acb_mul_onei(res, res);
+            acb_clear(t);
+        }
+        else
+        {
+            acb_sin_pi(res, z, prec + 4);
+            acb_inv(res, res, prec);
+        }
+    }
+}
+

acb/gamma.c

@@ -180,8 +180,8 @@ _acb_gamma(acb_t y, const acb_t x, slong prec, int inverse)
             /* gamma(x) = (rf(1-x, r) * pi) rgamma(1-x+r) csc(pi x) */
             acb_neg(v, v);
             acb_exp(v, v, wp);
-            acb_sin_pi(t, x, wp);   /* todo: write a csc_pi function */
-            acb_div(v, v, t, wp);
+            acb_csc_pi(t, x, wp);
+            acb_mul(v, v, t, wp);
             acb_mul(y, v, u, prec);
         }
     }

acb/test/t-csc_pi.c

@@ -0,0 +1,74 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+    Copyright (C) 2019 D.H.J. Polymath
+
+    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("csc_pi....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        acb_t x, a, b;
+        slong prec1, prec2;
+
+        prec1 = 2 + n_randint(state, 200);
+        prec2 = prec1 + 30;
+
+        acb_init(x);
+        acb_init(a);
+        acb_init(b);
+
+        acb_randtest_special(x, state, 1 + n_randint(state, 300), 100);
+        acb_randtest_special(a, state, 1 + n_randint(state, 300), 100);
+        acb_randtest_special(b, state, 1 + n_randint(state, 300), 100);
+
+        if (n_randint(state, 2))
+        {
+            acb_csc_pi(a, x, prec1);
+        }
+        else
+        {
+            acb_set(a, x);
+            acb_csc_pi(a, a, prec1);
+        }
+
+        acb_sin_pi(b, x, prec2);
+        acb_inv(b, b, prec2);
+
+        /* check consistency */
+        if (!acb_overlaps(a, b))
+        {
+            flint_printf("FAIL: overlap\n\n");
+            flint_printf("x = "); acb_printn(x, 20, 0); flint_printf("\n\n");
+            flint_printf("a = "); acb_printn(a, 20, 0); flint_printf("\n\n");
+            flint_printf("b = "); acb_printn(b, 20, 0); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        acb_clear(x);
+        acb_clear(a);
+        acb_clear(b);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

arb.h

@@ -520,6 +520,7 @@ void arb_acosh(arb_t z, const arb_t x, slong prec);
 
 void arb_sec(arb_t res, const arb_t x, slong prec);
 void arb_csc(arb_t res, const arb_t x, slong prec);
+void arb_csc_pi(arb_t res, const arb_t x, slong prec);
 void arb_sech(arb_t res, const arb_t x, slong prec);
 void arb_csch(arb_t res, const arb_t x, slong prec);
 

arb/csc_pi.c

@@ -0,0 +1,21 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+    Copyright (C) 2017 D.H.J. Polymath
+
+    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.h"
+
+void
+arb_csc_pi(arb_t res, const arb_t x, slong prec)
+{
+    arb_sin_pi(res, x, prec + 4);
+    arb_inv(res, res, prec);
+}
+

arb/gamma.c

@@ -860,8 +860,8 @@ _arb_gamma(arb_t y, const arb_t x, slong prec, int inverse)
             /* gamma(x) = (rf(1-x, r) * pi) rgamma(1-x+r) csc(pi x) */
             arb_neg(v, v);
             arb_exp(v, v, wp);
-            arb_sin_pi(t, x, wp);   /* todo: write a csc_pi function */
-            arb_div(v, v, t, wp);
+            arb_csc_pi(t, x, wp);
+            arb_mul(v, v, t, wp);
             arb_mul(y, v, u, prec);
         }
     }

arb/test/t-csc_pi.c

@@ -0,0 +1,74 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+    Copyright (C) 2019 D.H.J. Polymath
+
+    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.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("csc_pi....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        arb_t x, a, b;
+        slong prec1, prec2;
+
+        prec1 = 2 + n_randint(state, 200);
+        prec2 = prec1 + 30;
+
+        arb_init(x);
+        arb_init(a);
+        arb_init(b);
+
+        arb_randtest_special(x, state, 1 + n_randint(state, 300), 100);
+        arb_randtest_special(a, state, 1 + n_randint(state, 300), 100);
+        arb_randtest_special(b, state, 1 + n_randint(state, 300), 100);
+
+        if (n_randint(state, 2))
+        {
+            arb_csc_pi(a, x, prec1);
+        }
+        else
+        {
+            arb_set(a, x);
+            arb_csc_pi(a, a, prec1);
+        }
+
+        arb_sin_pi(b, x, prec2);
+        arb_inv(b, b, prec2);
+
+        /* check consistency */
+        if (!arb_overlaps(a, b))
+        {
+            flint_printf("FAIL: overlap\n\n");
+            flint_printf("x = "); arb_printn(x, 20, 0); flint_printf("\n\n");
+            flint_printf("a = "); arb_printn(a, 20, 0); flint_printf("\n\n");
+            flint_printf("b = "); arb_printn(b, 20, 0); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        arb_clear(x);
+        arb_clear(a);
+        arb_clear(b);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

doc/source/acb.rst

@@ -726,6 +726,11 @@ Trigonometric functions
 
     Computes `\csc(x) = 1 / \sin(z)`.
 
+.. function:: void acb_csc_pi(acb_t res, const acb_t z, slong prec)
+
+    Computes `\csc(\pi x) = 1 / \sin(\pi z)`. Evaluates the sine accurately
+    via :func:`acb_sin_pi`.
+
 .. function:: void acb_sinc(acb_t s, const acb_t z, slong prec)
 
     Sets `s = \operatorname{sinc}(x) = \sin(z) / z`.

doc/source/arb.rst

@@ -1109,6 +1109,10 @@ Trigonometric functions
 
     Computes `\csc(x) = 1 / \sin(x)`.
 
+.. function:: void arb_csc_pi(arb_t res, const arb_t x, slong prec)
+
+    Computes `\csc(\pi x) = 1 / \sin(\pi x)`.
+
 .. function:: void arb_sinc(arb_t z, const arb_t x, slong prec)
 
     Sets `z = \operatorname{sinc}(x) = \sin(x) / x`.