Skip to content

Commit

Permalink
analytic continuation algorithm for 2F1
Browse files Browse the repository at this point in the history
  • Loading branch information
@fredrik-johansson
fredrik-johansson committed Oct 18, 2015
1 parent 71b194d commit 2d28036
Show file tree
Hide file tree
Showing 3 changed files with 438 additions and 0 deletions.
4 changes: 4 additions & 0 deletions acb_hypgeom.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,6 +140,10 @@ void acb_hypgeom_chi(acb_t res, const acb_t z, long prec);

void acb_hypgeom_li(acb_t res, const acb_t z, int offset, long prec);

void acb_hypgeom_2f1_continuation(acb_t res0, acb_t res1,
const acb_t a, const acb_t b, const acb_t c, const acb_t z0,
const acb_t z1, const acb_t f0, const acb_t f1, long prec);

void acb_hypgeom_2f1_series_direct(acb_poly_t res, const acb_poly_t a, const acb_poly_t b,
const acb_poly_t c, const acb_poly_t z, int regularized, long len, long prec);
void acb_hypgeom_2f1_direct(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);
Expand Down
279 changes: 279 additions & 0 deletions acb_hypgeom/2f1_continuation.c
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,279 @@
/*=============================================================================
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) 2015 Fredrik Johansson
******************************************************************************/

#include "acb_hypgeom.h"

/* Differential equation for F(a,b,c,y+z):
(y+z)(y-1+z) F''(z) + ((y+z)(a+b+1) - c) F'(z) + a b F(z) = 0
Coefficients in the Taylor series are bounded by
A * binomial(N+k, k) * nu^k
using the Cauchy-Kovalevskaya majorant method.
See J. van der Hoeven, "Fast evaluation of holonomic functions near
and in regular singularities"
*/
static void
bound(mag_t A, mag_t nu, mag_t N,
const acb_t a, const acb_t b, const acb_t c, const acb_t y,
const acb_t f0, const acb_t f1)
{
mag_t M0, M1, t, u;
acb_t d;

acb_init(d);
mag_init(M0);
mag_init(M1);
mag_init(t);
mag_init(u);

/* nu = max(1/|y-1|, 1/|y|) = 1/min(|y-1|, |y|) */
acb_get_mag_lower(t, y);
acb_sub_ui(d, y, 1, MAG_BITS);
acb_get_mag_lower(u, d);
mag_min(t, t, u);
mag_one(u);
mag_div(nu, u, t);

/* M0 = 2 nu |ab| */
acb_get_mag(t, a);
acb_get_mag(u, b);
mag_mul(M0, t, u);
mag_mul(M0, M0, nu);
mag_mul_2exp_si(M0, M0, 1);

/* M1 = 2 nu |(a+b+1)y-c| + 2|a+b+1| */
acb_add(d, a, b, MAG_BITS);
acb_add_ui(d, d, 1, MAG_BITS);
acb_get_mag(t, d);
acb_mul(d, d, y, MAG_BITS);
acb_sub(d, d, c, MAG_BITS);
acb_get_mag(u, d);
mag_mul(u, u, nu);
mag_add(M1, t, u);
mag_mul_2exp_si(M1, M1, 1);

/* N = max(sqrt(2 M0), 2 M1) / nu */
mag_mul_2exp_si(M0, M0, 1);
mag_sqrt(M0, M0);
mag_mul_2exp_si(M1, M1, 1);
mag_max(N, M0, M1);
mag_div(N, N, nu);

/* A = max(|f0|, |f1| / (nu (N+1)) */
acb_get_mag(t, f0);
acb_get_mag(u, f1);
mag_div(u, u, nu);
mag_div(u, u, N); /* upper bound for dividing by N+1 */
mag_max(A, t, u);

acb_clear(d);
mag_clear(M0);
mag_clear(M1);
mag_clear(t);
mag_clear(u);
}

static void
mag_add_ui(mag_t y, const mag_t x, ulong k)
{
mag_t t;
mag_init(t); /* no need to free */
mag_set_ui(t, k);
mag_add(y, x, t);
}

