From e7f0fff4900c00b8bdfbab28a0ddfe15b6151fb8 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 5 Jun 2014 00:43:30 +0200 Subject: [PATCH] port gamma function power series methods --- acb.h | 13 ++ acb_poly.h | 7 +- acb_poly/gamma_series.c | 305 +++++++++++++++++++++++++ acb_poly/lgamma_series.c | 124 ++++++++++ acb_poly/rgamma_series.c | 156 +++++++++++++ acb_poly/rising_ui_series.c | 120 ++++++++++ acb_poly/test/t-gamma_series.c | 114 ++++++++++ acb_poly/test/t-lgamma_series.c | 115 ++++++++++ acb_poly/test/t-rgamma_series.c | 113 ++++++++++ acb_poly/test/t-rising_ui_series.c | 128 +++++++++++ arb.h | 9 +- arb_poly.h | 10 +- arb_poly/gamma_series.c | 350 +++++++++++++++++++++++++++++ arb_poly/lgamma_series.c | 137 +++++++++++ arb_poly/rgamma_series.c | 190 ++++++++++++++++ arb_poly/rising_ui_series.c | 119 ++++++++++ arb_poly/test/t-gamma_series.c | 121 ++++++++++ arb_poly/test/t-lgamma_series.c | 123 ++++++++++ arb_poly/test/t-rgamma_series.c | 120 ++++++++++ arb_poly/test/t-rising_ui_series.c | 128 +++++++++++ doc/source/acb.rst | 8 +- doc/source/acb_poly.rst | 29 +++ doc/source/arb_poly.rst | 32 +++ 23 files changed, 2556 insertions(+), 15 deletions(-) create mode 100644 acb_poly/gamma_series.c create mode 100644 acb_poly/lgamma_series.c create mode 100644 acb_poly/rgamma_series.c create mode 100644 acb_poly/rising_ui_series.c create mode 100644 acb_poly/test/t-gamma_series.c create mode 100644 acb_poly/test/t-lgamma_series.c create mode 100644 acb_poly/test/t-rgamma_series.c create mode 100644 acb_poly/test/t-rising_ui_series.c create mode 100644 arb_poly/gamma_series.c create mode 100644 arb_poly/lgamma_series.c create mode 100644 arb_poly/rgamma_series.c create mode 100644 arb_poly/rising_ui_series.c create mode 100644 arb_poly/test/t-gamma_series.c create mode 100644 arb_poly/test/t-lgamma_series.c create mode 100644 arb_poly/test/t-rgamma_series.c create mode 100644 arb_poly/test/t-rising_ui_series.c diff --git a/acb.h b/acb.h index 52c9a7fe..cd5bac0c 100644 --- a/acb.h +++ b/acb.h @@ -859,6 +859,19 @@ _acb_vec_add_error_arf_vec(acb_ptr res, arf_srcptr err, long len) acb_add_error_arf(res + i, err + i); } +static __inline__ void +_acb_vec_add_error_mag_vec(acb_ptr res, mag_srcptr err, long len) +{ + long i; + for (i = 0; i < len; i++) + { + mag_add(arb_radref(acb_realref(res + i)), + arb_radref(acb_realref(res + i)), err + i); + mag_add(arb_radref(acb_imagref(res + i)), + arb_radref(acb_imagref(res + i)), err + i); + } +} + static __inline__ void _acb_vec_indeterminate(acb_ptr vec, long len) { diff --git a/acb_poly.h b/acb_poly.h index 86ecb8c0..c8099c92 100644 --- a/acb_poly.h +++ b/acb_poly.h @@ -490,9 +490,6 @@ void _acb_poly_tan_series(acb_ptr g, acb_srcptr h, long hlen, long len, long pre void acb_poly_tan_series(acb_poly_t g, const acb_poly_t h, long n, long prec); -/* -TBD: - void _acb_poly_gamma_series(acb_ptr res, acb_srcptr h, long hlen, long len, long prec); void acb_poly_gamma_series(acb_poly_t res, const acb_poly_t f, long n, long prec); @@ -509,6 +506,10 @@ void _acb_poly_rising_ui_series(acb_ptr res, acb_srcptr f, long flen, ulong r, l void acb_poly_rising_ui_series(acb_poly_t res, const acb_poly_t f, ulong r, long trunc, long prec); + +/* +TBD: + void _acb_poly_zeta_series(acb_ptr res, acb_srcptr h, long hlen, const acb_t a, int deflate, long len, long prec); void acb_poly_zeta_series(acb_poly_t res, const acb_poly_t f, const acb_t a, int deflate, long n, long prec); diff --git a/acb_poly/gamma_series.c b/acb_poly/gamma_series.c new file mode 100644 index 00000000..ff8c63bb --- /dev/null +++ b/acb_poly/gamma_series.c @@ -0,0 +1,305 @@ +/*============================================================================= + + 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 "acb_poly.h" +#include "gamma.h" +#include "zeta.h" + +void acb_gamma_stirling_bound(mag_ptr err, const acb_t x, long k0, long knum, long n); + +void acb_gamma_stirling_choose_param(int * reflect, long * r, long * n, + const acb_t x, int use_reflect, int digamma, long prec); + +void arb_gamma_stirling_coeff(arb_t b, ulong k, int digamma, long prec); + +static void +bsplit(acb_ptr Q, acb_ptr T, const acb_t z, long a, long b, long num, long prec) +{ + if (b - a == 1) + { + arb_gamma_stirling_coeff(acb_realref(T), a, 0, prec); + arb_zero(acb_imagref(T)); + + if (a == 1) + { /* (z + t) */ + acb_set(Q, z); + if (num > 1) acb_one(Q + 1); + if (num > 2) acb_zero(Q + 2); + } + else + { /* (z + t)^2 */ + acb_mul(Q, z, z, prec); /* TODO: precompute */ + if (num > 1) acb_mul_2exp_si(Q + 1, z, 1); + if (num > 2) acb_one(Q + 2); + } + } + else + { + long m, n1, n2, q1len, q2len, t1len, t2len, qlen, tlen, alloc; + acb_ptr Q1, T1, Q2, T2; + + m = a + (b - a) / 2; + + n1 = m - a; + n2 = b - m; + q1len = FLINT_MIN(2 * n1 + 1, num); + t1len = FLINT_MIN(2 * n1 - 1, num); + q2len = FLINT_MIN(2 * n2 + 1, num); + t2len = FLINT_MIN(2 * n2 - 1, num); + qlen = FLINT_MIN(q1len + q2len - 1, num); + tlen = FLINT_MIN(t1len + q2len - 1, num); + + alloc = q1len + q2len + t1len + t2len; + Q1 = _acb_vec_init(alloc); + Q2 = Q1 + q1len; + T1 = Q2 + q2len; + T2 = T1 + t1len; + + bsplit(Q1, T1, z, a, m, num, prec); + bsplit(Q2, T2, z, m, b, num, prec); + + _acb_poly_mullow(Q, Q2, q2len, Q1, q1len, qlen, prec); + _acb_poly_mullow(T, Q2, q2len, T1, t1len, tlen, prec); + _acb_poly_add(T, T, tlen, T2, t2len, prec); + + _acb_vec_clear(Q1, alloc); + } +} + +void +_acb_poly_mullow_cpx(acb_ptr res, acb_srcptr src, long len, const acb_t c, long trunc, long prec) +{ + long i; + + if (len < trunc) + acb_set(res + len, src + len - 1); + + for (i = len - 1; i > 0; i--) + { + acb_mul(res + i, src + i, c, prec); + acb_add(res + i, res + i, src + i - 1, prec); + } + + acb_mul(res, src, c, prec); +} + +void +_acb_poly_log_cpx_series(acb_ptr res, const acb_t c, long num, long prec) +{ + long i; + + for (i = 0; i < num; i++) + { + if (i == 0) + acb_log(res + i, c, prec); + else if (i == 1) + acb_inv(res + i, c, prec); + else + acb_mul(res + i, res + i - 1, res + 1, prec); + } + + for (i = 2; i < num; i++) + { + acb_div_ui(res + i, res + i, i, prec); + + if (i % 2 == 0) + acb_neg(res + i, res + i); + } +} + +void +_acb_poly_gamma_stirling_eval(acb_ptr res, const acb_t z, long n, long num, long prec) +{ + long tlen, qlen; + acb_ptr T, Q; + mag_ptr err; + acb_t c; + + T = _acb_vec_init(num); + Q = _acb_vec_init(num); + err = _mag_vec_init(num); + acb_init(c); + + acb_gamma_stirling_bound(err, z, 0, num, n); + + if (n <= 1) + { + _acb_vec_zero(res, num); + } + else + { + qlen = FLINT_MIN(2 * (n - 1) + 1, num); + tlen = FLINT_MIN(2 * (n - 1) - 1, num); + bsplit(Q, T, z, 1, n, num, prec); + _acb_poly_div_series(res, T, tlen, Q, qlen, num, prec); + } + + /* ((z-1/2) + t) * log(z+t) */ + _acb_poly_log_cpx_series(T, z, num, prec); + acb_one(c); + acb_mul_2exp_si(c, c, -1); + acb_sub(c, z, c, prec); + _acb_poly_mullow_cpx(T, T, num, c, num, prec); + + /* constant term */ + arb_const_log_sqrt2pi(acb_realref(c), prec); + arb_zero(acb_imagref(c)); + acb_add(T, T, c, prec); + + /* subtract (z+t) */ + acb_sub(T, T, z, prec); + if (num > 1) + acb_sub_ui(T + 1, T + 1, 1, prec); + + _acb_vec_add(res, res, T, num, prec); + + _acb_vec_add_error_mag_vec(res, err, num); + + _acb_vec_clear(T, num); + _acb_vec_clear(Q, num); + _mag_vec_clear(err, num); + acb_clear(c); +} + +void +_acb_poly_gamma_series(acb_ptr res, acb_srcptr h, long hlen, long len, long prec) +{ + int reflect; + long i, rflen, r, n, wp; + acb_ptr t, u, v; + acb_struct f[2]; + + hlen = FLINT_MIN(hlen, len); + wp = prec + FLINT_BIT_COUNT(prec); + + t = _acb_vec_init(len); + u = _acb_vec_init(len); + v = _acb_vec_init(len); + acb_init(f); + acb_init(f + 1); + + /* TODO: use real code at real numbers */ + if (0) + { + } + else + { + /* otherwise use Stirling series */ + acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); + + /* gamma(h) = (rf(1-h, r) * pi) / (gamma(1-h+r) sin(pi h)), h = h0 + t*/ + if (reflect) + { + /* u = 1/gamma(r+1-h) */ + acb_sub_ui(f, h, r + 1, wp); + acb_neg(f, f); + _acb_poly_gamma_stirling_eval(t, f, n, len, wp); + _acb_vec_neg(t, t, len); + _acb_poly_exp_series(u, t, len, len, wp); + for (i = 1; i < len; i += 2) + acb_neg(u + i, u + i); + + /* v = 1/sin(pi x) */ + acb_const_pi(f + 1, wp); + acb_mul(f, h, f + 1, wp); + _acb_poly_sin_series(t, f, 2, len, wp); + _acb_poly_inv_series(v, t, len, len, wp); + + _acb_poly_mullow(t, u, len, v, len, len, wp); + + /* rf(1-h,r) * pi */ + if (r == 0) + { + rflen = 1; + acb_const_pi(u, wp); + } + else + { + acb_sub_ui(f, h, 1, wp); + acb_neg(f, f); + acb_set_si(f + 1, -1); + rflen = FLINT_MIN(len, r + 1); + _acb_poly_rising_ui_series(u, f, FLINT_MIN(2, len), r, rflen, wp); + acb_const_pi(v, wp); + _acb_vec_scalar_mul(u, u, rflen, v, wp); + } + + /* multiply by rising factorial */ + _acb_poly_mullow(v, t, len, u, rflen, len, wp); + } + else + { + /* gamma(h) = gamma(h+r) / rf(h,r) */ + if (r == 0) + { + acb_add_ui(f, h, r, wp); + _acb_poly_gamma_stirling_eval(t, f, n, len, wp); + _acb_poly_exp_series(v, t, len, len, wp); + } + else + { + /* TODO: div_series may be better (once it has a good basecase), + if the rising factorial is short */ + acb_set(f, h); + acb_one(f + 1); + rflen = FLINT_MIN(len, r + 1); + _acb_poly_rising_ui_series(u, f, FLINT_MIN(2, len), r, rflen, wp); + _acb_poly_inv_series(t, u, rflen, len, wp); + + acb_add_ui(f, h, r, wp); + _acb_poly_gamma_stirling_eval(v, f, n, len, wp); + _acb_poly_exp_series(u, v, len, len, wp); + + _acb_poly_mullow(v, u, len, t, len, len, wp); + } + } + } + + /* compose with nonconstant part */ + acb_zero(t); + _acb_vec_set(t + 1, h + 1, hlen - 1); + _acb_poly_compose_series(res, v, len, t, hlen, len, prec); + + acb_clear(f); + acb_clear(f + 1); + _acb_vec_clear(t, len); + _acb_vec_clear(u, len); + _acb_vec_clear(v, len); +} + +void +acb_poly_gamma_series(acb_poly_t res, const acb_poly_t f, long n, long prec) +{ + acb_poly_fit_length(res, n); + + if (f->length == 0 || n == 0) + _acb_vec_indeterminate(res->coeffs, n); + else + _acb_poly_gamma_series(res->coeffs, f->coeffs, f->length, n, prec); + + _acb_poly_set_length(res, n); + _acb_poly_normalise(res); +} + diff --git a/acb_poly/lgamma_series.c b/acb_poly/lgamma_series.c new file mode 100644 index 00000000..06c18e2f --- /dev/null +++ b/acb_poly/lgamma_series.c @@ -0,0 +1,124 @@ +/*============================================================================= + + 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 "acb_poly.h" +#include "gamma.h" +#include "zeta.h" + +void +_acb_log_rising_correct_branch(acb_t t, + const acb_t t_wrong, const acb_t z, ulong r, long prec); + +void acb_gamma_stirling_choose_param(int * reflect, long * r, long * n, + const acb_t x, int use_reflect, int digamma, long prec); + +void +_acb_poly_gamma_stirling_eval(acb_ptr res, const acb_t z, long n, long num, long prec); + +static __inline__ void +_log_rising_ui_series(acb_ptr t, const acb_t x, long r, long len, long prec) +{ + acb_struct f[2]; + long rflen; + + acb_init(f); + acb_init(f + 1); + + acb_set(f, x); + acb_one(f + 1); + + rflen = FLINT_MIN(len, r + 1); + _acb_poly_rising_ui_series(t, f, FLINT_MIN(2, len), r, rflen, prec); + _acb_poly_log_series(t, t, rflen, len, prec); + + _acb_log_rising_correct_branch(t, t, x, r, prec); + + acb_clear(f); + acb_clear(f + 1); +} + +void +_acb_poly_lgamma_series(acb_ptr res, acb_srcptr h, long hlen, long len, long prec) +{ + int reflect; + long r, n, wp; + acb_t zr; + acb_ptr t, u; + + hlen = FLINT_MIN(hlen, len); + wp = prec + FLINT_BIT_COUNT(prec); + + t = _acb_vec_init(len); + u = _acb_vec_init(len); + acb_init(zr); + + /* TODO: use real code at real numbers */ + if (0) + { + } + else if (len <= 2) + { + acb_lgamma(u, h, wp); + if (len == 2) + acb_digamma(u + 1, h, wp); + } + else + { + /* otherwise use Stirling series */ + acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 0, 0, wp); + acb_add_ui(zr, h, r, wp); + _acb_poly_gamma_stirling_eval(u, zr, n, len, wp); + + if (r != 0) + { + _log_rising_ui_series(t, h, r, len, wp); + _acb_vec_sub(u, u, t, len, wp); + } + } + + /* compose with nonconstant part */ + acb_zero(t); + _acb_vec_set(t + 1, h + 1, hlen - 1); + _acb_poly_compose_series(res, u, len, t, hlen, len, prec); + + acb_clear(zr); + _acb_vec_clear(t, len); + _acb_vec_clear(u, len); +} + +void +acb_poly_lgamma_series(acb_poly_t res, const acb_poly_t f, long n, long prec) +{ + acb_poly_fit_length(res, n); + + if (f->length == 0 || n == 0) + _acb_vec_indeterminate(res->coeffs, n); + else + _acb_poly_lgamma_series(res->coeffs, f->coeffs, f->length, n, prec); + + _acb_poly_set_length(res, n); + _acb_poly_normalise(res); +} + diff --git a/acb_poly/rgamma_series.c b/acb_poly/rgamma_series.c new file mode 100644 index 00000000..09f0afc0 --- /dev/null +++ b/acb_poly/rgamma_series.c @@ -0,0 +1,156 @@ +/*============================================================================= + + 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 "acb_poly.h" +#include "gamma.h" +#include "zeta.h" + +void acb_gamma_stirling_choose_param(int * reflect, long * r, long * n, + const acb_t x, int use_reflect, int digamma, long prec); + +void +_acb_poly_gamma_stirling_eval(acb_ptr res, const acb_t z, long n, long num, long prec); + +void +_acb_poly_rgamma_series(acb_ptr res, acb_srcptr h, long hlen, long len, long prec) +{ + int reflect; + long i, rflen, r, n, wp; + acb_ptr t, u, v; + acb_struct f[2]; + + hlen = FLINT_MIN(hlen, len); + wp = prec + FLINT_BIT_COUNT(prec); + + t = _acb_vec_init(len); + u = _acb_vec_init(len); + v = _acb_vec_init(len); + acb_init(f); + acb_init(f + 1); + + /* TODO: use real code at real numbers */ + if (0) + { + } + else + { + /* otherwise use Stirling series */ + acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); + + /* rgamma(h) = (gamma(1-h+r) sin(pi h)) / (rf(1-h, r) * pi), h = h0 + t*/ + if (reflect) + { + /* u = gamma(r+1-h) */ + acb_sub_ui(f, h, r + 1, wp); + acb_neg(f, f); + _acb_poly_gamma_stirling_eval(t, f, n, len, wp); + _acb_poly_exp_series(u, t, len, len, wp); + for (i = 1; i < len; i += 2) + acb_neg(u + i, u + i); + + /* v = sin(pi x) */ + acb_const_pi(f + 1, wp); + acb_mul(f, h, f + 1, wp); + _acb_poly_sin_series(v, f, 2, len, wp); + + _acb_poly_mullow(t, u, len, v, len, len, wp); + + /* rf(1-h,r) * pi */ + if (r == 0) + { + acb_const_pi(u, wp); + _acb_vec_scalar_div(v, t, len, u, wp); + } + else + { + acb_sub_ui(f, h, 1, wp); + acb_neg(f, f); + acb_set_si(f + 1, -1); + rflen = FLINT_MIN(len, r + 1); + _acb_poly_rising_ui_series(v, f, FLINT_MIN(2, len), r, rflen, wp); + acb_const_pi(u, wp); + _acb_vec_scalar_mul(v, v, rflen, u, wp); + + /* divide by rising factorial */ + /* TODO: might better to use div_series, when it has a good basecase */ + _acb_poly_inv_series(u, v, rflen, len, wp); + _acb_poly_mullow(v, t, len, u, len, len, wp); + } + } + else + { + /* rgamma(h) = rgamma(h+r) rf(h,r) */ + if (r == 0) + { + acb_add_ui(f, h, r, wp); + _acb_poly_gamma_stirling_eval(t, f, n, len, wp); + _acb_vec_neg(t, t, len); + _acb_poly_exp_series(v, t, len, len, wp); + } + else + { + acb_set(f, h); + acb_one(f + 1); + rflen = FLINT_MIN(len, r + 1); + _acb_poly_rising_ui_series(t, f, FLINT_MIN(2, len), r, rflen, wp); + + acb_add_ui(f, h, r, wp); + _acb_poly_gamma_stirling_eval(v, f, n, len, wp); + _acb_vec_neg(v, v, len); + _acb_poly_exp_series(u, v, len, len, wp); + + _acb_poly_mullow(v, u, len, t, rflen, len, wp); + } + } + } + + /* compose with nonconstant part */ + acb_zero(t); + _acb_vec_set(t + 1, h + 1, hlen - 1); + _acb_poly_compose_series(res, v, len, t, hlen, len, prec); + + acb_clear(f); + acb_clear(f + 1); + _acb_vec_clear(t, len); + _acb_vec_clear(u, len); + _acb_vec_clear(v, len); +} + +void +acb_poly_rgamma_series(acb_poly_t res, const acb_poly_t f, long n, long prec) +{ + if (f->length == 0 || n == 0) + { + acb_poly_zero(res); + } + else + { + acb_poly_fit_length(res, n); + _acb_poly_rgamma_series(res->coeffs, f->coeffs, f->length, n, prec); + _acb_poly_set_length(res, n); + _acb_poly_normalise(res); + } +} + diff --git a/acb_poly/rising_ui_series.c b/acb_poly/rising_ui_series.c new file mode 100644 index 00000000..f748765b --- /dev/null +++ b/acb_poly/rising_ui_series.c @@ -0,0 +1,120 @@ +/*============================================================================= + + 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, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "acb_poly.h" +#include "gamma.h" + +static void +_acb_poly_rising_ui_series_bsplit(acb_ptr res, + acb_srcptr f, long flen, ulong a, ulong b, + long trunc, long prec) +{ + flen = FLINT_MIN(flen, trunc); + + if (b - a == 1) + { + acb_add_ui(res, f, a, prec); + _acb_vec_set(res + 1, f + 1, flen - 1); + } + else + { + acb_ptr L, R; + long len1, len2; + + long m = a + (b - a) / 2; + + len1 = poly_pow_length(flen, m - a, trunc); + len2 = poly_pow_length(flen, b - m, trunc); + + L = _acb_vec_init(len1 + len2); + R = L + len1; + + _acb_poly_rising_ui_series_bsplit(L, f, flen, a, m, trunc, prec); + _acb_poly_rising_ui_series_bsplit(R, f, flen, m, b, trunc, prec); + + _acb_poly_mullow(res, L, len1, R, len2, + FLINT_MIN(trunc, len1 + len2 - 1), prec); + + _acb_vec_clear(L, len1 + len2); + } +} + +void +_acb_poly_rising_ui_series(acb_ptr res, + acb_srcptr f, long flen, ulong r, + long trunc, long prec) +{ + if (trunc == 1 || flen == 1) + { + acb_rising_ui(res, f, r, prec); + _acb_vec_zero(res + 1, trunc - 1); + } + else if (trunc == 2) + { + acb_rising2_ui(res, res + 1, f, r, prec); + acb_mul(res + 1, res + 1, f + 1, prec); + } + else + { + _acb_poly_rising_ui_series_bsplit(res, f, flen, 0, r, trunc, prec); + } +} + +void +acb_poly_rising_ui_series(acb_poly_t res, const acb_poly_t f, ulong r, long trunc, long prec) +{ + long len; + + if (f->length == 0 && r != 0) + { + acb_poly_zero(res); + return; + } + + if (r == 0) + { + acb_poly_one(res); + return; + } + + len = poly_pow_length(f->length, r, trunc); + + if (f == res) + { + acb_poly_t tmp; + acb_poly_init(tmp); + acb_poly_rising_ui_series(tmp, f, r, len, prec); + acb_poly_swap(tmp, res); + acb_poly_clear(tmp); + } + else + { + acb_poly_fit_length(res, len); + _acb_poly_rising_ui_series(res->coeffs, f->coeffs, f->length, r, len, prec); + _acb_poly_set_length(res, len); + _acb_poly_normalise(res); + } +} + diff --git a/acb_poly/test/t-gamma_series.c b/acb_poly/test/t-gamma_series.c new file mode 100644 index 00000000..65b13b4a --- /dev/null +++ b/acb_poly/test/t-gamma_series.c @@ -0,0 +1,114 @@ +/*============================================================================= + + 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, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "acb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("gamma_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long m, n1, n2, rbits1, rbits2, rbits3; + acb_poly_t a, b, c, d; + + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + rbits3 = 2 + n_randint(state, 200); + + m = 1 + n_randint(state, 25); + n1 = 1 + n_randint(state, 25); + n2 = 1 + n_randint(state, 25); + + acb_poly_init(a); + acb_poly_init(b); + acb_poly_init(c); + acb_poly_init(d); + + acb_poly_randtest(a, state, m, rbits1, 3); + + acb_poly_gamma_series(b, a, n1, rbits2); + acb_poly_gamma_series(c, a, n2, rbits3); + + acb_poly_set(d, b); + acb_poly_truncate(d, FLINT_MIN(n1, n2)); + acb_poly_truncate(c, FLINT_MIN(n1, n2)); + + if (!acb_poly_overlaps(c, d)) + { + printf("FAIL\n\n"); + printf("n1 = %ld, n2 = %ld, bits2 = %ld, bits3 = %ld\n", n1, n2, rbits2, rbits3); + + printf("a = "); acb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); acb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); acb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + /* check gamma(a) * a = gamma(a+1) */ + acb_poly_mullow(c, b, a, n1, rbits2); + + acb_poly_set(d, a); + acb_add_ui(d->coeffs, d->coeffs, 1, rbits2); + acb_poly_gamma_series(d, d, n1, rbits2); + + if (!acb_poly_overlaps(c, d)) + { + printf("FAIL (functional equation, n1 = %ld)\n\n", n1); + + printf("a = "); acb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); acb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); acb_poly_printd(c, 15); printf("\n\n"); + printf("d = "); acb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + acb_poly_gamma_series(a, a, n1, rbits2); + if (!acb_poly_overlaps(a, b)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + acb_poly_clear(a); + acb_poly_clear(b); + acb_poly_clear(c); + acb_poly_clear(d); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/acb_poly/test/t-lgamma_series.c b/acb_poly/test/t-lgamma_series.c new file mode 100644 index 00000000..ccbb11d9 --- /dev/null +++ b/acb_poly/test/t-lgamma_series.c @@ -0,0 +1,115 @@ +/*============================================================================= + + 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, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "acb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("lgamma_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 500; iter++) + { + long m, n1, n2, rbits1, rbits2, rbits3; + acb_poly_t a, b, c, d; + + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + rbits3 = 2 + n_randint(state, 200); + + m = 1 + n_randint(state, 30); + n1 = 1 + n_randint(state, 30); + n2 = 1 + n_randint(state, 30); + + acb_poly_init(a); + acb_poly_init(b); + acb_poly_init(c); + acb_poly_init(d); + + acb_poly_randtest(a, state, m, rbits1, 3); + + acb_poly_lgamma_series(b, a, n1, rbits2); + acb_poly_lgamma_series(c, a, n2, rbits3); + + acb_poly_set(d, b); + acb_poly_truncate(d, FLINT_MIN(n1, n2)); + acb_poly_truncate(c, FLINT_MIN(n1, n2)); + + if (!acb_poly_overlaps(c, d)) + { + printf("FAIL\n\n"); + printf("n1 = %ld, n2 = %ld, bits2 = %ld, bits3 = %ld\n", n1, n2, rbits2, rbits3); + + printf("a = "); acb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); acb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); acb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + /* check loggamma(a) + log(a) = loggamma(a+1) */ + acb_poly_log_series(c, a, n1, rbits2); + acb_poly_add(c, b, c, rbits2); + + acb_poly_set(d, a); + acb_add_ui(d->coeffs, d->coeffs, 1, rbits2); + acb_poly_lgamma_series(d, d, n1, rbits2); + + if (!acb_poly_overlaps(c, d)) + { + printf("FAIL (functional equation)\n\n"); + + printf("a = "); acb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); acb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); acb_poly_printd(c, 15); printf("\n\n"); + printf("d = "); acb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + acb_poly_lgamma_series(a, a, n1, rbits2); + if (!acb_poly_overlaps(a, b)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + acb_poly_clear(a); + acb_poly_clear(b); + acb_poly_clear(c); + acb_poly_clear(d); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/acb_poly/test/t-rgamma_series.c b/acb_poly/test/t-rgamma_series.c new file mode 100644 index 00000000..a2d9ac16 --- /dev/null +++ b/acb_poly/test/t-rgamma_series.c @@ -0,0 +1,113 @@ +/*============================================================================= + + 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, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "acb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("rgamma_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long m, n1, n2, rbits1, rbits2, rbits3; + acb_poly_t a, b, c, d; + + rbits1 = 2 + n_randint(state, 200); + rbits2 = 2 + n_randint(state, 200); + rbits3 = 2 + n_randint(state, 200); + + m = 1 + n_randint(state, 30); + n1 = 1 + n_randint(state, 30); + n2 = 1 + n_randint(state, 30); + + acb_poly_init(a); + acb_poly_init(b); + acb_poly_init(c); + acb_poly_init(d); + + acb_poly_randtest(a, state, m, rbits1, 3); + + acb_poly_rgamma_series(b, a, n1, rbits2); + acb_poly_rgamma_series(c, a, n2, rbits3); + + acb_poly_set(d, b); + acb_poly_truncate(d, FLINT_MIN(n1, n2)); + acb_poly_truncate(c, FLINT_MIN(n1, n2)); + + if (!acb_poly_overlaps(c, d)) + { + printf("FAIL\n\n"); + printf("n1 = %ld, n2 = %ld, bits2 = %ld, bits3 = %ld\n", n1, n2, rbits2, rbits3); + + printf("a = "); acb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); acb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); acb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + /* check rgamma(a) = a gamma(a+1) */ + acb_poly_set(d, a); + acb_add_ui(d->coeffs, d->coeffs, 1, rbits2); + acb_poly_rgamma_series(d, d, n1, rbits2); + acb_poly_mullow(c, d, a, n1, rbits2); + + if (!acb_poly_overlaps(b, c)) + { + printf("FAIL (functional equation, n1 = %ld)\n\n", n1); + + printf("a = "); acb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); acb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); acb_poly_printd(c, 15); printf("\n\n"); + printf("d = "); acb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + acb_poly_rgamma_series(a, a, n1, rbits2); + if (!acb_poly_overlaps(a, b)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + acb_poly_clear(a); + acb_poly_clear(b); + acb_poly_clear(c); + acb_poly_clear(d); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/acb_poly/test/t-rising_ui_series.c b/acb_poly/test/t-rising_ui_series.c new file mode 100644 index 00000000..00bb905f --- /dev/null +++ b/acb_poly/test/t-rising_ui_series.c @@ -0,0 +1,128 @@ +/*============================================================================= + + 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 "acb_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("rising_ui_series...."); + fflush(stdout); + + flint_randinit(state); + + /* check rf(f, a) * rf(f + a, b) = rf(f, a + b) */ + for (iter = 0; iter < 1000; iter++) + { + long bits, trunc; + ulong a, b; + acb_poly_t f, g, h1, h2, h1h2, h3; + + bits = 2 + n_randint(state, 200); + trunc = 1 + n_randint(state, 20); + a = n_randint(state, 10); + b = n_randint(state, 10); + + acb_poly_init(f); + acb_poly_init(g); + acb_poly_init(h1); + acb_poly_init(h2); + acb_poly_init(h1h2); + acb_poly_init(h3); + + acb_poly_randtest(f, state, 1 + n_randint(state, 20), bits, 4); + acb_poly_set(g, f); + + /* g = f + 1 */ + if (g->length == 0) + { + acb_poly_fit_length(g, 1); + acb_set_ui(g->coeffs, a); + _acb_poly_set_length(g, 1); + _acb_poly_normalise(g); + } + else + { + acb_add_ui(g->coeffs, g->coeffs, a, bits); + _acb_poly_normalise(g); + } + + acb_poly_rising_ui_series(h1, f, a, trunc, bits); + acb_poly_rising_ui_series(h2, g, b, trunc, bits); + acb_poly_rising_ui_series(h3, f, a + b, trunc, bits); + + acb_poly_mullow(h1h2, h1, h2, trunc, bits); + + if (!acb_poly_overlaps(h1h2, h3)) + { + printf("FAIL\n\n"); + printf("bits = %ld\n", bits); + printf("trunc = %ld\n", trunc); + printf("a = %lu\n", a); + printf("b = %lu\n", a); + + printf("f = "); acb_poly_printd(f, 15); printf("\n\n"); + printf("g = "); acb_poly_printd(g, 15); printf("\n\n"); + printf("h1 = "); acb_poly_printd(h1, 15); printf("\n\n"); + printf("h2 = "); acb_poly_printd(h2, 15); printf("\n\n"); + printf("h1h2 = "); acb_poly_printd(h1h2, 15); printf("\n\n"); + printf("h3 = "); acb_poly_printd(h3, 15); printf("\n\n"); + + abort(); + } + + acb_poly_rising_ui_series(f, f, a, trunc, bits); + + if (!acb_poly_equal(f, h1)) + { + printf("FAIL (aliasing)\n\n"); + + printf("bits = %ld\n", bits); + printf("trunc = %ld\n", trunc); + printf("a = %lu\n", a); + printf("b = %lu\n", a); + + printf("f = "); acb_poly_printd(f, 15); printf("\n\n"); + printf("h1 = "); acb_poly_printd(h1, 15); printf("\n\n"); + + abort(); + } + + acb_poly_clear(f); + acb_poly_clear(g); + acb_poly_clear(h1); + acb_poly_clear(h2); + acb_poly_clear(h1h2); + acb_poly_clear(h3); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arb.h b/arb.h index 6b4e2ca9..210ce5c5 100644 --- a/arb.h +++ b/arb.h @@ -854,7 +854,6 @@ _arb_vec_set_powers(arb_ptr xs, const arb_t x, long len, long prec) } } -/* TODO: mag version ? */ static __inline__ void _arb_vec_add_error_arf_vec(arb_ptr res, arf_srcptr err, long len) { @@ -863,6 +862,14 @@ _arb_vec_add_error_arf_vec(arb_ptr res, arf_srcptr err, long len) arb_add_error_arf(res + i, err + i); } +static __inline__ void +_arb_vec_add_error_mag_vec(arb_ptr res, mag_srcptr err, long len) +{ + long i; + for (i = 0; i < len; i++) + mag_add(arb_radref(res + i), arb_radref(res + i), err + i); +} + static __inline__ void _arb_vec_indeterminate(arb_ptr vec, long len) { diff --git a/arb_poly.h b/arb_poly.h index 0669d84f..58c7d24d 100644 --- a/arb_poly.h +++ b/arb_poly.h @@ -546,25 +546,21 @@ void arb_poly_evaluate2_acb_rectangular(acb_t y, acb_t z, const arb_poly_t f, co void _arb_poly_evaluate2_acb(acb_t y, acb_t z, arb_srcptr f, long len, const acb_t x, long prec); void arb_poly_evaluate2_acb(acb_t y, acb_t z, const arb_poly_t f, const acb_t x, long prec); -/* -TBD: - void _arb_poly_gamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long prec); - void arb_poly_gamma_series(arb_poly_t res, const arb_poly_t f, long n, long prec); void _arb_poly_rgamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long prec); - void arb_poly_rgamma_series(arb_poly_t res, const arb_poly_t f, long n, long prec); void _arb_poly_lgamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long prec); - void arb_poly_lgamma_series(arb_poly_t res, const arb_poly_t f, long n, long prec); void _arb_poly_rising_ui_series(arb_ptr res, arb_srcptr f, long flen, ulong r, long trunc, long prec); - void arb_poly_rising_ui_series(arb_poly_t res, const arb_poly_t f, ulong r, long trunc, long prec); +/* +TBD: + void _arb_poly_zeta_series(arb_ptr res, arb_srcptr h, long hlen, const arb_t a, int deflate, long len, long prec); void arb_poly_zeta_series(arb_poly_t res, const arb_poly_t f, const arb_t a, int deflate, long n, long prec); diff --git a/arb_poly/gamma_series.c b/arb_poly/gamma_series.c new file mode 100644 index 00000000..ebd5db2e --- /dev/null +++ b/arb_poly/gamma_series.c @@ -0,0 +1,350 @@ +/*============================================================================= + + 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 "arb_poly.h" +#include "gamma.h" + +void arb_gamma_stirling_bound(mag_ptr err, const arb_t x, long k0, long knum, long n); + +void arb_gamma_stirling_choose_param(int * reflect, long * r, long * n, + const arb_t x, int use_reflect, int digamma, long prec); + +void arb_gamma_stirling_coeff(arb_t b, ulong k, int digamma, long prec); + +long +arf_get_si(const arf_t x, arf_rnd_t rnd) +{ + fmpr_t t; + long v; + fmpr_init(t); + arf_get_fmpr(t, x); + v = fmpr_get_si(t, rnd); + fmpr_clear(t); + return v; +} + +void +_arb_poly_lgamma_series_at_one(arb_ptr u, long len, long prec) +{ + fmprb_ptr t; + long i; + + t = _fmprb_vec_init(len); + gamma_lgamma_series_at_one(t, len, prec); + for (i = 0; i < len; i++) + arb_set_fmprb(u + i, t + i); + + _fmprb_vec_clear(t, len); +} + + +static void +bsplit(arb_ptr Q, arb_ptr T, const arb_t z, long a, long b, long num, long prec) +{ + if (b - a == 1) + { + arb_gamma_stirling_coeff(T, a, 0, prec); + + if (a == 1) + { /* (z + t) */ + arb_set(Q, z); + if (num > 1) arb_one(Q + 1); + if (num > 2) arb_zero(Q + 2); + } + else + { /* (z + t)^2 */ + arb_mul(Q, z, z, prec); /* TODO: precompute */ + if (num > 1) arb_mul_2exp_si(Q + 1, z, 1); + if (num > 2) arb_one(Q + 2); + } + } + else + { + long m, n1, n2, q1len, q2len, t1len, t2len, qlen, tlen, alloc; + arb_ptr Q1, T1, Q2, T2; + + m = a + (b - a) / 2; + + n1 = m - a; + n2 = b - m; + q1len = FLINT_MIN(2 * n1 + 1, num); + t1len = FLINT_MIN(2 * n1 - 1, num); + q2len = FLINT_MIN(2 * n2 + 1, num); + t2len = FLINT_MIN(2 * n2 - 1, num); + qlen = FLINT_MIN(q1len + q2len - 1, num); + tlen = FLINT_MIN(t1len + q2len - 1, num); + + alloc = q1len + q2len + t1len + t2len; + Q1 = _arb_vec_init(alloc); + Q2 = Q1 + q1len; + T1 = Q2 + q2len; + T2 = T1 + t1len; + + bsplit(Q1, T1, z, a, m, num, prec); + bsplit(Q2, T2, z, m, b, num, prec); + + _arb_poly_mullow(Q, Q2, q2len, Q1, q1len, qlen, prec); + _arb_poly_mullow(T, Q2, q2len, T1, t1len, tlen, prec); + _arb_poly_add(T, T, tlen, T2, t2len, prec); + + _arb_vec_clear(Q1, alloc); + } +} + +void +_arb_poly_mullow_cpx(arb_ptr res, arb_srcptr src, long len, const arb_t c, long trunc, long prec) +{ + long i; + + if (len < trunc) + arb_set(res + len, src + len - 1); + + for (i = len - 1; i > 0; i--) + { + arb_mul(res + i, src + i, c, prec); + arb_add(res + i, res + i, src + i - 1, prec); + } + + arb_mul(res, src, c, prec); +} + +void +_arb_poly_log_cpx_series(arb_ptr res, const arb_t c, long num, long prec) +{ + long i; + + for (i = 0; i < num; i++) + { + if (i == 0) + arb_log(res + i, c, prec); + else if (i == 1) + arb_inv(res + i, c, prec); + else + arb_mul(res + i, res + i - 1, res + 1, prec); + } + + for (i = 2; i < num; i++) + { + arb_div_ui(res + i, res + i, i, prec); + + if (i % 2 == 0) + arb_neg(res + i, res + i); + } +} + +void +_arb_poly_gamma_stirling_eval(arb_ptr res, const arb_t z, long n, long num, long prec) +{ + long tlen, qlen; + arb_ptr T, Q; + mag_ptr err; + arb_t c; + + T = _arb_vec_init(num); + Q = _arb_vec_init(num); + err = _mag_vec_init(num); + arb_init(c); + + arb_gamma_stirling_bound(err, z, 0, num, n); + + if (n <= 1) + { + _arb_vec_zero(res, num); + } + else + { + qlen = FLINT_MIN(2 * (n - 1) + 1, num); + tlen = FLINT_MIN(2 * (n - 1) - 1, num); + bsplit(Q, T, z, 1, n, num, prec); + _arb_poly_div_series(res, T, tlen, Q, qlen, num, prec); + } + + /* ((z-1/2) + t) * log(z+t) */ + _arb_poly_log_cpx_series(T, z, num, prec); + arb_one(c); + arb_mul_2exp_si(c, c, -1); + arb_sub(c, z, c, prec); + _arb_poly_mullow_cpx(T, T, num, c, num, prec); + + /* constant term */ + arb_const_log_sqrt2pi(c, prec); + arb_add(T, T, c, prec); + + /* subtract (z+t) */ + arb_sub(T, T, z, prec); + if (num > 1) + arb_sub_ui(T + 1, T + 1, 1, prec); + + _arb_vec_add(res, res, T, num, prec); + + _arb_vec_add_error_mag_vec(res, err, num); + + _arb_vec_clear(T, num); + _arb_vec_clear(Q, num); + _mag_vec_clear(err, num); + arb_clear(c); +} + +void +_arb_poly_gamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long prec) +{ + int reflect; + long i, rflen, r, n, wp; + arb_ptr t, u, v; + arb_struct f[2]; + + hlen = FLINT_MIN(hlen, len); + wp = prec + FLINT_BIT_COUNT(prec); + + t = _arb_vec_init(len); + u = _arb_vec_init(len); + v = _arb_vec_init(len); + arb_init(f); + arb_init(f + 1); + + /* use zeta values at small integers */ + if (arb_is_int(h) && (arf_cmpabs_ui(arb_midref(h), prec / 2) < 0)) + { + r = arf_get_si(arb_midref(h), FMPR_RND_DOWN); + + if (r <= 0) + { + _arb_vec_indeterminate(v, len); + } + else if (r == 1) + { + _arb_poly_lgamma_series_at_one(u, len, wp); + _arb_poly_exp_series(v, u, len, len, wp); + } + else + { + _arb_poly_lgamma_series_at_one(u, len, wp); + _arb_poly_exp_series(t, u, len, len, wp); + arb_one(f); + arb_one(f + 1); + rflen = FLINT_MIN(len, r); + _arb_poly_rising_ui_series(u, f, FLINT_MIN(2, len), r - 1, rflen, wp); + _arb_poly_mullow(v, t, len, u, rflen, len, wp); + } + } + else + { + /* otherwise use Stirling series */ + arb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); + + /* gamma(h) = (rf(1-h, r) * pi) / (gamma(1-h+r) sin(pi h)), h = h0 + t*/ + if (reflect) + { + /* u = 1/gamma(r+1-h) */ + arb_sub_ui(f, h, r + 1, wp); + arb_neg(f, f); + _arb_poly_gamma_stirling_eval(t, f, n, len, wp); + _arb_vec_neg(t, t, len); + _arb_poly_exp_series(u, t, len, len, wp); + for (i = 1; i < len; i += 2) + arb_neg(u + i, u + i); + + /* v = 1/sin(pi x) */ + arb_const_pi(f + 1, wp); + arb_mul(f, h, f + 1, wp); + _arb_poly_sin_series(t, f, 2, len, wp); + _arb_poly_inv_series(v, t, len, len, wp); + + _arb_poly_mullow(t, u, len, v, len, len, wp); + + /* rf(1-h,r) * pi */ + if (r == 0) + { + rflen = 1; + arb_const_pi(u, wp); + } + else + { + arb_sub_ui(f, h, 1, wp); + arb_neg(f, f); + arb_set_si(f + 1, -1); + rflen = FLINT_MIN(len, r + 1); + _arb_poly_rising_ui_series(u, f, FLINT_MIN(2, len), r, rflen, wp); + arb_const_pi(v, wp); + _arb_vec_scalar_mul(u, u, rflen, v, wp); + } + + /* multiply by rising factorial */ + _arb_poly_mullow(v, t, len, u, rflen, len, wp); + } + else + { + /* gamma(h) = gamma(h+r) / rf(h,r) */ + if (r == 0) + { + arb_add_ui(f, h, r, wp); + _arb_poly_gamma_stirling_eval(t, f, n, len, wp); + _arb_poly_exp_series(v, t, len, len, wp); + } + else + { + /* TODO: div_series may be better (once it has a good basecase), + if the rising factorial is short */ + arb_set(f, h); + arb_one(f + 1); + rflen = FLINT_MIN(len, r + 1); + _arb_poly_rising_ui_series(u, f, FLINT_MIN(2, len), r, rflen, wp); + _arb_poly_inv_series(t, u, rflen, len, wp); + + arb_add_ui(f, h, r, wp); + _arb_poly_gamma_stirling_eval(v, f, n, len, wp); + _arb_poly_exp_series(u, v, len, len, wp); + + _arb_poly_mullow(v, u, len, t, len, len, wp); + } + } + } + + /* compose with nonconstant part */ + arb_zero(t); + _arb_vec_set(t + 1, h + 1, hlen - 1); + _arb_poly_compose_series(res, v, len, t, hlen, len, prec); + + arb_clear(f); + arb_clear(f + 1); + _arb_vec_clear(t, len); + _arb_vec_clear(u, len); + _arb_vec_clear(v, len); +} + +void +arb_poly_gamma_series(arb_poly_t res, const arb_poly_t f, long n, long prec) +{ + arb_poly_fit_length(res, n); + + if (f->length == 0 || n == 0) + _arb_vec_indeterminate(res->coeffs, n); + else + _arb_poly_gamma_series(res->coeffs, f->coeffs, f->length, n, prec); + + _arb_poly_set_length(res, n); + _arb_poly_normalise(res); +} + diff --git a/arb_poly/lgamma_series.c b/arb_poly/lgamma_series.c new file mode 100644 index 00000000..4050c67f --- /dev/null +++ b/arb_poly/lgamma_series.c @@ -0,0 +1,137 @@ +/*============================================================================= + + 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 "arb_poly.h" +#include "gamma.h" +#include "zeta.h" + +long arf_get_si(const arf_t x, arf_rnd_t rnd); + +void _arb_poly_lgamma_series_at_one(arb_ptr u, long len, long prec); + +void arb_gamma_stirling_choose_param(int * reflect, long * r, long * n, + const arb_t x, int use_reflect, int digamma, long prec); + +void _arb_poly_gamma_stirling_eval(arb_ptr res, const arb_t z, long n, long num, long prec); + +static __inline__ void +_log_rising_ui_series(arb_ptr t, const arb_t x, long r, long len, long prec) +{ + arb_struct f[2]; + long rflen; + + arb_init(f); + arb_init(f + 1); + arb_set(f, x); + arb_one(f + 1); + + rflen = FLINT_MIN(len, r + 1); + _arb_poly_rising_ui_series(t, f, FLINT_MIN(2, len), r, rflen, prec); + _arb_poly_log_series(t, t, rflen, len, prec); + + arb_clear(f); + arb_clear(f + 1); +} + +void +_arb_poly_lgamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long prec) +{ + int reflect; + long r, n, wp; + arb_t zr; + arb_ptr t, u; + + hlen = FLINT_MIN(hlen, len); + wp = prec + FLINT_BIT_COUNT(prec); + + t = _arb_vec_init(len); + u = _arb_vec_init(len); + arb_init(zr); + + /* use zeta values at small integers */ + if (arb_is_int(h) && (arf_cmpabs_ui(arb_midref(h), prec / 2) < 0)) + { + r = arf_get_si(arb_midref(h), FMPR_RND_DOWN); + + if (r <= 0) + { + _arb_vec_indeterminate(res, len); + } + else + { + _arb_poly_lgamma_series_at_one(u, len, wp); + + if (r != 1) + { + arb_one(zr); + _log_rising_ui_series(t, zr, r - 1, len, wp); + _arb_vec_add(u, u, t, len, wp); + } + } + } + else if (len <= 2) + { + arb_lgamma(u, h, wp); + if (len == 2) + arb_digamma(u + 1, h, wp); + } + else + { + /* otherwise use Stirling series */ + arb_gamma_stirling_choose_param(&reflect, &r, &n, h, 0, 0, wp); + arb_add_ui(zr, h, r, wp); + _arb_poly_gamma_stirling_eval(u, zr, n, len, wp); + + if (r != 0) + { + _log_rising_ui_series(t, h, r, len, wp); + _arb_vec_sub(u, u, t, len, wp); + } + } + + /* compose with nonconstant part */ + arb_zero(t); + _arb_vec_set(t + 1, h + 1, hlen - 1); + _arb_poly_compose_series(res, u, len, t, hlen, len, prec); + + arb_clear(zr); + _arb_vec_clear(t, len); + _arb_vec_clear(u, len); +} + +void +arb_poly_lgamma_series(arb_poly_t res, const arb_poly_t f, long n, long prec) +{ + arb_poly_fit_length(res, n); + + if (f->length == 0 || n == 0) + _arb_vec_indeterminate(res->coeffs, n); + else + _arb_poly_lgamma_series(res->coeffs, f->coeffs, f->length, n, prec); + + _arb_poly_set_length(res, n); + _arb_poly_normalise(res); +} + diff --git a/arb_poly/rgamma_series.c b/arb_poly/rgamma_series.c new file mode 100644 index 00000000..dcede0ea --- /dev/null +++ b/arb_poly/rgamma_series.c @@ -0,0 +1,190 @@ +/*============================================================================= + + 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 "arb_poly.h" +#include "gamma.h" +#include "zeta.h" + +long arf_get_si(const arf_t x, arf_rnd_t rnd); + +void _arb_poly_lgamma_series_at_one(arb_ptr u, long len, long prec); + +void arb_gamma_stirling_choose_param(int * reflect, long * r, long * n, + const arb_t x, int use_reflect, int digamma, long prec); + +void _arb_poly_gamma_stirling_eval(arb_ptr res, const arb_t z, long n, long num, long prec); + + +void +_arb_poly_rgamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long prec) +{ + int reflect; + long i, rflen, r, n, wp; + arb_ptr t, u, v; + arb_struct f[2]; + + hlen = FLINT_MIN(hlen, len); + wp = prec + FLINT_BIT_COUNT(prec); + + t = _arb_vec_init(len); + u = _arb_vec_init(len); + v = _arb_vec_init(len); + arb_init(f); + arb_init(f + 1); + + /* use zeta values at small integers */ + if (arb_is_int(h) && (arf_cmpabs_ui(arb_midref(h), prec / 2) < 0)) + { + r = arf_get_si(arb_midref(h), FMPR_RND_DOWN); + + _arb_poly_lgamma_series_at_one(u, len, wp); + + _arb_vec_neg(u, u, len); + _arb_poly_exp_series(t, u, len, len, wp); + + if (r == 1) + { + _arb_vec_swap(v, t, len); + } + else if (r <= 0) + { + arb_set(f, h); + arb_one(f + 1); + rflen = FLINT_MIN(len, 2 - r); + _arb_poly_rising_ui_series(u, f, FLINT_MIN(2, len), 1 - r, rflen, wp); + _arb_poly_mullow(v, t, len, u, rflen, len, wp); + } + else + { + arb_one(f); + arb_one(f + 1); + rflen = FLINT_MIN(len, r); + _arb_poly_rising_ui_series(v, f, FLINT_MIN(2, len), r - 1, rflen, wp); + + /* TODO: use div_series? */ + _arb_poly_inv_series(u, v, rflen, len, wp); + _arb_poly_mullow(v, t, len, u, len, len, wp); + } + } + else + { + /* otherwise use Stirling series */ + arb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp); + + /* rgamma(h) = (gamma(1-h+r) sin(pi h)) / (rf(1-h, r) * pi), h = h0 + t*/ + if (reflect) + { + /* u = gamma(r+1-h) */ + arb_sub_ui(f, h, r + 1, wp); + arb_neg(f, f); + _arb_poly_gamma_stirling_eval(t, f, n, len, wp); + _arb_poly_exp_series(u, t, len, len, wp); + for (i = 1; i < len; i += 2) + arb_neg(u + i, u + i); + + /* v = sin(pi x) */ + arb_const_pi(f + 1, wp); + arb_mul(f, h, f + 1, wp); + _arb_poly_sin_series(v, f, 2, len, wp); + + _arb_poly_mullow(t, u, len, v, len, len, wp); + + /* rf(1-h,r) * pi */ + if (r == 0) + { + arb_const_pi(u, wp); + _arb_vec_scalar_div(v, t, len, u, wp); + } + else + { + arb_sub_ui(f, h, 1, wp); + arb_neg(f, f); + arb_set_si(f + 1, -1); + rflen = FLINT_MIN(len, r + 1); + _arb_poly_rising_ui_series(v, f, FLINT_MIN(2, len), r, rflen, wp); + arb_const_pi(u, wp); + _arb_vec_scalar_mul(v, v, rflen, u, wp); + + /* divide by rising factorial */ + /* TODO: might better to use div_series, when it has a good basecase */ + _arb_poly_inv_series(u, v, rflen, len, wp); + _arb_poly_mullow(v, t, len, u, len, len, wp); + } + } + else + { + /* rgamma(h) = rgamma(h+r) rf(h,r) */ + if (r == 0) + { + arb_add_ui(f, h, r, wp); + _arb_poly_gamma_stirling_eval(t, f, n, len, wp); + _arb_vec_neg(t, t, len); + _arb_poly_exp_series(v, t, len, len, wp); + } + else + { + arb_set(f, h); + arb_one(f + 1); + rflen = FLINT_MIN(len, r + 1); + _arb_poly_rising_ui_series(t, f, FLINT_MIN(2, len), r, rflen, wp); + + arb_add_ui(f, h, r, wp); + _arb_poly_gamma_stirling_eval(v, f, n, len, wp); + _arb_vec_neg(v, v, len); + _arb_poly_exp_series(u, v, len, len, wp); + + _arb_poly_mullow(v, u, len, t, rflen, len, wp); + } + } + } + + /* compose with nonconstant part */ + arb_zero(t); + _arb_vec_set(t + 1, h + 1, hlen - 1); + _arb_poly_compose_series(res, v, len, t, hlen, len, prec); + + arb_clear(f); + arb_clear(f + 1); + _arb_vec_clear(t, len); + _arb_vec_clear(u, len); + _arb_vec_clear(v, len); +} + +void +arb_poly_rgamma_series(arb_poly_t res, const arb_poly_t f, long n, long prec) +{ + if (f->length == 0 || n == 0) + { + arb_poly_zero(res); + } + else + { + arb_poly_fit_length(res, n); + _arb_poly_rgamma_series(res->coeffs, f->coeffs, f->length, n, prec); + _arb_poly_set_length(res, n); + _arb_poly_normalise(res); + } +} + diff --git a/arb_poly/rising_ui_series.c b/arb_poly/rising_ui_series.c new file mode 100644 index 00000000..e6e3013a --- /dev/null +++ b/arb_poly/rising_ui_series.c @@ -0,0 +1,119 @@ +/*============================================================================= + + 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, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +static void +_arb_poly_rising_ui_series_bsplit(arb_ptr res, + arb_srcptr f, long flen, ulong a, ulong b, + long trunc, long prec) +{ + flen = FLINT_MIN(flen, trunc); + + if (b - a == 1) + { + arb_add_ui(res, f, a, prec); + _arb_vec_set(res + 1, f + 1, flen - 1); + } + else + { + arb_ptr L, R; + long len1, len2; + + long m = a + (b - a) / 2; + + len1 = poly_pow_length(flen, m - a, trunc); + len2 = poly_pow_length(flen, b - m, trunc); + + L = _arb_vec_init(len1 + len2); + R = L + len1; + + _arb_poly_rising_ui_series_bsplit(L, f, flen, a, m, trunc, prec); + _arb_poly_rising_ui_series_bsplit(R, f, flen, m, b, trunc, prec); + + _arb_poly_mullow(res, L, len1, R, len2, + FLINT_MIN(trunc, len1 + len2 - 1), prec); + + _arb_vec_clear(L, len1 + len2); + } +} + +void +_arb_poly_rising_ui_series(arb_ptr res, + arb_srcptr f, long flen, ulong r, + long trunc, long prec) +{ + if (trunc == 1 || flen == 1) + { + arb_rising_ui(res, f, r, prec); + _arb_vec_zero(res + 1, trunc - 1); + } + else if (trunc == 2) + { + arb_rising2_ui(res, res + 1, f, r, prec); + arb_mul(res + 1, res + 1, f + 1, prec); + } + else + { + _arb_poly_rising_ui_series_bsplit(res, f, flen, 0, r, trunc, prec); + } +} + +void +arb_poly_rising_ui_series(arb_poly_t res, const arb_poly_t f, ulong r, long trunc, long prec) +{ + long len; + + if (f->length == 0 && r != 0) + { + arb_poly_zero(res); + return; + } + + if (r == 0) + { + arb_poly_one(res); + return; + } + + len = poly_pow_length(f->length, r, trunc); + + if (f == res) + { + arb_poly_t tmp; + arb_poly_init(tmp); + arb_poly_rising_ui_series(tmp, f, r, len, prec); + arb_poly_swap(tmp, res); + arb_poly_clear(tmp); + } + else + { + arb_poly_fit_length(res, len); + _arb_poly_rising_ui_series(res->coeffs, f->coeffs, f->length, r, len, prec); + _arb_poly_set_length(res, len); + _arb_poly_normalise(res); + } +} + diff --git a/arb_poly/test/t-gamma_series.c b/arb_poly/test/t-gamma_series.c new file mode 100644 index 00000000..041e9553 --- /dev/null +++ b/arb_poly/test/t-gamma_series.c @@ -0,0 +1,121 @@ +/*============================================================================= + + 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, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("gamma_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 2000; iter++) + { + long m, n1, n2, qbits, rbits1, rbits2, rbits3; + fmpq_poly_t A; + arb_poly_t a, b, c, d; + + qbits = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 400); + rbits2 = 2 + n_randint(state, 400); + rbits3 = 2 + n_randint(state, 400); + + m = 1 + n_randint(state, 25); + n1 = 1 + n_randint(state, 25); + n2 = 1 + n_randint(state, 25); + + fmpq_poly_init(A); + arb_poly_init(a); + arb_poly_init(b); + arb_poly_init(c); + arb_poly_init(d); + + fmpq_poly_randtest_not_zero(A, state, m, qbits); + arb_poly_set_fmpq_poly(a, A, rbits1); + + arb_poly_gamma_series(b, a, n1, rbits2); + arb_poly_gamma_series(c, a, n2, rbits3); + + arb_poly_set(d, b); + arb_poly_truncate(d, FLINT_MIN(n1, n2)); + arb_poly_truncate(c, FLINT_MIN(n1, n2)); + + if (!arb_poly_overlaps(c, d)) + { + printf("FAIL\n\n"); + printf("n1 = %ld, n2 = %ld, bits2 = %ld, bits3 = %ld\n", n1, n2, rbits2, rbits3); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("a = "); arb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); arb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); arb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + /* check gamma(a) * a = gamma(a+1) */ + arb_poly_mullow(c, b, a, n1, rbits2); + + arb_poly_set(d, a); + arb_add_ui(d->coeffs, d->coeffs, 1, rbits2); + arb_poly_gamma_series(d, d, n1, rbits2); + + if (!arb_poly_overlaps(c, d)) + { + printf("FAIL (functional equation, n1 = %ld)\n\n", n1); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("a = "); arb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); arb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); arb_poly_printd(c, 15); printf("\n\n"); + printf("d = "); arb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + arb_poly_gamma_series(a, a, n1, rbits2); + if (!arb_poly_overlaps(a, b)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + arb_poly_clear(a); + arb_poly_clear(b); + arb_poly_clear(c); + arb_poly_clear(d); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb_poly/test/t-lgamma_series.c b/arb_poly/test/t-lgamma_series.c new file mode 100644 index 00000000..bc4d5db6 --- /dev/null +++ b/arb_poly/test/t-lgamma_series.c @@ -0,0 +1,123 @@ +/*============================================================================= + + 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, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("lgamma_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + long m, n1, n2, qbits, rbits1, rbits2, rbits3; + fmpq_poly_t A; + arb_poly_t a, b, c, d; + + qbits = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 400); + rbits2 = 2 + n_randint(state, 400); + rbits3 = 2 + n_randint(state, 400); + + m = 1 + n_randint(state, 30); + n1 = 1 + n_randint(state, 30); + n2 = 1 + n_randint(state, 30); + + fmpq_poly_init(A); + arb_poly_init(a); + arb_poly_init(b); + arb_poly_init(c); + arb_poly_init(d); + + fmpq_poly_randtest_not_zero(A, state, m, qbits); + fmpz_abs(A->coeffs, A->coeffs); + arb_poly_set_fmpq_poly(a, A, rbits1); + + arb_poly_lgamma_series(b, a, n1, rbits2); + arb_poly_lgamma_series(c, a, n2, rbits3); + + arb_poly_set(d, b); + arb_poly_truncate(d, FLINT_MIN(n1, n2)); + arb_poly_truncate(c, FLINT_MIN(n1, n2)); + + if (!arb_poly_overlaps(c, d)) + { + printf("FAIL\n\n"); + printf("n1 = %ld, n2 = %ld, bits2 = %ld, bits3 = %ld\n", n1, n2, rbits2, rbits3); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("a = "); arb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); arb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); arb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + /* check loggamma(a) + log(a) = loggamma(a+1) */ + arb_poly_log_series(c, a, n1, rbits2); + arb_poly_add(c, b, c, rbits2); + + arb_poly_set(d, a); + arb_add_ui(d->coeffs, d->coeffs, 1, rbits2); + arb_poly_lgamma_series(d, d, n1, rbits2); + + if (!arb_poly_overlaps(c, d)) + { + printf("FAIL (functional equation)\n\n"); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("a = "); arb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); arb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); arb_poly_printd(c, 15); printf("\n\n"); + printf("d = "); arb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + arb_poly_lgamma_series(a, a, n1, rbits2); + if (!arb_poly_overlaps(a, b)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + arb_poly_clear(a); + arb_poly_clear(b); + arb_poly_clear(c); + arb_poly_clear(d); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb_poly/test/t-rgamma_series.c b/arb_poly/test/t-rgamma_series.c new file mode 100644 index 00000000..15882eb5 --- /dev/null +++ b/arb_poly/test/t-rgamma_series.c @@ -0,0 +1,120 @@ +/*============================================================================= + + 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, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("rgamma_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 2000; iter++) + { + long m, n1, n2, qbits, rbits1, rbits2, rbits3; + fmpq_poly_t A; + arb_poly_t a, b, c, d; + + qbits = 2 + n_randint(state, 200); + rbits1 = 2 + n_randint(state, 400); + rbits2 = 2 + n_randint(state, 400); + rbits3 = 2 + n_randint(state, 400); + + m = 1 + n_randint(state, 30); + n1 = 1 + n_randint(state, 30); + n2 = 1 + n_randint(state, 30); + + fmpq_poly_init(A); + arb_poly_init(a); + arb_poly_init(b); + arb_poly_init(c); + arb_poly_init(d); + + fmpq_poly_randtest_not_zero(A, state, m, qbits); + arb_poly_set_fmpq_poly(a, A, rbits1); + + arb_poly_rgamma_series(b, a, n1, rbits2); + arb_poly_rgamma_series(c, a, n2, rbits3); + + arb_poly_set(d, b); + arb_poly_truncate(d, FLINT_MIN(n1, n2)); + arb_poly_truncate(c, FLINT_MIN(n1, n2)); + + if (!arb_poly_overlaps(c, d)) + { + printf("FAIL\n\n"); + printf("n1 = %ld, n2 = %ld, bits2 = %ld, bits3 = %ld\n", n1, n2, rbits2, rbits3); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("a = "); arb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); arb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); arb_poly_printd(c, 15); printf("\n\n"); + + abort(); + } + + /* check rgamma(a) = a gamma(a+1) */ + arb_poly_set(d, a); + arb_add_ui(d->coeffs, d->coeffs, 1, rbits2); + arb_poly_rgamma_series(d, d, n1, rbits2); + arb_poly_mullow(c, d, a, n1, rbits2); + + if (!arb_poly_overlaps(b, c)) + { + printf("FAIL (functional equation, n1 = %ld)\n\n", n1); + + printf("A = "); fmpq_poly_print(A); printf("\n\n"); + printf("a = "); arb_poly_printd(a, 15); printf("\n\n"); + printf("b = "); arb_poly_printd(b, 15); printf("\n\n"); + printf("c = "); arb_poly_printd(c, 15); printf("\n\n"); + printf("d = "); arb_poly_printd(d, 15); printf("\n\n"); + + abort(); + } + + arb_poly_rgamma_series(a, a, n1, rbits2); + if (!arb_poly_overlaps(a, b)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + fmpq_poly_clear(A); + arb_poly_clear(a); + arb_poly_clear(b); + arb_poly_clear(c); + arb_poly_clear(d); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb_poly/test/t-rising_ui_series.c b/arb_poly/test/t-rising_ui_series.c new file mode 100644 index 00000000..87f28245 --- /dev/null +++ b/arb_poly/test/t-rising_ui_series.c @@ -0,0 +1,128 @@ +/*============================================================================= + + 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_poly.h" + + +int main() +{ + long iter; + flint_rand_t state; + + printf("rising_ui_series...."); + fflush(stdout); + + flint_randinit(state); + + /* check rf(f, a) * rf(f + a, b) = rf(f, a + b) */ + for (iter = 0; iter < 1000; iter++) + { + long bits, trunc; + ulong a, b; + arb_poly_t f, g, h1, h2, h1h2, h3; + + bits = 2 + n_randint(state, 200); + trunc = 1 + n_randint(state, 20); + a = n_randint(state, 10); + b = n_randint(state, 10); + + arb_poly_init(f); + arb_poly_init(g); + arb_poly_init(h1); + arb_poly_init(h2); + arb_poly_init(h1h2); + arb_poly_init(h3); + + arb_poly_randtest(f, state, 1 + n_randint(state, 20), bits, 4); + arb_poly_set(g, f); + + /* g = f + 1 */ + if (g->length == 0) + { + arb_poly_fit_length(g, 1); + arb_set_ui(g->coeffs, a); + _arb_poly_set_length(g, 1); + _arb_poly_normalise(g); + } + else + { + arb_add_ui(g->coeffs, g->coeffs, a, bits); + _arb_poly_normalise(g); + } + + arb_poly_rising_ui_series(h1, f, a, trunc, bits); + arb_poly_rising_ui_series(h2, g, b, trunc, bits); + arb_poly_rising_ui_series(h3, f, a + b, trunc, bits); + + arb_poly_mullow(h1h2, h1, h2, trunc, bits); + + if (!arb_poly_overlaps(h1h2, h3)) + { + printf("FAIL\n\n"); + printf("bits = %ld\n", bits); + printf("trunc = %ld\n", trunc); + printf("a = %lu\n", a); + printf("b = %lu\n", a); + + printf("f = "); arb_poly_printd(f, 15); printf("\n\n"); + printf("g = "); arb_poly_printd(g, 15); printf("\n\n"); + printf("h1 = "); arb_poly_printd(h1, 15); printf("\n\n"); + printf("h2 = "); arb_poly_printd(h2, 15); printf("\n\n"); + printf("h1h2 = "); arb_poly_printd(h1h2, 15); printf("\n\n"); + printf("h3 = "); arb_poly_printd(h3, 15); printf("\n\n"); + + abort(); + } + + arb_poly_rising_ui_series(f, f, a, trunc, bits); + + if (!arb_poly_equal(f, h1)) + { + printf("FAIL (aliasing)\n\n"); + + printf("bits = %ld\n", bits); + printf("trunc = %ld\n", trunc); + printf("a = %lu\n", a); + printf("b = %lu\n", a); + + printf("f = "); arb_poly_printd(f, 15); printf("\n\n"); + printf("h1 = "); arb_poly_printd(h1, 15); printf("\n\n"); + + abort(); + } + + arb_poly_clear(f); + arb_poly_clear(g); + arb_poly_clear(h1); + arb_poly_clear(h2); + arb_poly_clear(h1h2); + arb_poly_clear(h3); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/doc/source/acb.rst b/doc/source/acb.rst index 4dfbb70e..04041486 100644 --- a/doc/source/acb.rst +++ b/doc/source/acb.rst @@ -439,13 +439,13 @@ Rising factorials Rising factorials ------------------------------------------------------------------------------- -.. function:: void arb_rising_ui_bs(arb_t z, const arb_t x, ulong n, long prec) +.. function:: void acb_rising_ui_bs(acb_t z, const acb_t x, ulong n, long prec) -.. function:: void arb_rising_ui_rs(arb_t z, const arb_t x, ulong n, ulong step, long prec) +.. function:: void acb_rising_ui_rs(acb_t z, const acb_t x, ulong n, ulong step, long prec) -.. function:: void arb_rising_ui_rec(arb_t z, const arb_t x, ulong n, long prec) +.. function:: void acb_rising_ui_rec(acb_t z, const acb_t x, ulong n, long prec) -.. function:: void arb_rising_ui(arb_t z, const arb_t x, ulong n, long prec) +.. function:: void acb_rising_ui(acb_t z, const acb_t x, ulong n, long prec) Computes the rising factorial `z = x (x+1) (x+2) \cdots (x+n-1)`. diff --git a/doc/source/acb_poly.rst b/doc/source/acb_poly.rst index eeb31dd1..31606ef0 100644 --- a/doc/source/acb_poly.rst +++ b/doc/source/acb_poly.rst @@ -646,6 +646,35 @@ Special functions The underscore version does not support aliasing, and requires the lengths to be nonzero. +.. function:: void _acb_poly_gamma_series(acb_ptr res, acb_srcptr h, long hlen, long n, long prec) + +.. function:: void acb_poly_gamma_series(acb_poly_t res, const acb_poly_t h, long n, long prec) + +.. function:: void _acb_poly_rgamma_series(acb_ptr res, acb_srcptr h, long hlen, long n, long prec) + +.. function:: void acb_poly_rgamma_series(acb_poly_t res, const acb_poly_t h, long n, long prec) + +.. function:: void _acb_poly_lgamma_series(acb_ptr res, acb_srcptr h, long hlen, long n, long prec) + +.. function:: void acb_poly_lgamma_series(acb_poly_t res, const acb_poly_t h, long n, long prec) + + Sets *res* to the series expansion of `\Gamma(h(x))`, `1/\Gamma(h(x))`, + or `\log \Gamma(h(x))`, truncated to length *n*. + + These functions first generate the Taylor series at the constant + term of *h*, and then call :func:`_acb_poly_compose_series`. + The Taylor coefficients are generated using Stirling's series. + + The underscore methods support aliasing of the input and output + arrays, and require that *hlen* and *n* are greater than zero. + +.. function:: void _acb_poly_rising_ui_series(acb_ptr res, acb_srcptr f, long flen, ulong r, long trunc, long prec) + +.. function:: void acb_poly_rising_ui_series(acb_poly_t res, const acb_poly_t f, ulong r, long trunc, long prec) + + Sets *res* to the rising factorial `(f) (f+1) (f+2) \cdots (f+r-1)`, truncated + to length *trunc*. The underscore method assumes that *flen*, *r* and *trunc* + are at least 1, and does not support aliasing. Uses binary splitting. Root-finding ------------------------------------------------------------------------------- diff --git a/doc/source/arb_poly.rst b/doc/source/arb_poly.rst index 345acf5c..3ac8f1e3 100644 --- a/doc/source/arb_poly.rst +++ b/doc/source/arb_poly.rst @@ -796,6 +796,38 @@ Powers and special functions The underscore version does not support aliasing, and requires the lengths to be nonzero. +.. function:: void _arb_poly_gamma_series(arb_ptr res, arb_srcptr h, long hlen, long n, long prec) + +.. function:: void arb_poly_gamma_series(arb_poly_t res, const arb_poly_t h, long n, long prec) + +.. function:: void _arb_poly_rgamma_series(arb_ptr res, arb_srcptr h, long hlen, long n, long prec) + +.. function:: void arb_poly_rgamma_series(arb_poly_t res, const arb_poly_t h, long n, long prec) + +.. function:: void _arb_poly_lgamma_series(arb_ptr res, arb_srcptr h, long hlen, long n, long prec) + +.. function:: void arb_poly_lgamma_series(arb_poly_t res, const arb_poly_t h, long n, long prec) + + Sets *res* to the series expansion of `\Gamma(h(x))`, `1/\Gamma(h(x))`, + or `\log \Gamma(h(x))`, truncated to length *n*. + + These functions first generate the Taylor series at the constant + term of *h*, and then call :func:`_arb_poly_compose_series`. + The Taylor coefficients are generated using the Riemann zeta function + if the constant term of *h* is a small integer, + and with Stirling's series otherwise. + + The underscore methods support aliasing of the input and output + arrays, and require that *hlen* and *n* are greater than zero. + +.. function:: void _arb_poly_rising_ui_series(arb_ptr res, arb_srcptr f, long flen, ulong r, long trunc, long prec) + +.. function:: void arb_poly_rising_ui_series(arb_poly_t res, const arb_poly_t f, ulong r, long trunc, long prec) + + Sets *res* to the rising factorial `(f) (f+1) (f+2) \cdots (f+r-1)`, truncated + to length *trunc*. The underscore method assumes that *flen*, *r* and *trunc* + are at least 1, and does not support aliasing. Uses binary splitting. + Root-finding -------------------------------------------------------------------------------