Skip to content

Commit

Permalink
initial code for hypergeometric 2F1 (some cases still missing)
Browse files Browse the repository at this point in the history
  • Loading branch information
@fredrik-johansson
fredrik-johansson committed Oct 17, 2015
1 parent 3ba474d commit 71b194d
Show file tree
Hide file tree
Showing 7 changed files with 742 additions and 0 deletions.
8 changes: 8 additions & 0 deletions acb_hypgeom.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,6 +140,14 @@ 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_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);
void acb_hypgeom_2f1_pfaff(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);
void acb_hypgeom_2f1_inf(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);
void acb_hypgeom_2f1_inf_limit(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);
void acb_hypgeom_2f1(acb_t res, const acb_t a, const acb_t b, const acb_t c, const acb_t z, int regularized, long prec);

#ifdef __cplusplus
}
#endif
Expand Down
132 changes: 132 additions & 0 deletions acb_hypgeom/2f1.c
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
/*=============================================================================
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"

void
acb_hypgeom_2f1(acb_t res, const acb_t a, const acb_t b,
const acb_t c, const acb_t z, int regularized, long prec)
{
double x, y;

if (!acb_is_finite(a) || !acb_is_finite(b) || !acb_is_finite(c) || !acb_is_finite(z))
{
acb_indeterminate(res);
return;
}

if (acb_is_zero(z))
{
if (regularized)
acb_rgamma(res, c, prec);
else
acb_one(res);
return;
}

/* polynomial */
if (acb_is_int(a) && arf_sgn(arb_midref(acb_realref(a))) <= 0 &&
arf_cmpabs_ui(arb_midref(acb_realref(a)), prec) < 0)
{
acb_hypgeom_2f1_direct(res, a, b, c, z, regularized, prec);
return;
}

/* polynomial */
if (acb_is_int(b) && arf_sgn(arb_midref(acb_realref(b))) <= 0 &&
arf_cmpabs_ui(arb_midref(acb_realref(b)), prec) < 0)
{
acb_hypgeom_2f1_direct(res, a, b, c, z, regularized, prec);
return;
}

/* special value at z = 1 */
if (acb_is_one(z))
{
acb_t t, u, v;

acb_init(t);
acb_init(u);
acb_init(v);

acb_sub(t, c, a, prec);
acb_sub(u, c, b, prec);
acb_sub(v, t, b, prec);

if (arb_is_positive(acb_realref(v)))
{
acb_rgamma(t, t, prec);
acb_rgamma(u, u, prec);
acb_mul(t, t, u, prec);
acb_gamma(v, v, prec);
acb_mul(t, t, v, prec);

if (!regularized)
{
acb_gamma(v, c, prec);
acb_mul(t, t, v, prec);
}

acb_set(res, t);
}
else
{
acb_indeterminate(res);
}

acb_clear(t);
acb_clear(u);
acb_clear(v);

return;
}

x = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
y = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);

x = FLINT_MAX(FLINT_MIN(x, 1e10), -1e10);
y = FLINT_MAX(FLINT_MIN(y, 1e10), -1e10);

if (x*x + y*y < 0.7) /* series in z */
{
acb_hypgeom_2f1_direct(res, a, b, c, z, regularized, prec);
}
else if (x < 0.7 && (x*x + y*y)/((x-1)*(x-1) + y*y) < 0.7) /* series in z/(z-1) */
{
acb_hypgeom_2f1_pfaff(res, a, b, c, z, regularized, prec);
}
else if (x*x + y*y > 1.1) /* series in 1/z */
{
acb_hypgeom_2f1_inf(res, a, b, c, z, regularized, prec);
}
else
{
if (x*x + y*y < 1.0)
acb_hypgeom_2f1_direct(res, a, b, c, z, regularized, prec);
else
acb_hypgeom_2f1_inf(res, a, b, c, z, regularized, prec);
}
}

134 changes: 134 additions & 0 deletions acb_hypgeom/2f1_direct.c
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,134 @@
/*=============================================================================
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"

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)
{
/* 2F1R(a,b,-n,z) = (a)_(n+1) * (b)_(n+1) * z^(n+1) / (n+1)!
* 2F1(a+n+1, b+n+1, n+2, z) */
if (regularized && acb_is_int(c)
&& arf_sgn(arb_midref(acb_realref(c))) <= 0)
{
acb_t n, n1, t, u, v;
acb_ptr aa;
int p, q;

acb_init(n);
acb_init(n1);
acb_init(t);
acb_init(u);
acb_init(v);
aa = _acb_vec_init(4);

acb_neg(n, c);
acb_add_ui(n1, n, 1, prec);

acb_add(aa, a, n1, prec);
acb_add(aa + 1, b, n1, prec);
acb_add_ui(aa + 2, n1, 1, prec);

if (acb_is_one(aa))
{
p = q = 1;
acb_swap(aa, aa + 1);
}
else if (acb_is_one(aa + 1))
{
p = q = 1;
}
else
{
p = q = 2;
acb_one(aa + 3);
}

acb_hypgeom_pfq_direct(t, aa, p, aa + 2, q, z, -1, prec);

/* z^(n+1) */
acb_pow(u, z, n1, prec);
acb_mul(t, t, u, prec);

/* todo: totally wrong??? */
acb_rising(u, a, n1, prec);
acb_mul(t, t, u, prec);

acb_rising(u, b, n1, prec);
acb_mul(t, t, u, prec);

/* 1/(n+1)! */
acb_rgamma(u, aa + 2, prec);
acb_mul(t, t, u, prec);

acb_set(res, t);

acb_clear(n);
acb_clear(n1);
acb_clear(t);
acb_clear(u);
acb_clear(v);
_acb_vec_clear(aa, 4);
}
else
{
acb_ptr aa;
int p, q;

aa = _acb_vec_init(4);

acb_set(aa + 2, c);

if (acb_is_one(a))
{
p = q = 1;
acb_set(aa, b);
}
else if (acb_is_one(b))
{
p = q = 1;
acb_set(aa, a);
}
else
{
p = q = 2;
acb_set(aa, a);
acb_set(aa + 1, b);
acb_one(aa + 3);
}

acb_hypgeom_pfq_direct(res, aa, p, aa + 2, q, z, -1, prec);

if (regularized)
{
acb_rgamma(aa + 2, aa + 2, prec);
acb_mul(res, res, aa + 2, prec);
}

_acb_vec_clear(aa, 4);
}
}

Loading

0 comments on commit 71b194d

Please sign in to comment.