/*
F(x) = c0 + c1 x + c2 x^2 + c3 x^3 + [...]
F'(x) = c1 + 2 c2 x + 3 c3 x^2 + 4 c4 x^3 + [...]
*/
static void
evaluate_sum(acb_t res, acb_t res1,
const acb_t a, const acb_t b, const acb_t c, const acb_t y,
const acb_t x, const acb_t f0, const acb_t f1, long num, long prec)
{
acb_t s, s2, w, d, e, xpow, ck, cknext;
long k;

acb_init(s);
acb_init(s2);
acb_init(w);
acb_init(d);
acb_init(e);
acb_init(xpow);
acb_init(ck);
acb_init(cknext);

/* d = (y-1)*y */
acb_sub_ui(d, y, 1, prec);
acb_mul(d, d, y, prec);
acb_one(xpow);

for (k = 0; k < num; k++)
{
if (k == 0)
{
acb_set(ck, f0);
acb_set(cknext, f1);
}
else
{
acb_add_ui(w, b, k-1, prec);
acb_mul(w, w, ck, prec);
acb_add_ui(e, a, k-1, prec);
acb_mul(w, w, e, prec);

acb_add(e, a, b, prec);
acb_add_ui(e, e, 2*(k+1)-3, prec);
acb_mul(e, e, y, prec);
acb_sub(e, e, c, prec);
acb_sub_ui(e, e, k-1, prec);
acb_mul_ui(e, e, k, prec);
acb_addmul(w, e, cknext, prec);

acb_mul_ui(e, d, k+1, prec);
acb_mul_ui(e, e, k, prec);
acb_div(w, w, e, prec);
acb_neg(w, w);

acb_set(ck, cknext);
acb_set(cknext, w);
}

acb_addmul(s, ck, xpow, prec);
acb_mul_ui(w, cknext, k+1, prec);
acb_addmul(s2, w, xpow, prec);

acb_mul(xpow, xpow, x, prec);
}

acb_set(res, s);
acb_set(res1, s2);

acb_clear(s);
acb_clear(s2);
acb_clear(w);
acb_clear(d);
acb_clear(e);
acb_clear(xpow);
acb_clear(ck);
acb_clear(cknext);
}

void
acb_hypgeom_2f1_continuation(acb_t res, acb_t res1,
const acb_t a, const acb_t b, const acb_t c, const acb_t y,
const acb_t z, const acb_t f0, const acb_t f1, long prec)
{
mag_t A, nu, N, w, err, err1, R, T, goal;
acb_t x;
long j, k;

mag_init(A);
mag_init(nu);
mag_init(N);
mag_init(err);
mag_init(err1);
mag_init(w);
mag_init(R);
mag_init(T);
mag_init(goal);
acb_init(x);

bound(A, nu, N, a, b, c, y, f0, f1);

acb_sub(x, z, y, prec);

/* |T(k)| <= A * binomial(N+k, k) * nu^k * |x|^k */
acb_get_mag(w, x);
mag_mul(w, w, nu); /* w = nu |x| */
mag_mul_2exp_si(goal, A, -prec-2);

/* bound for T(0) */
mag_set(T, A);
mag_inf(R);

for (k = 1; k < 100 * prec; k++)
{
/* T(k) = T(k) * R(k), R(k) = (N+k)/k * w = (1 + N/k) w */
mag_div_ui(R, N, k);
mag_add_ui(R, R, 1);
mag_mul(R, R, w);

/* T(k) */
mag_mul(T, T, R);

if (mag_cmp(T, goal) <= 0 && mag_cmp_2exp_si(R, 0) < 0)
break;
}

/* T(k) [1 + R + R^2 + R^3 + ...] */
mag_geom_series(err, R, 0);
mag_mul(err, T, err);

/* Now compute T, R for the derivative */
/* Coefficients are A * (k+1) * binomial(N+k+1, k+1) */
mag_add_ui(T, N, 1);
mag_mul(T, T, A);
mag_inf(R);

for (j = 1; j <= k; j++)
{
mag_add_ui(R, N, k + 1);
mag_div_ui(R, R, k);
mag_mul(R, R, w);
mag_mul(T, T, R);
}

mag_geom_series(err1, R, 0);
mag_mul(err1, T, err1);

if (mag_is_inf(err))
{
acb_indeterminate(res);
acb_indeterminate(res1);
}
else
{
evaluate_sum(res, res1, a, b, c, y, x, f0, f1, k, prec);

acb_add_error_mag(res, err);
acb_add_error_mag(res1, err1);
}

mag_clear(A);
mag_clear(nu);
mag_clear(N);
mag_clear(err);
mag_clear(err1);
mag_clear(w);
mag_clear(R);
mag_clear(T);
mag_clear(goal);
acb_clear(x);
}

Loading

0 comments on commit 2d28036

Please sign in to comment.