wrap airy functions for arb_hypgeom

arb_hypgeom.h

@@ -68,6 +68,13 @@ void arb_hypgeom_li(arb_t res, const arb_t z, int offset, slong prec);
 void _arb_hypgeom_li_series(arb_ptr g, arb_srcptr h, slong hlen, int offset, slong len, slong prec);
 void arb_hypgeom_li_series(arb_poly_t g, const arb_poly_t h, int offset, slong len, slong prec);
 
+void arb_hypgeom_airy(arb_t ai, arb_t aip, arb_t bi, arb_t bip, const arb_t z, slong prec);
+void arb_hypgeom_airy_jet(arb_ptr ai, arb_ptr bi, const arb_t z, slong len, slong prec);
+void arb_hypgeom_airy_series(arb_poly_t ai, arb_poly_t ai_prime,
+    arb_poly_t bi, arb_poly_t bi_prime, const arb_poly_t z, slong len, slong prec);
+void _arb_hypgeom_airy_series(arb_ptr ai, arb_ptr ai_prime,
+    arb_ptr bi, arb_ptr bi_prime, arb_srcptr z, slong zlen, slong len, slong prec);
+
 #ifdef __cplusplus
 }
 #endif

arb_hypgeom/airy.c

@@ -0,0 +1,46 @@
+/*
+    Copyright (C) 2015 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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_hypgeom.h"
+#include "acb_hypgeom.h"
+
+void
+arb_hypgeom_airy(arb_t ai, arb_t aip, arb_t bi, arb_t bip, const arb_t z, slong prec)
+{
+    acb_struct tmp[5];
+
+    if (ai != NULL) acb_init(tmp);
+    if (aip != NULL) acb_init(tmp + 1);
+    if (bi != NULL) acb_init(tmp + 2);
+    if (bip != NULL) acb_init(tmp + 3);
+
+    acb_init(tmp + 4);
+    acb_set_arb(tmp + 4, z);
+
+    acb_hypgeom_airy(ai ? tmp : NULL,
+                     aip ? tmp + 1 : NULL,
+                     bi ? tmp + 2 : NULL,
+                     bip ? tmp + 3 : NULL,
+                     tmp + 4, prec);
+
+    if (ai != NULL) arb_set(ai, acb_realref(tmp));
+    if (aip != NULL) arb_set(aip, acb_realref(tmp + 1));
+    if (bi != NULL) arb_set(bi, acb_realref(tmp + 2));
+    if (bip != NULL) arb_set(bip, acb_realref(tmp + 3));
+
+    if (ai != NULL) acb_clear(tmp);
+    if (aip != NULL) acb_clear(tmp + 1);
+    if (bi != NULL) acb_clear(tmp + 2);
+    if (bip != NULL) acb_clear(tmp + 3);
+
+    acb_clear(tmp + 4);
+}
+

arb_hypgeom/airy_jet.c

@@ -0,0 +1,50 @@
+/*
+    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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_hypgeom.h"
+
+static void
+airy_recurrence(arb_ptr ai, const arb_t z, slong len, slong prec)
+{
+    slong k;
+
+    if (len >= 3)
+    {
+        arb_mul(ai + 2, ai, z, prec);
+        arb_mul_2exp_si(ai + 2, ai + 2, -1);
+    }
+
+    for (k = 3; k < len; k++)
+    {
+        arb_mul(ai + k, ai + k - 2, z, prec);
+        arb_add(ai + k, ai + k, ai + k - 3, prec);
+        arb_div_ui(ai + k, ai + k, k * (k - 1), prec);
+    }
+}
+
+void
+arb_hypgeom_airy_jet(arb_ptr ai, arb_ptr bi, const arb_t z, slong len, slong prec)
+{
+    if (len <= 0)
+        return;
+
+    if (len == 1)
+    {
+        arb_hypgeom_airy(ai, NULL, bi, NULL, z, prec);
+        return;
+    }
+
+    arb_hypgeom_airy(ai, ai ? (ai + 1) : NULL, bi, bi ? (bi + 1) : NULL, z, prec);
+
+    if (ai != NULL) airy_recurrence(ai, z, len, prec);
+    if (bi != NULL) airy_recurrence(bi, z, len, prec);
+}
+

arb_hypgeom/airy_series.c

@@ -0,0 +1,127 @@
+/*
+    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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_hypgeom.h"
+
+void
+_arb_hypgeom_airy_series(arb_ptr ai, arb_ptr ai_prime,
+    arb_ptr bi, arb_ptr bi_prime, arb_srcptr z, slong zlen, slong len, slong prec)
+{
+    arb_ptr t, u, v;
+    slong tlen = len + ((ai_prime != NULL) || (bi_prime != NULL));
+
+    zlen = FLINT_MIN(zlen, len);
+
+    if (zlen <= 0)
+        return;
+
+    if (zlen == 1)
+    {
+        arb_hypgeom_airy(ai, ai_prime, bi, bi_prime, z, prec);
+        return;
+    }
+
+    t = _arb_vec_init(tlen);
+    u = _arb_vec_init(tlen);
+    v = _arb_vec_init(len);
+
+    arb_hypgeom_airy_jet((ai || ai_prime) ? t : NULL,
+                         (bi || bi_prime) ? u : NULL, z, tlen, prec);
+
+    /* compose with nonconstant part */
+    arb_zero(v);
+    _arb_vec_set(v + 1, z + 1, zlen - 1);
+
+    if (ai != NULL) _arb_poly_compose_series(ai, t, len, v, zlen, len, prec);
+    if (bi != NULL) _arb_poly_compose_series(bi, u, len, v, zlen, len, prec);
+
+    /* todo: use chain rule to avoid compositions for derivatives? */
+    if (ai_prime != NULL)
+    {
+        _arb_poly_derivative(t, t, len + 1, prec);
+        _arb_poly_compose_series(ai_prime, t, len, v, zlen, len, prec);
+    }
+
+    if (bi_prime != NULL)
+    {
+        _arb_poly_derivative(u, u, len + 1, prec);
+        _arb_poly_compose_series(bi_prime, u, len, v, zlen, len, prec);
+    }
+
+    _arb_vec_clear(t, tlen);
+    _arb_vec_clear(u, tlen);
+    _arb_vec_clear(v, len);
+}
+
+void
+arb_hypgeom_airy_series(arb_poly_t ai, arb_poly_t ai_prime,
+    arb_poly_t bi, arb_poly_t bi_prime, const arb_poly_t z, slong len, slong prec)
+{
+    if (len == 0)
+    {
+        if (ai != NULL) arb_poly_zero(ai);
+        if (ai_prime != NULL) arb_poly_zero(ai_prime);
+        if (bi != NULL) arb_poly_zero(bi);
+        if (bi_prime != NULL) arb_poly_zero(bi_prime);
+        return;
+    }
+
+    if (z->length <= 1)
+        len = 1;
+
+    if (ai != NULL) arb_poly_fit_length(ai, len);
+    if (ai_prime != NULL) arb_poly_fit_length(ai_prime, len);
+    if (bi != NULL) arb_poly_fit_length(bi, len);
+    if (bi_prime != NULL) arb_poly_fit_length(bi_prime, len);
+
+    if (z->length == 0)
+    {
+        arb_t t;
+        arb_init(t);
+        _arb_hypgeom_airy_series(
+            ai ? ai->coeffs : NULL, ai_prime ? ai_prime->coeffs : NULL,
+            bi ? bi->coeffs : NULL, bi_prime ? bi_prime->coeffs : NULL,
+            t, 1, len, prec);
+        arb_clear(t);
+    }
+    else
+    {
+        _arb_hypgeom_airy_series(
+            ai ? ai->coeffs : NULL, ai_prime ? ai_prime->coeffs : NULL,
+            bi ? bi->coeffs : NULL, bi_prime ? bi_prime->coeffs : NULL,
+            z->coeffs, z->length, len, prec);
+    }
+
+    if (ai != NULL)
+    {
+        _arb_poly_set_length(ai, len);
+        _arb_poly_normalise(ai);
+    }
+
+    if (ai_prime != NULL)
+    {
+        _arb_poly_set_length(ai_prime, len);
+        _arb_poly_normalise(ai_prime);
+    }
+
+    if (bi != NULL)
+    {
+        _arb_poly_set_length(bi, len);
+        _arb_poly_normalise(bi);
+    }
+
+    if (bi_prime != NULL)
+    {
+        _arb_poly_set_length(bi_prime, len);
+        _arb_poly_normalise(bi_prime);
+    }
+}
+

