diff --git a/doc/source/fmprb.rst b/doc/source/fmprb.rst index 7f451365..6e6f9907 100644 --- a/doc/source/fmprb.rst +++ b/doc/source/fmprb.rst @@ -480,6 +480,24 @@ Elementary functions Sets `s = \sin \pi x`, `c = \cos \pi x`. +.. function:: void fmprb_sin_pi_fmpq(fmprb_t s, const fmpq_t x, long prec) + +.. function:: void fmprb_cos_pi_fmpq(fmprb_t c, const fmpq_t x, long prec) + +.. function:: void fmprb_sin_cos_pi_fmpq(fmprb_t s, fmprb_t c, const fmpq_t x, long prec) + + Sets `s = \sin \pi x`, `c = \cos \pi x` where `x` is a rational + number (whose numerator and denominator are assumed to be reduced). + We first use trigonometric symmetries to reduce the argument to the + octant `[0, 1/4]`. If the reduced argument is one of + `0, 1/4, 1/6, 1/8, 1/5, 1/10`, we use an explicit algebraic + formula. Otherwise + we multiply by a numerical approximation of `\pi` + and evaluate the trigonometric function the usual way. + Since the argument has been reduced to the first octant, this + gives full accuracy even if the original argument is close to + some root other the origin. + .. function:: void fmprb_atan(fmprb_t z, const fmprb_t x, long prec) Sets `z = \tan^{-1} x`. Letting `d = \max(0, |m| - r)`, diff --git a/fmprb.h b/fmprb.h index afbb88e7..38b35b56 100644 --- a/fmprb.h +++ b/fmprb.h @@ -275,6 +275,10 @@ void fmprb_sin_pi(fmprb_t s, const fmprb_t x, long prec); void fmprb_cos_pi(fmprb_t c, const fmprb_t x, long prec); void fmprb_sin_cos_pi(fmprb_t s, fmprb_t c, const fmprb_t x, long prec); +void fmprb_sin_cos_pi_fmpq(fmprb_t s, fmprb_t c, const fmpq_t x, long prec); +void fmprb_sin_pi_fmpq(fmprb_t s, const fmpq_t x, long prec); +void fmprb_cos_pi_fmpq(fmprb_t c, const fmpq_t x, long prec); + void fmprb_sinh(fmprb_t z, const fmprb_t x, long prec); void fmprb_cosh(fmprb_t z, const fmprb_t x, long prec); void fmprb_sinh_cosh(fmprb_t s, fmprb_t c, const fmprb_t x, long prec); diff --git a/fmprb/sin_cos_pi_fmpq.c b/fmprb/sin_cos_pi_fmpq.c new file mode 100644 index 00000000..942832e9 --- /dev/null +++ b/fmprb/sin_cos_pi_fmpq.c @@ -0,0 +1,289 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmprb.h" + +/* evaluates sin(pi v/w), cos(pi v/w) assuming + (for best performance and accuracy) + v, w reduced and 0 <= v/w <= 1/4. */ + +void +_fmprb_sin_cos_pi_fmpq_oct(fmprb_t s, fmprb_t c, + const fmpz_t v, const fmpz_t w, long prec) +{ + if (*v == 0) + { + fmprb_zero(s); + fmprb_one(c); + } + else if (*v == 1 && *w == 4) + { + fmprb_sqrt_ui(s, 2, prec); + fmprb_mul_2exp_si(s, s, -1); + fmprb_set(c, s); + } + else if (*v == 1 && *w == 6) + { + fmprb_one(s); + fmprb_mul_2exp_si(s, s, -1); + fmprb_sqrt_ui(c, 3, prec); + fmprb_mul_2exp_si(c, c, -1); + } + else if (*v == 1 && *w == 8) + { + fmprb_sqrt_ui(s, 2, prec); + fmprb_add_ui(c, s, 2, prec); + fmprb_sqrt(c, c, prec); + fmprb_mul_2exp_si(c, c, -1); + fmprb_sub_ui(s, s, 2, prec); + fmprb_neg(s, s); + fmprb_sqrt(s, s, prec); + fmprb_mul_2exp_si(s, s, -1); + } + else if (*v == 1 && *w == 5) + { + fmprb_sqrt_ui(s, 5, prec); + fmprb_add_ui(c, s, 1, prec); + fmprb_mul_2exp_si(c, c, -2); + fmprb_sub_ui(s, s, 5, prec); + fmprb_neg(s, s); + fmprb_mul_2exp_si(s, s, -3); + fmprb_sqrt(s, s, prec); + } + else if (*v == 1 && *w == 10) + { + fmprb_sqrt_ui(c, 5, prec); + fmprb_sub_ui(s, c, 1, prec); + fmprb_mul_2exp_si(s, s, -2); + fmprb_add_ui(c, c, 5, prec); + fmprb_mul_2exp_si(c, c, -3); + fmprb_sqrt(c, c, prec); + } + else + { + fmprb_const_pi(s, prec); + fmprb_mul_fmpz(s, s, v, prec); + fmprb_div_fmpz(s, s, w, prec); + fmprb_sin_cos(s, c, s, prec); + } +} + +void +_fmprb_sin_pi_fmpq_oct(fmprb_t s, const fmpz_t v, const fmpz_t w, long prec) +{ + if (*v == 0) + { + fmprb_zero(s); + } + else if (*v == 1 && *w == 4) + { + fmprb_sqrt_ui(s, 2, prec); + fmprb_mul_2exp_si(s, s, -1); + } + else if (*v == 1 && *w == 6) + { + fmprb_one(s); + fmprb_mul_2exp_si(s, s, -1); + } + else if (*v == 1 && *w == 8) + { + fmprb_sqrt_ui(s, 2, prec); + fmprb_sub_ui(s, s, 2, prec); + fmprb_neg(s, s); + fmprb_sqrt(s, s, prec); + fmprb_mul_2exp_si(s, s, -1); + } + else if (*v == 1 && *w == 5) + { + fmprb_sqrt_ui(s, 5, prec); + fmprb_sub_ui(s, s, 5, prec); + fmprb_neg(s, s); + fmprb_mul_2exp_si(s, s, -3); + fmprb_sqrt(s, s, prec); + } + else if (*v == 1 && *w == 10) + { + fmprb_sqrt_ui(s, 5, prec); + fmprb_sub_ui(s, s, 1, prec); + fmprb_mul_2exp_si(s, s, -2); + } + else + { + fmprb_const_pi(s, prec); + fmprb_mul_fmpz(s, s, v, prec); + fmprb_div_fmpz(s, s, w, prec); + fmprb_sin(s, s, prec); + } +} + +void +_fmprb_cos_pi_fmpq_oct(fmprb_t c, const fmpz_t v, const fmpz_t w, long prec) +{ + if (*v == 0) + { + fmprb_one(c); + } + else if (*v == 1 && *w == 4) + { + fmprb_sqrt_ui(c, 2, prec); + fmprb_mul_2exp_si(c, c, -1); + } + else if (*v == 1 && *w == 6) + { + fmprb_sqrt_ui(c, 3, prec); + fmprb_mul_2exp_si(c, c, -1); + } + else if (*v == 1 && *w == 8) + { + fmprb_sqrt_ui(c, 2, prec); + fmprb_add_ui(c, c, 2, prec); + fmprb_sqrt(c, c, prec); + fmprb_mul_2exp_si(c, c, -1); + } + else if (*v == 1 && *w == 5) + { + fmprb_sqrt_ui(c, 5, prec); + fmprb_add_ui(c, c, 1, prec); + fmprb_mul_2exp_si(c, c, -2); + } + else if (*v == 1 && *w == 10) + { + fmprb_sqrt_ui(c, 5, prec); + fmprb_add_ui(c, c, 5, prec); + fmprb_mul_2exp_si(c, c, -3); + fmprb_sqrt(c, c, prec); + } + else + { + fmprb_const_pi(c, prec); + fmprb_mul_fmpz(c, c, v, prec); + fmprb_div_fmpz(c, c, w, prec); + fmprb_cos(c, c, prec); + } +} + +unsigned int +reduce_octant(fmpz_t v, fmpz_t w, const fmpq_t x) +{ + const fmpz * p = fmpq_numref(x); + const fmpz * q = fmpq_denref(x); + fmpz_t u; + unsigned int octant; + mp_bitcnt_t vval, wval; + + fmpz_init(u); + + fmpz_mul_2exp(w, p, 2); + fmpz_fdiv_qr(w, v, w, q); + octant = fmpz_fdiv_ui(w, 8); + fmpz_mul_2exp(w, q, 2); + + if (octant % 2 != 0) + fmpz_sub(v, q, v); + + vval = fmpz_val2(v); + wval = fmpz_val2(w); + vval = FLINT_MIN(vval, wval); + + if (vval != 0) + { + fmpz_tdiv_q_2exp(v, v, vval); + fmpz_tdiv_q_2exp(w, w, vval); + } + + fmpz_clear(u); + + return octant; +} + +void +fmprb_sin_cos_pi_fmpq(fmprb_t s, fmprb_t c, const fmpq_t x, long prec) +{ + fmpz_t v, w; + unsigned int octant; + + fmpz_init(v); + fmpz_init(w); + + octant = reduce_octant(v, w, x); + + if ((octant + 1) % 4 < 2) + _fmprb_sin_cos_pi_fmpq_oct(s, c, v, w, prec); + else + _fmprb_sin_cos_pi_fmpq_oct(c, s, v, w, prec); + + if ((octant + 6) % 8 < 4) fmprb_neg(c, c); + if (octant >= 4) fmprb_neg(s, s); + + fmpz_clear(v); + fmpz_clear(w); +} + +void +fmprb_sin_pi_fmpq(fmprb_t s, const fmpq_t x, long prec) +{ + fmpz_t v, w; + unsigned int octant; + + fmpz_init(v); + fmpz_init(w); + + octant = reduce_octant(v, w, x); + + if ((octant + 1) % 4 < 2) + _fmprb_sin_pi_fmpq_oct(s, v, w, prec); + else + _fmprb_cos_pi_fmpq_oct(s, v, w, prec); + + if (octant >= 4) + fmprb_neg(s, s); + + fmpz_clear(v); + fmpz_clear(w); +} + +void +fmprb_cos_pi_fmpq(fmprb_t c, const fmpq_t x, long prec) +{ + fmpz_t v, w; + unsigned int octant; + + fmpz_init(v); + fmpz_init(w); + + octant = reduce_octant(v, w, x); + + if ((octant + 1) % 4 < 2) + _fmprb_cos_pi_fmpq_oct(c, v, w, prec); + else + _fmprb_sin_pi_fmpq_oct(c, v, w, prec); + + if ((octant + 6) % 8 < 4) + fmprb_neg(c, c); + + fmpz_clear(v); + fmpz_clear(w); +} + diff --git a/fmprb/test/t-cos_pi_fmpq.c b/fmprb/test/t-cos_pi_fmpq.c new file mode 100644 index 00000000..5c0de92b --- /dev/null +++ b/fmprb/test/t-cos_pi_fmpq.c @@ -0,0 +1,78 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmprb.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("cos_pi_fmpq...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + fmprb_t c1, c2; + fmpq_t x; + long prec; + + prec = 2 + n_randint(state, 1000); + + fmprb_init(c1); + fmprb_init(c2); + fmpq_init(x); + + fmpq_randtest(x, state, 1 + n_randint(state, 200)); + + fmprb_cos_pi_fmpq(c1, x, prec); + + fmprb_const_pi(c2, prec); + fmprb_mul_fmpz(c2, c2, fmpq_numref(x), prec); + fmprb_div_fmpz(c2, c2, fmpq_denref(x), prec); + fmprb_cos(c2, c2, prec); + + if (!fmprb_overlaps(c1, c2)) + { + printf("FAIL: overlap\n\n"); + printf("x = "); fmpq_print(x); printf("\n\n"); + printf("c1 = "); fmprb_printd(c1, 15); printf("\n\n"); + printf("c2 = "); fmprb_printd(c2, 15); printf("\n\n"); + abort(); + } + + fmprb_clear(c1); + fmprb_clear(c2); + fmpq_clear(x); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/fmprb/test/t-sin_cos_pi_fmpq.c b/fmprb/test/t-sin_cos_pi_fmpq.c new file mode 100644 index 00000000..3a53875d --- /dev/null +++ b/fmprb/test/t-sin_cos_pi_fmpq.c @@ -0,0 +1,84 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmprb.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("sin_cos_pi_fmpq...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + fmprb_t s1, c1, s2, c2; + fmpq_t x; + long prec; + + prec = 2 + n_randint(state, 1000); + + fmprb_init(s1); + fmprb_init(c1); + fmprb_init(s2); + fmprb_init(c2); + fmpq_init(x); + + fmpq_randtest(x, state, 1 + n_randint(state, 200)); + + fmprb_sin_cos_pi_fmpq(s1, c1, x, prec); + + fmprb_const_pi(s2, prec); + fmprb_mul_fmpz(s2, s2, fmpq_numref(x), prec); + fmprb_div_fmpz(s2, s2, fmpq_denref(x), prec); + fmprb_sin_cos(s2, c2, s2, prec); + + if (!fmprb_overlaps(s1, s2) || !fmprb_overlaps(c1, c2)) + { + printf("FAIL: overlap\n\n"); + printf("x = "); fmpq_print(x); printf("\n\n"); + printf("s1 = "); fmprb_printd(s1, 15); printf("\n\n"); + printf("c1 = "); fmprb_printd(c1, 15); printf("\n\n"); + printf("s2 = "); fmprb_printd(s2, 15); printf("\n\n"); + printf("c2 = "); fmprb_printd(c2, 15); printf("\n\n"); + abort(); + } + + fmprb_clear(s1); + fmprb_clear(c1); + fmprb_clear(s2); + fmprb_clear(c2); + fmpq_clear(x); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/fmprb/test/t-sin_pi_fmpq.c b/fmprb/test/t-sin_pi_fmpq.c new file mode 100644 index 00000000..0b418d8f --- /dev/null +++ b/fmprb/test/t-sin_pi_fmpq.c @@ -0,0 +1,78 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmprb.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("sin_pi_fmpq...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + fmprb_t s1, s2; + fmpq_t x; + long prec; + + prec = 2 + n_randint(state, 1000); + + fmprb_init(s1); + fmprb_init(s2); + fmpq_init(x); + + fmpq_randtest(x, state, 1 + n_randint(state, 200)); + + fmprb_sin_pi_fmpq(s1, x, prec); + + fmprb_const_pi(s2, prec); + fmprb_mul_fmpz(s2, s2, fmpq_numref(x), prec); + fmprb_div_fmpz(s2, s2, fmpq_denref(x), prec); + fmprb_sin(s2, s2, prec); + + if (!fmprb_overlaps(s1, s2)) + { + printf("FAIL: overlap\n\n"); + printf("x = "); fmpq_print(x); printf("\n\n"); + printf("s1 = "); fmprb_printd(s1, 15); printf("\n\n"); + printf("s2 = "); fmprb_printd(s2, 15); printf("\n\n"); + abort(); + } + + fmprb_clear(s1); + fmprb_clear(s2); + fmpq_clear(x); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} +