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

use recurrence relations for Legendre polynomials of some orders

Browse source
This commit is contained in:
Fredrik Johansson 2015-11-01 15:52:49 +01:00
parent 27a838e66c
commit d84efb61ca
6 changed files with 292 additions and 3 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
acb_hypgeom.h
View file

@ -160,6 +160,7 @@ int acb_hypgeom_2f1_choose(const acb_t z);
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);
void acb_hypgeom_legendre_p_uiui_rec(acb_t res, ulong n, ulong m, const acb_t z, long prec);
void acb_hypgeom_legendre_p(acb_t res, const acb_t n, const acb_t m, const acb_t z, int type, long prec);
void acb_hypgeom_legendre_q(acb_t res, const acb_t n, const acb_t m, const acb_t z, int type, long prec);
void acb_hypgeom_jacobi_p(acb_t res, const acb_t n, const acb_t a, const acb_t b, const acb_t z, long prec);

59
acb_hypgeom/legendre_p.c
View file

@ -31,6 +31,65 @@ acb_hypgeom_legendre_p(acb_t res, const acb_t n, const acb_t m,
{
acb_t a, b, c, w;
if (!acb_is_finite(z))
{
acb_indeterminate(res);
return;
}
if (acb_is_int(n) && acb_is_int(m) && arb_is_nonnegative(acb_realref(n))
&& arb_is_nonnegative(acb_realref(m)) && type == 0)
{
arf_srcptr nn = arb_midref(acb_realref(n));
arf_srcptr mm = arb_midref(acb_realref(m));
if (arf_cmpabs(mm, nn) > 0)
{
acb_zero(res);
return;
}
if (arf_cmpabs_2exp_si(nn, FLINT_BITS - 1) < 0 &&
arf_cmpabs_2exp_si(mm, FLINT_BITS - 1) < 0)
{
long nnn, mmm;
nnn = arf_get_si(nn, ARF_RND_DOWN);
mmm = arf_get_si(mm, ARF_RND_DOWN);
/* we will probably lose all accuracy... */
if (nnn - mmm > 2 * prec)
{
acb_indeterminate(res);
} /* hypergeometric series is faster at high precision, but avoid at 1 */
else if (prec < 500 + nnn * 10.0 || (nnn - mmm) < 10 ||
(arb_contains_si(acb_realref(z), 1) &&
arb_contains_zero(acb_imagref(z))))
{
if (mmm == 0)
{
acb_hypgeom_legendre_p_uiui_rec(res, nnn, mmm, z, prec);
}
else
{
acb_init(a);
acb_init(b);
acb_mul(a, z, z, prec);
acb_sub_ui(a, a, 1, prec);
acb_neg(a, a);
acb_mul_2exp_si(b, m, -1);
acb_pow(a, a, b, prec);
acb_hypgeom_legendre_p_uiui_rec(res, nnn, mmm, z, prec);
acb_mul(res, res, a, prec);
acb_clear(a);
acb_clear(b);
}
return;
}
}
}
acb_init(a);
acb_init(b);
acb_init(c);

89
acb_hypgeom/legendre_p_uiui_rec.c Normal file
View file

@ -0,0 +1,89 @@
/*=============================================================================
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_legendre_p_uiui_rec(acb_t res, ulong n, ulong m, const acb_t z, long prec)
{
acb_t t, u, v;
long k;
if (!acb_is_finite(z))
{
acb_indeterminate(res);
return;
}
if (m > n)
{
acb_zero(res);
return;
}
if ((n - m) / 4 > prec)
{
acb_indeterminate(res);
return;
}
acb_init(t);
acb_init(u);
acb_init(v);
/* t = p(m,m) = (-1)^m (2m-1)!! */
if (m == 0)
arb_one(acb_realref(t));
else
arb_fac2_ui(acb_realref(t), 2 * m - 1, prec);
if (m % 2)
arb_neg(acb_realref(t), acb_realref(t));
if (n > m)
{
/* t = p(m+1,m) = z(2m+1)p(m,m), u = p(m,m) */
acb_mul_ui(u, t, 2 * m + 1, prec);
acb_mul(u, u, z, prec);
acb_swap(t, u);
for (k = m + 2; k <= n; k++)
{
/* t, u = ((2*k-1)*z*t - (k+m-1)*u) / (k-m), t */
acb_mul(v, t, z, prec);
acb_mul_ui(v, v, 2 * k - 1, prec);
acb_submul_ui(v, u, k + m - 1, prec);
acb_div_ui(u, v, k - m, prec);
acb_swap(t, u);
}
}
acb_set(res, t);
acb_clear(t);
acb_clear(u);
acb_clear(v);
}

126
acb_hypgeom/test/t-legendre_p.c Normal file
View file

@ -0,0 +1,126 @@
/*=============================================================================
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"
int main()
{
long iter;
flint_rand_t state;
printf("legendre_p....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 2000; iter++)
{
acb_t n, na, nb, m, z, res1, res2, res3, t, u;
long prec1, prec2, ebits;
int type;
acb_init(n);
acb_init(na);
acb_init(nb);
acb_init(m);
acb_init(z);
acb_init(res1);
acb_init(res2);
acb_init(res3);
acb_init(t);
acb_init(u);
prec1 = 2 + n_randint(state, 300);
prec2 = 2 + n_randint(state, 300);
ebits = 1 + n_randint(state, 10);
if (n_randint(state, 2))
{
acb_set_si(m, n_randint(state, 20) - 10);
acb_set_si(n, n_randint(state, 20) - 10);
}
else
{
acb_randtest_param(n, state, 1 + n_randint(state, 400), ebits);
acb_randtest_param(m, state, 1 + n_randint(state, 400), ebits);
}
acb_randtest_param(z, state, 1 + n_randint(state, 400), ebits);
acb_sub_ui(na, n, 1, prec2);
acb_add_ui(nb, n, 1, prec2);
type = n_randint(state, 2);
acb_hypgeom_legendre_p(res1, n, m, z, type, prec1);
acb_hypgeom_legendre_p(res2, na, m, z, type, prec2);
acb_hypgeom_legendre_p(res3, nb, m, z, type, prec2);
acb_add(t, n, m, prec2);
acb_mul(t, t, res2, prec2);
acb_sub(u, n, m, prec2);
acb_add_ui(u, u, 1, prec2);
acb_mul(u, u, res3, prec2);
acb_add(t, t, u, prec2);
acb_mul_2exp_si(u, n, 1);
acb_add_ui(u, u, 1, prec2);
acb_mul(u, u, z, prec2);
acb_mul(u, u, res1, prec2);
if (!acb_overlaps(t, u))
{
printf("FAIL: consistency\n\n");
printf("iter = %ld, prec1 = %ld, prec2 = %ld\n\n", iter, prec1, prec2);
printf("type = %d\n\n", type);
printf("n = "); acb_printd(n, 30); printf("\n\n");
printf("m = "); acb_printd(m, 30); printf("\n\n");
printf("z = "); acb_printd(z, 30); printf("\n\n");
printf("res1 = "); acb_printd(res1, 30); printf("\n\n");
printf("res2 = "); acb_printd(res2, 30); printf("\n\n");
printf("res3 = "); acb_printd(res3, 30); printf("\n\n");
printf("t = "); acb_printd(t, 30); printf("\n\n");
printf("u = "); acb_printd(u, 30); printf("\n\n");
abort();
}
acb_clear(n);
acb_clear(na);
acb_clear(nb);
acb_clear(m);
acb_clear(z);
acb_clear(res1);
acb_clear(res2);
acb_clear(res3);
acb_clear(t);
acb_clear(u);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

15
acb_hypgeom/test/t-legendre_q.c
View file

@ -91,7 +91,7 @@ int main()
acb_clear(res2);
}
for (iter = 0; iter < 1000; iter++)
for (iter = 0; iter < 2000; iter++)
{
acb_t n, m, z, res1, res2, t, u;
long prec1, prec2, ebits;
@ -109,8 +109,17 @@ int main()
prec2 = 2 + n_randint(state, 300);
ebits = 1 + n_randint(state, 10);
acb_randtest_param(n, state, 1 + n_randint(state, 400), ebits);
acb_randtest_param(m, state, 1 + n_randint(state, 400), ebits);
if (n_randint(state, 2))
{
acb_set_si(m, n_randint(state, 20) - 10);
acb_set_si(n, n_randint(state, 20) - 10);
}
else
{
acb_randtest_param(n, state, 1 + n_randint(state, 400), ebits);
acb_randtest_param(m, state, 1 + n_randint(state, 400), ebits);
}
acb_randtest_param(z, state, 1 + n_randint(state, 400), ebits);
type = n_randint(state, 2);

5
doc/source/acb_hypgeom.rst
View file

@ -782,3 +782,8 @@ Orthogonal polynomials and functions
.. [WQ3c] http://functions.wolfram.com/07.12.26.0003.01
.. [WQ3d] http://functions.wolfram.com/07.12.26.0088.01
.. function:: void acb_hypgeom_legendre_p_uiui_rec(acb_t res, ulong n, ulong m, const acb_t z, long prec)
For nonnegative integer *n* and *m*, uses recurrence relations to evaluate
`(1-z^2)^{-m/2} P_n^m(z)` which is a polynomial in *z*.