hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

port gamma function power series methods

Browse source
This commit is contained in:
Fredrik Johansson 2014-06-05 00:43:30 +02:00
parent 66d5b62262
commit e7f0fff490
23 changed files with 2556 additions and 15 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

13
acb.h
View file

@ -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); 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 static __inline__ void
_acb_vec_indeterminate(acb_ptr vec, long len) _acb_vec_indeterminate(acb_ptr vec, long len)
{ {

7
acb_poly.h
View file

@ -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); 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_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); 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); 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_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); void acb_poly_zeta_series(acb_poly_t res, const acb_poly_t f, const acb_t a, int deflate, long n, long prec);

305
acb_poly/gamma_series.c Normal file
View file

@ -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);
}

124
acb_poly/lgamma_series.c Normal file
View file

@ -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);
}

156
acb_poly/rgamma_series.c Normal file
View file

@ -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);
}
}

120
acb_poly/rising_ui_series.c Normal file
View file

@ -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);
}
}

114
acb_poly/test/t-gamma_series.c Normal file
View file

@ -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;
}

115
acb_poly/test/t-lgamma_series.c Normal file
View file

@ -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;
}

113
acb_poly/test/t-rgamma_series.c Normal file
View file

@ -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;
}

128
acb_poly/test/t-rising_ui_series.c Normal file
View file

@ -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;
}

9
arb.h
View file

@ -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 static __inline__ void
_arb_vec_add_error_arf_vec(arb_ptr res, arf_srcptr err, long len) _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); 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 static __inline__ void
_arb_vec_indeterminate(arb_ptr vec, long len) _arb_vec_indeterminate(arb_ptr vec, long len)
{ {

10
arb_poly.h
View file

@ -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, 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); 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_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_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_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_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_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_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_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); 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_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); void arb_poly_zeta_series(arb_poly_t res, const arb_poly_t f, const arb_t a, int deflate, long n, long prec);

350
arb_poly/gamma_series.c Normal file
View file

@ -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);
}

137
arb_poly/lgamma_series.c Normal file
View file

@ -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);
}

190
arb_poly/rgamma_series.c Normal file
View file

@ -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);
}
}

119
arb_poly/rising_ui_series.c Normal file
View file

@ -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);
}
}

121
arb_poly/test/t-gamma_series.c Normal file
View file

@ -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;
}

123
arb_poly/test/t-lgamma_series.c Normal file
View file

@ -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;
}

120
arb_poly/test/t-rgamma_series.c Normal file
View file

@ -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;
}

128
arb_poly/test/t-rising_ui_series.c Normal file
View file

@ -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;
}

8
doc/source/acb.rst
View file

@ -439,13 +439,13 @@ Rising factorials
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)`. Computes the rising factorial `z = x (x+1) (x+2) \cdots (x+n-1)`.

29
doc/source/acb_poly.rst
View file

@ -646,6 +646,35 @@ Special functions
The underscore version does not support aliasing, and requires The underscore version does not support aliasing, and requires
the lengths to be nonzero. 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 Root-finding
------------------------------------------------------------------------------- -------------------------------------------------------------------------------

32
doc/source/arb_poly.rst
View file

@ -796,6 +796,38 @@ Powers and special functions
The underscore version does not support aliasing, and requires The underscore version does not support aliasing, and requires
the lengths to be nonzero. 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 Root-finding
------------------------------------------------------------------------------- -------------------------------------------------------------------------------