From b15e14d6e74dc8f82a324ac0e5ec8209cf4173bc Mon Sep 17 00:00:00 2001
From: p15-git-acc <37548430+p15-git-acc@users.noreply.github.com>
Date: Fri, 16 Aug 2019 15:56:47 -0500
Subject: [PATCH 1/2] add csc_pi for arb_t and acb_t types
---
acb.h | 7 +++++
acb/gamma.c | 4 +--
acb/test/t-csc_pi.c | 74 +++++++++++++++++++++++++++++++++++++++++++++
arb.h | 1 +
arb/csc_pi.c | 21 +++++++++++++
arb/gamma.c | 4 +--
arb/test/t-csc_pi.c | 74 +++++++++++++++++++++++++++++++++++++++++++++
doc/source/acb.rst | 5 +++
doc/source/arb.rst | 4 +++
9 files changed, 190 insertions(+), 4 deletions(-)
create mode 100644 acb/test/t-csc_pi.c
create mode 100644 arb/csc_pi.c
create mode 100644 arb/test/t-csc_pi.c
diff --git a/acb.h b/acb.h
index 84dc347c..0d747dcc 100644
--- a/acb.h
+++ b/acb.h
@@ -723,6 +723,13 @@ 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);
+ACB_INLINE void
+acb_csc_pi(acb_t y, const acb_t x, slong prec)
+{
+ acb_sin_pi(y, x, prec);
+ acb_inv(y, y, 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/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 fb80313a..0d7dd2c1 100644
--- a/arb.h
+++ b/arb.h
@@ -510,6 +510,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 ba59d171..6696dcf9 100644
--- a/doc/source/acb.rst
+++ b/doc/source/acb.rst
@@ -716,6 +716,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 6598e1bc..6d65f98f 100644
--- a/doc/source/arb.rst
+++ b/doc/source/arb.rst
@@ -1053,6 +1053,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`.
From 2e42836203ebef0c75dbc561ba0a8989f104206c Mon Sep 17 00:00:00 2001
From: p15-git-acc <37548430+p15-git-acc@users.noreply.github.com>
Date: Thu, 10 Sep 2020 11:46:52 -0500
Subject: [PATCH 2/2] avoid inverting sine in acb_csc_pi
---
acb.h | 8 +-----
acb/csc_pi.c | 70 ++++++++++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 71 insertions(+), 7 deletions(-)
create mode 100644 acb/csc_pi.c
diff --git a/acb.h b/acb.h
index 0d747dcc..3dc90a8e 100644
--- a/acb.h
+++ b/acb.h
@@ -722,13 +722,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);
-
-ACB_INLINE void
-acb_csc_pi(acb_t y, const acb_t x, slong prec)
-{
- acb_sin_pi(y, x, prec);
- acb_inv(y, y, 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);
+ }
+ }
+}
+