doc/source/arb_hypgeom.rst

@@ -191,3 +191,31 @@ Exponential and trigonometric integrals
     Computes the logarithmic integral (optionally the offset version)
     of the power series *z*, truncated to length *len*.
 
+Airy functions
+-------------------------------------------------------------------------------
+
+.. function:: void arb_hypgeom_airy(arb_t ai, arb_t ai_prime, arb_t bi, arb_t bi_prime, const arb_t z, slong prec)
+
+    Computes the Airy functions `(\operatorname{Ai}(z), \operatorname{Ai}'(z), \operatorname{Bi}(z), \operatorname{Bi}'(z))`
+    simultaneously. Any of the four function values can be omitted by passing
+    *NULL* for the unwanted output variables, speeding up the evaluation.
+
+.. function:: void arb_hypgeom_airy_jet(arb_ptr ai, arb_ptr bi, const arb_t z, slong len, slong prec)
+
+    Writes to *ai* and *bi* the respective Taylor expansions of the Airy functions
+    at the point *z*, truncated to length *len*.
+    Either of the outputs can be *NULL* to avoid computing that function.
+    The variable *z* is not allowed to be aliased with the outputs.
+    To simplify the implementation, this method does not compute the
+    series expansions of the primed versions directly; these are
+    easily obtained by computing one extra coefficient and differentiating
+    the output with :func:`_arb_poly_derivative`.
+
+.. function:: void _arb_hypgeom_airy_series(arb_ptr ai, arb_ptr ai_prime, arb_ptr bi, arb_ptr bi_prime, arb_srcptr z, slong zlen, slong len, slong prec)
+
+.. function:: void arb_hypgeom_airy_series(arb_poly_t ai, arb_poly_t ai_prime, arb_poly_t bi, arb_poly_t bi_prime, const arb_poly_t z, slong len, slong prec)
+
+    Computes the Airy functions evaluated at the power series *z*,
+    truncated to length *len*. As with the other Airy methods, any of the
+    outputs can be *NULL*.
+