diff --git a/Makefile.in b/Makefile.in index bfadbe01..0421296e 100644 --- a/Makefile.in +++ b/Makefile.in @@ -13,7 +13,7 @@ QUIET_AR = @echo ' ' AR ' ' $@; AT=@ BUILD_DIRS = fmpr arf mag arb arb_mat arb_poly arb_calc acb acb_mat acb_poly \ - acb_calc acb_hypgeom acb_modular acb_dirichlet bernoulli hypgeom \ + acb_calc acb_hypgeom acb_modular acb_dirichlet arb_hypgeom bernoulli hypgeom \ fmpz_extras bool_mat partitions \ $(EXTRA_BUILD_DIRS) diff --git a/arb_hypgeom.h b/arb_hypgeom.h new file mode 100644 index 00000000..9f3dbe2a --- /dev/null +++ b/arb_hypgeom.h @@ -0,0 +1,43 @@ +/* + Copyright (C) 2016 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#ifndef ARB_HYPGEOM_H +#define ARB_HYPGEOM_H + +#include "arb.h" +#include "arb_poly.h" + +#ifdef __cplusplus +extern "C" { +#endif + +void arb_hypgeom_erf(arb_t res, const arb_t z, slong prec); +void _arb_hypgeom_erf_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec); +void arb_hypgeom_erf_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec); + +void arb_hypgeom_erfc(arb_t res, const arb_t z, slong prec); +void _arb_hypgeom_erfc_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec); +void arb_hypgeom_erfc_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec); + +void arb_hypgeom_erfi(arb_t res, const arb_t z, slong prec); +void _arb_hypgeom_erfi_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec); +void arb_hypgeom_erfi_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec); + +void arb_hypgeom_fresnel(arb_t res1, arb_t res2, const arb_t z, int normalized, slong prec); +void _arb_hypgeom_fresnel_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, int normalized, slong len, slong prec); +void arb_hypgeom_fresnel_series(arb_poly_t s, arb_poly_t c, const arb_poly_t h, int normalized, slong len, slong prec); + +#ifdef __cplusplus +} +#endif + +#endif + diff --git a/arb_hypgeom/erf_series.c b/arb_hypgeom/erf_series.c new file mode 100644 index 00000000..d44b04d0 --- /dev/null +++ b/arb_hypgeom/erf_series.c @@ -0,0 +1,75 @@ +/* + Copyright (C) 2016 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "arb_hypgeom.h" + +void +_arb_hypgeom_erf_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec) +{ + arb_t c; + arb_init(c); + + arb_hypgeom_erf(c, h, prec); + + hlen = FLINT_MIN(hlen, len); + + if (hlen == 1) + { + _arb_vec_zero(g + 1, len - 1); + } + else + { + arb_ptr t, u; + slong ulen; + + t = _arb_vec_init(len); + u = _arb_vec_init(len); + + /* erf(h(x)) = integral(h'(x) exp(-h(x)^2)) * (2/sqrt(pi)) */ + ulen = FLINT_MIN(len, 2 * hlen - 1); + + _arb_poly_mullow(u, h, hlen, h, hlen, ulen, prec); + _arb_vec_neg(u, u, ulen); + _arb_poly_exp_series(u, u, ulen, len, prec); + _arb_poly_derivative(t, h, hlen, prec); + _arb_poly_mullow(g, u, len, t, hlen - 1, len, prec); + _arb_poly_integral(g, g, len, prec); + + arb_const_sqrt_pi(t, prec); + arb_inv(t, t, prec); + arb_mul_2exp_si(t, t, 1); + _arb_vec_scalar_mul(g, g, len, t, prec); + + _arb_vec_clear(t, len); + _arb_vec_clear(u, len); + } + + arb_swap(g, c); + arb_clear(c); +} + +void +arb_hypgeom_erf_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec) +{ + slong hlen = h->length; + + if (hlen == 0 || len == 0) + { + arb_poly_zero(g); + return; + } + + arb_poly_fit_length(g, len); + _arb_hypgeom_erf_series(g->coeffs, h->coeffs, hlen, len, prec); + _arb_poly_set_length(g, len); + _arb_poly_normalise(g); +} + diff --git a/arb_hypgeom/erfc_series.c b/arb_hypgeom/erfc_series.c new file mode 100644 index 00000000..296de963 --- /dev/null +++ b/arb_hypgeom/erfc_series.c @@ -0,0 +1,82 @@ +/* + Copyright (C) 2016 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "arb_hypgeom.h" + +void +_arb_hypgeom_erfc_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec) +{ + arb_t c; + arb_init(c); + + arb_hypgeom_erfc(c, h, prec); + + hlen = FLINT_MIN(hlen, len); + + if (hlen == 1) + { + _arb_vec_zero(g + 1, len - 1); + } + else + { + arb_ptr t, u; + slong ulen; + + t = _arb_vec_init(len); + u = _arb_vec_init(len); + + /* erfc(h(x)) = integral(-h'(x) exp(-h(x)^2)) * (2/sqrt(pi)) */ + ulen = FLINT_MIN(len, 2 * hlen - 1); + + _arb_poly_mullow(u, h, hlen, h, hlen, ulen, prec); + _arb_vec_neg(u, u, ulen); + _arb_poly_exp_series(u, u, ulen, len, prec); + _arb_poly_derivative(t, h, hlen, prec); + _arb_poly_mullow(g, u, len, t, hlen - 1, len, prec); + _arb_poly_integral(g, g, len, prec); + + arb_const_sqrt_pi(t, prec); + arb_inv(t, t, prec); + arb_mul_2exp_si(t, t, 1); + _arb_vec_scalar_mul(g, g, len, t, prec); + _arb_vec_neg(g, g, len); + + _arb_vec_clear(t, len); + _arb_vec_clear(u, len); + } + + arb_swap(g, c); + arb_clear(c); +} + +void +arb_hypgeom_erfc_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec) +{ + slong hlen = h->length; + + if (len == 0) + { + arb_poly_zero(g); + return; + } + + if (hlen == 0) + { + arb_poly_one(g); + return; + } + + arb_poly_fit_length(g, len); + _arb_hypgeom_erfc_series(g->coeffs, h->coeffs, hlen, len, prec); + _arb_poly_set_length(g, len); + _arb_poly_normalise(g); +} + diff --git a/arb_hypgeom/erfi_series.c b/arb_hypgeom/erfi_series.c new file mode 100644 index 00000000..2b17686f --- /dev/null +++ b/arb_hypgeom/erfi_series.c @@ -0,0 +1,74 @@ +/* + Copyright (C) 2016 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "arb_hypgeom.h" + +void +_arb_hypgeom_erfi_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec) +{ + arb_t c; + arb_init(c); + + arb_hypgeom_erfi(c, h, prec); + + hlen = FLINT_MIN(hlen, len); + + if (hlen == 1) + { + _arb_vec_zero(g + 1, len - 1); + } + else + { + arb_ptr t, u; + slong ulen; + + t = _arb_vec_init(len); + u = _arb_vec_init(len); + + /* erfi(h(x)) = integral(h'(x) exp(h(x)^2)) * (2/sqrt(pi)) */ + ulen = FLINT_MIN(len, 2 * hlen - 1); + + _arb_poly_mullow(u, h, hlen, h, hlen, ulen, prec); + _arb_poly_exp_series(u, u, ulen, len, prec); + _arb_poly_derivative(t, h, hlen, prec); + _arb_poly_mullow(g, u, len, t, hlen - 1, len, prec); + _arb_poly_integral(g, g, len, prec); + + arb_const_sqrt_pi(t, prec); + arb_inv(t, t, prec); + arb_mul_2exp_si(t, t, 1); + _arb_vec_scalar_mul(g, g, len, t, prec); + + _arb_vec_clear(t, len); + _arb_vec_clear(u, len); + } + + arb_swap(g, c); + arb_clear(c); +} + +void +arb_hypgeom_erfi_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec) +{ + slong hlen = h->length; + + if (hlen == 0 || len == 0) + { + arb_poly_zero(g); + return; + } + + arb_poly_fit_length(g, len); + _arb_hypgeom_erfi_series(g->coeffs, h->coeffs, hlen, len, prec); + _arb_poly_set_length(g, len); + _arb_poly_normalise(g); +} + diff --git a/arb_hypgeom/fresnel_series.c b/arb_hypgeom/fresnel_series.c new file mode 100644 index 00000000..622cf7be --- /dev/null +++ b/arb_hypgeom/fresnel_series.c @@ -0,0 +1,112 @@ +/* + Copyright (C) 2016 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "arb_hypgeom.h" + +void +_arb_hypgeom_fresnel_series(arb_ptr s, arb_ptr c, + arb_srcptr h, slong hlen, int normalized, slong len, slong prec) +{ + hlen = FLINT_MIN(hlen, len); + + if (hlen == 1) + { + arb_hypgeom_fresnel(s, c, h, normalized, prec); + + if (s != NULL) _arb_vec_zero(s + 1, len - 1); + if (c != NULL) _arb_vec_zero(c + 1, len - 1); + } + else + { + arb_t s0, c0; + arb_ptr t, u, v; + slong ulen; + + arb_init(s0); + arb_init(c0); + + arb_hypgeom_fresnel((s != NULL) ? s0 : NULL, + (c != NULL) ? c0 : NULL, h, normalized, prec); + + t = _arb_vec_init(len); + u = _arb_vec_init(len); + v = _arb_vec_init(len); + + /* normalized: */ + /* C(h(x)) = integral(h'(x) cos(-(pi/2) h(x)^2)) */ + /* S(h(x)) = -integral(h'(x) sin(-(pi/2) h(x)^2)) */ + ulen = FLINT_MIN(len, 2 * hlen - 1); + + _arb_poly_mullow(u, h, hlen, h, hlen, ulen, prec); + _arb_vec_neg(u, u, ulen); + + if (normalized) + { + _arb_vec_scalar_mul_2exp_si(u, u, ulen, -1); + _arb_poly_sin_cos_pi_series(u, v, u, ulen, len, prec); + } + else + { + _arb_poly_sin_cos_series(u, v, u, ulen, len, prec); + } + + _arb_poly_derivative(t, h, hlen, prec); + + if (s != NULL) + { + _arb_poly_mullow(s, u, len, t, hlen - 1, len, prec); + _arb_poly_integral(s, s, len, prec); + _arb_vec_neg(s, s, len); + arb_swap(s, s0); + } + + if (c != NULL) + { + _arb_poly_mullow(c, v, len, t, hlen - 1, len, prec); + _arb_poly_integral(c, c, len, prec); + arb_swap(c, c0); + } + + _arb_vec_clear(t, len); + _arb_vec_clear(u, len); + _arb_vec_clear(v, len); + + arb_clear(s0); + arb_clear(c0); + } +} + +void +arb_hypgeom_fresnel_series(arb_poly_t s, arb_poly_t c, + const arb_poly_t h, int normalized, slong len, slong prec) +{ + slong hlen = h->length; + + if (hlen == 0 || len == 0) + { + if (s != NULL) arb_poly_zero(s); + if (c != NULL) arb_poly_zero(c); + return; + } + + if (s != NULL) arb_poly_fit_length(s, len); + if (c != NULL) arb_poly_fit_length(c, len); + + _arb_hypgeom_fresnel_series((s != NULL) ? s->coeffs : NULL, + (c != NULL) ? c->coeffs : NULL, + h->coeffs, hlen, normalized, len, prec); + + if (s != NULL) _arb_poly_set_length(s, len); + if (c != NULL) _arb_poly_set_length(c, len); + if (s != NULL) _arb_poly_normalise(s); + if (c != NULL) _arb_poly_normalise(c); +} + diff --git a/arb_hypgeom/wrappers.c b/arb_hypgeom/wrappers.c new file mode 100644 index 00000000..bbb1d23a --- /dev/null +++ b/arb_hypgeom/wrappers.c @@ -0,0 +1,90 @@ +/* + Copyright (C) 2016 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "arb_hypgeom.h" +#include "acb_hypgeom.h" + +void +arb_hypgeom_erf(arb_t res, const arb_t z, slong prec) +{ + if (!arb_is_finite(z)) + { + arb_indeterminate(res); + } + else + { + acb_t t; + acb_init(t); + arb_set(acb_realref(t), z); + acb_hypgeom_erf(t, t, prec); + arb_swap(res, acb_realref(t)); + acb_clear(t); + } +} + +void +arb_hypgeom_erfc(arb_t res, const arb_t z, slong prec) +{ + if (!arb_is_finite(z)) + { + arb_indeterminate(res); + } + else + { + acb_t t; + acb_init(t); + arb_set(acb_realref(t), z); + acb_hypgeom_erfc(t, t, prec); + arb_swap(res, acb_realref(t)); + acb_clear(t); + } +} + +void +arb_hypgeom_erfi(arb_t res, const arb_t z, slong prec) +{ + if (!arb_is_finite(z)) + { + arb_indeterminate(res); + } + else + { + acb_t t; + acb_init(t); + arb_set(acb_realref(t), z); + acb_hypgeom_erfi(t, t, prec); + arb_swap(res, acb_realref(t)); + acb_clear(t); + } +} + +void +arb_hypgeom_fresnel(arb_t res1, arb_t res2, const arb_t z, int normalized, slong prec) +{ + if (!arb_is_finite(z)) + { + if (res1 != NULL) arb_indeterminate(res1); + if (res2 != NULL) arb_indeterminate(res2); + } + else + { + acb_t t, u; + acb_init(t); + acb_init(u); + arb_set(acb_realref(t), z); + acb_hypgeom_fresnel(res1 ? t : NULL, res2 ? u : NULL, t, normalized, prec); + if (res1 != NULL) arb_swap(res1, acb_realref(t)); + if (res2 != NULL) arb_swap(res2, acb_realref(u)); + acb_clear(t); + acb_clear(u); + } +} +