diff --git a/acb.h b/acb.h index 54cbf12f..c0a3068d 100644 --- a/acb.h +++ b/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); diff --git a/acb/csc_pi.c b/acb/csc_pi.c new file mode 100644 index 00000000..1d8c95f1 --- /dev/null +++ b/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 . +*/ + +#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); + } + } +} + diff --git a/acb/gamma.c b/acb/gamma.c index 21f707f5..21d62a90 100644 --- a/acb/gamma.c +++ b/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); } } diff --git a/acb/test/t-csc_pi.c b/acb/test/t-csc_pi.c new file mode 100644 index 00000000..78abf586 --- /dev/null +++ b/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 . +*/ + +#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; +} + diff --git a/arb.h b/arb.h index 3b757135..d82aaffd 100644 --- a/arb.h +++ b/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); diff --git a/arb/csc_pi.c b/arb/csc_pi.c new file mode 100644 index 00000000..3690d55c --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb/gamma.c b/arb/gamma.c index da531814..eb3cad3b 100644 --- a/arb/gamma.c +++ b/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); } } diff --git a/arb/test/t-csc_pi.c b/arb/test/t-csc_pi.c new file mode 100644 index 00000000..d7c11af5 --- /dev/null +++ b/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 . +*/ + +#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; +} + diff --git a/doc/source/acb.rst b/doc/source/acb.rst index a083dec1..d012bb0a 100644 --- a/doc/source/acb.rst +++ b/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`. diff --git a/doc/source/arb.rst b/doc/source/arb.rst index feea1222..dc8f7eb1 100644 --- a/doc/source/arb.rst +++ b/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`.