From 81236fa0937c29ba0b7f69849c0430ac4da99a87 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Tue, 17 Apr 2012 12:45:59 +0200 Subject: [PATCH] enhancements to log/exp series; general zeta_ui; series expansion of log gamma --- arb.h | 3 + arb/log_ui.c | 2 + arb/test/t-zeta_ui.c | 81 ++++++++++++++++++++++++++ arb/zeta_ui.c | 104 ++++++++++++++++++++++++++++++++++ arb/zeta_ui_mpfr.c | 42 ++++++++++++++ arb_poly.h | 10 +++- arb_poly/exp_series.c | 34 ++++++----- arb_poly/gamma_log_series.c | 68 ++++++++++++++++++++++ arb_poly/log_series.c | 56 +++++------------- arb_poly/set_coeff_same_exp.c | 63 ++++++++++++++++++++ 10 files changed, 403 insertions(+), 60 deletions(-) create mode 100644 arb/test/t-zeta_ui.c create mode 100644 arb/zeta_ui.c create mode 100644 arb/zeta_ui_mpfr.c create mode 100644 arb_poly/gamma_log_series.c create mode 100644 arb_poly/set_coeff_same_exp.c diff --git a/arb.h b/arb.h index bb955b48..f27e4304 100644 --- a/arb.h +++ b/arb.h @@ -148,7 +148,10 @@ void _arb_rad_add_ufloat(arb_t y, const ufloat_t err); void arb_const_pi_chudnovsky(arb_t x); void arb_const_euler_brent_mcmillan(arb_t x); void arb_const_zeta3_bsplit(arb_t x); + void arb_zeta_ui_bsplit(arb_t x, ulong s); +void arb_zeta_ui_mpfr(arb_t x, ulong n); +void arb_zeta_ui(arb_t x, ulong n); /* fmpz extras */ diff --git a/arb/log_ui.c b/arb/log_ui.c index 0db5ab49..38f04194 100644 --- a/arb/log_ui.c +++ b/arb/log_ui.c @@ -63,4 +63,6 @@ arb_log_ui(arb_t x, ulong n) } arb_set_mpfr(x, t, 1); + + mpfr_clear(t); } diff --git a/arb/test/t-zeta_ui.c b/arb/test/t-zeta_ui.c new file mode 100644 index 00000000..5d0fc05a --- /dev/null +++ b/arb/test/t-zeta_ui.c @@ -0,0 +1,81 @@ +/*============================================================================= + + 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 "arb.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("zeta_ui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + arb_t r; + ulong n; + mpfr_t s; + long effective_prec; + + arb_init(r, 1 + n_randint(state, 1 << n_randint(state, 14))); + mpfr_init2(s, arb_prec(r) + 100); + arb_randtest(r, state, 10); + + do { n = n_randint(state, 1 << n_randint(state, 10)); } while (n == 1); + + arb_zeta_ui(r, n); + mpfr_zeta_ui(s, n, MPFR_RNDN); + + if (!arb_contains_mpfr(r, s)) + { + printf("FAIL: containment\n\n"); + printf("n = %lu\n\n", n); + printf("r = "); arb_debug(r); printf("\n\n"); + printf("s = "); mpfr_printf("%.275Rf\n", s); printf("\n\n"); + abort(); + } + + effective_prec = fmpz_bits(arb_midref(r)) - fmpz_bits(arb_radref(r)); + + if (!fmpz_is_zero(arb_radref(r)) && effective_prec < arb_prec(r) - 4) + { + printf("FAIL: poor accuracy\n\n"); + printf("r = "); arb_debug(r); printf("\n\n"); + abort(); + } + + arb_clear(r); + mpfr_clear(s); + } + + flint_randclear(state); + _fmpz_cleanup(); + mpfr_free_cache(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arb/zeta_ui.c b/arb/zeta_ui.c new file mode 100644 index 00000000..c1b2f46f --- /dev/null +++ b/arb/zeta_ui.c @@ -0,0 +1,104 @@ +/*============================================================================= + + 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 "arb.h" +#include "arith.h" + +static void +mpfr_zeta_ui_bernoulli(mpfr_t x, ulong n) +{ + fmpq_t b; + mpfr_t t; + mpz_t f; + + if (n % 2) + abort(); + + fmpq_init(b); + mpfr_init2(t, mpfr_get_prec(x) + FLINT_BIT_COUNT(n) + 2); + mpz_init(f); + + bernoulli_number(b, n); + + fmpq_get_mpfr(x, b, MPFR_RNDD); + mpfr_const_pi(t, MPFR_RNDD); + mpfr_mul_2exp(t, t, 1, MPFR_RNDD); + mpfr_pow_ui(t, t, n, MPFR_RNDD); + mpz_fac_ui(f, n); + mpfr_div_z(t, t, f, MPFR_RNDD); + mpfr_mul(x, x, t, MPFR_RNDD); + mpfr_abs(x, x, MPFR_RNDD); + mpfr_div_2exp(x, x, 1, MPFR_RNDD); + + mpfr_clear(t); + fmpq_clear(b); + mpz_clear(f); +} + +void +arb_zeta_ui(arb_t x, ulong n) +{ + long prec = arb_prec(x); + + /* asymptotic case */ + if (n == 0 || n > 0.7 * prec) + { + arb_zeta_ui_mpfr(x, n); + } + else if (n == 3) + { + arb_const_zeta3_bsplit(x); + } + /* small even n */ + else if ((n % 2 == 0) && (n < 40 + 0.11*prec)) + { + mpfr_t t; + mpfr_init2(t, prec + 10); + mpfr_zeta_ui_bernoulli(t, n); + arb_set_mpfr(x, t, 100); /* XXX */ + mpfr_clear(t); + } + /* small odd n, extremely high precision */ + /* FIXME: arb: n < prec * 0.0006 */ + else if (n < prec * 0.0006) + { + arb_zeta_ui_bsplit(x, n); + } + /* large n */ + else if (prec > 20 && n > 0.4 * pow(prec, 0.8)) + { + mpfr_t t; + mpfr_init2(t, prec + 10); + _zeta_inv_euler_product(t, n, 0); + mpfr_si_div(t, 1, t, MPFR_RNDD); + arb_set_mpfr(x, t, 100); /* XXX */ + mpfr_clear(t); + } + /* fallback */ + else + { + arb_zeta_ui_mpfr(x, n); + } +} diff --git a/arb/zeta_ui_mpfr.c b/arb/zeta_ui_mpfr.c new file mode 100644 index 00000000..0d67dca6 --- /dev/null +++ b/arb/zeta_ui_mpfr.c @@ -0,0 +1,42 @@ +/*============================================================================= + + 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 "arb.h" + +void +arb_zeta_ui_mpfr(arb_t x, ulong n) +{ + mpfr_t t; + + mpfr_init2(t, FLINT_MAX(arb_prec(x), FLINT_BITS)); + + if (n == 1) + mpfr_const_euler(t, MPFR_RNDN); + else + mpfr_zeta_ui(t, n, MPFR_RNDN); + + arb_set_mpfr(x, t, 1); + mpfr_clear(t); +} diff --git a/arb_poly.h b/arb_poly.h index 083c3b8c..059327cf 100644 --- a/arb_poly.h +++ b/arb_poly.h @@ -113,7 +113,13 @@ void arb_poly_inv_series(arb_poly_t z, const arb_poly_t x, long n); void arb_poly_derivative(arb_poly_t z, const arb_poly_t x); void arb_poly_integral(arb_poly_t z, const arb_poly_t x); -void arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n); -void arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n); +void arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n, int exact_const); +void arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n, int exact_const); + +void _arb_poly_set_coeff_same_exp(arb_poly_t z, long i, const arb_t c); + + +void arb_poly_gamma_log_series(arb_poly_t z, ulong n); + #endif diff --git a/arb_poly/exp_series.c b/arb_poly/exp_series.c index 1c614ba0..d45dbc06 100644 --- a/arb_poly/exp_series.c +++ b/arb_poly/exp_series.c @@ -26,7 +26,7 @@ #include "arb_poly.h" void -arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n) +arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n, int exact_const) { long a[FLINT_BITS]; long i; @@ -52,19 +52,22 @@ arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n) arb_poly_init(c, arb_poly_prec(z)); - /* first coefficient */ - arb_init(r, arb_poly_prec(z)); - fmpz_set(arb_midref(r), arb_poly_coeffs(x)); - fmpz_set(arb_radref(r), arb_poly_radref(x)); - fmpz_set(arb_expref(r), arb_poly_expref(x)); + if (!exact_const) + { + /* first coefficient */ + arb_init(r, arb_poly_prec(z)); + fmpz_set(arb_midref(r), arb_poly_coeffs(x)); + fmpz_set(arb_radref(r), arb_poly_radref(x)); + fmpz_set(arb_expref(r), arb_poly_expref(x)); - arb_exp(r, r); - _arb_poly_fit_length(c, 1); - fmpz_set(arb_poly_coeffs(c), arb_midref(r)); - fmpz_set(arb_poly_radref(c), arb_radref(r)); - fmpz_set(arb_poly_expref(c), arb_expref(r)); - c->length = 1; - arb_clear(r); + arb_exp(r, r); + _arb_poly_fit_length(c, 1); + fmpz_set(arb_poly_coeffs(c), arb_midref(r)); + fmpz_set(arb_poly_radref(c), arb_radref(r)); + fmpz_set(arb_poly_expref(c), arb_expref(r)); + c->length = 1; + arb_clear(r); + } arb_poly_set_si(z, 1); @@ -72,14 +75,15 @@ arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n) { n = a[i]; - arb_poly_log_series(t, z, n); + arb_poly_log_series(t, z, n, 1); arb_poly_neg(t, t); arb_poly_add(t, t, v); arb_poly_mullow(t, z, t, n); arb_poly_add(z, z, t); } - arb_poly_mul(z, z, c); + if (!exact_const) + arb_poly_mul(z, z, c); arb_poly_clear(t); arb_poly_clear(v); diff --git a/arb_poly/gamma_log_series.c b/arb_poly/gamma_log_series.c new file mode 100644 index 00000000..e182da00 --- /dev/null +++ b/arb_poly/gamma_log_series.c @@ -0,0 +1,68 @@ +/*============================================================================= + + This file is part of FLINT. + + FLINT 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. + + FLINT 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 FLINT; 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 "arb_poly.h" + +/* series expansion for log(gamma(1-x)) at x = 0 */ + +void +arb_poly_gamma_log_series(arb_poly_t z, ulong n) +{ + arb_t c; + long i, prec; + + prec = arb_poly_prec(z) + FLINT_BIT_COUNT(n); + + arb_init(c, prec); + + arb_poly_zero(z); + _arb_poly_fit_length(z, n); + z->length = n; + + fmpz_set_si(arb_poly_expref(z), -prec); + + for (i = 0; i < n; i++) + { + if (i == 0) + { + arb_zero(c); + } + else if (i == 1) + { + arb_const_euler_brent_mcmillan(c); + } + else + { + /* printf("%ld of %ld\n", i, n); */ + arb_zeta_ui(c, i); + fmpz_tdiv_q_ui(arb_midref(c), arb_midref(c), i); + fmpz_add_ui(arb_radref(c), arb_radref(c), 1); + } + + _arb_poly_set_coeff_same_exp(z, i, c); + } + + arb_clear(c); +} diff --git a/arb_poly/log_series.c b/arb_poly/log_series.c index cfb626f8..f4677f4b 100644 --- a/arb_poly/log_series.c +++ b/arb_poly/log_series.c @@ -26,43 +26,7 @@ #include "arb_poly.h" void -_arb_poly_set_coeff_same_exp(arb_poly_t z, const arb_t c) -{ - long e1, e2; - fmpz_t rad; - - e1 = z->exp; - e2 = c->exp; - - if (z->length < 1) - { - _arb_poly_fit_length(z, 1); - fmpz_set(arb_expref(z), arb_expref(c)); - z->length = 1; - } - - fmpz_init(rad); - - if (e1 > e2) - { - fmpz_tdiv_q_2exp(z->coeffs, arb_midref(c), e1 - e2); - fmpz_cdiv_q_2exp(rad, arb_radref(c), e1 - e2); - fmpz_add_ui(rad, rad, 1); - } - else - { - fmpz_mul_2exp(z->coeffs, arb_midref(c), e2 - e1); - fmpz_mul_2exp(rad, arb_radref(c), e2 - e1); - } - - if (fmpz_cmp(arb_poly_radref(z), rad) < 0) - fmpz_set(arb_poly_radref(z), rad); - - fmpz_clear(rad); -} - -void -arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n) +arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n, int exact_const) { arb_poly_t t, u; arb_t c; @@ -72,18 +36,24 @@ arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n) arb_init(c, arb_poly_prec(z)); /* first term */ - fmpz_set(arb_midref(c), x->coeffs); - fmpz_set(arb_radref(c), arb_poly_radref(x)); - fmpz_set(arb_expref(c), arb_poly_expref(x)); - - arb_log(c, c); + if (!exact_const) + { + fmpz_set(arb_midref(c), x->coeffs); + fmpz_set(arb_radref(c), arb_poly_radref(x)); + fmpz_set(arb_expref(c), arb_poly_expref(x)); + arb_log(c, c); + } arb_poly_inv_series(t, x, n); + arb_poly_derivative(u, x); arb_poly_mullow(t, t, u, n - 1); arb_poly_integral(z, t); - _arb_poly_set_coeff_same_exp(z, c); + if (!exact_const) + { + _arb_poly_set_coeff_same_exp(z, 0, c); + } arb_clear(c); arb_poly_clear(t); diff --git a/arb_poly/set_coeff_same_exp.c b/arb_poly/set_coeff_same_exp.c new file mode 100644 index 00000000..6e717a87 --- /dev/null +++ b/arb_poly/set_coeff_same_exp.c @@ -0,0 +1,63 @@ +/*============================================================================= + + This file is part of FLINT. + + FLINT 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. + + FLINT 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 FLINT; 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 "arb_poly.h" + +void +_arb_poly_set_coeff_same_exp(arb_poly_t z, long i, const arb_t c) +{ + long e1, e2, j; + fmpz_t rad; + + e1 = z->exp; + e2 = c->exp; + + if (z->length < i + 1) + { + _arb_poly_fit_length(z, i + 1); + for (j = z->length; j < i; j++) + fmpz_zero(z->coeffs + j); + z->length = i + 1; + } + + fmpz_init(rad); + + if (e1 > e2) + { + fmpz_tdiv_q_2exp(z->coeffs + i, arb_midref(c), e1 - e2); + fmpz_cdiv_q_2exp(rad, arb_radref(c), e1 - e2); + fmpz_add_ui(rad, rad, 1); + } + else + { + fmpz_mul_2exp(z->coeffs + i, arb_midref(c), e2 - e1); + fmpz_mul_2exp(rad, arb_radref(c), e2 - e1); + } + + if (fmpz_cmp(arb_poly_radref(z), rad) < 0) + fmpz_set(arb_poly_radref(z), rad); + + fmpz_clear(rad); +}