add derivative, integral, in-progress series log, exp

arb_poly.h

@@ -110,4 +110,10 @@ void arb_poly_mullow(arb_poly_t z, const arb_poly_t x, const arb_poly_t y, long
 
 void arb_poly_inv_series(arb_poly_t z, const arb_poly_t x, long n);
 
+void arb_poly_derivative(arb_poly_t z, const arb_poly_t x);
+void arb_poly_integral(arb_poly_t z, const arb_poly_t x);
+
+void arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n);
+void arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n);
+
 #endif

arb_poly/derivative.c

@@ -0,0 +1,50 @@
+/*=============================================================================
+
+    This file is part of FLINT.
+
+    FLINT 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.
+
+    FLINT 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 FLINT; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+void
+arb_poly_derivative(arb_poly_t z, const arb_poly_t x)
+{
+    long i, len;
+
+    len = x->length;
+
+    if (len <= 1)
+    {
+        arb_poly_zero(z);
+        return;
+    }
+
+    _arb_poly_fit_length(z, len - 1);
+
+    for (i = 1; i < len; i++)
+        fmpz_mul_ui(arb_poly_coeffs(z) + i - 1, arb_poly_coeffs(x) + i, i);
+
+    z->length = len - 1;
+
+    fmpz_mul_ui(arb_poly_radref(z), arb_poly_radref(x), len - 1);
+    fmpz_set(arb_poly_expref(z), arb_poly_expref(x));
+}

arb_poly/exp_series.c

@@ -0,0 +1,74 @@
+/*=============================================================================
+
+    This file is part of FLINT.
+
+    FLINT 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.
+
+    FLINT 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 FLINT; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+void
+arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n)
+{
+    long a[FLINT_BITS];
+    long i, m;
+    arb_poly_t t;
+    long xlen;
+
+    xlen = x->length;
+
+    if (z == x)
+    {
+        arb_poly_t t;
+        arb_poly_init(t, arb_poly_prec(x));
+        arb_poly_set(t, x);
+        arb_poly_exp_series(z, t, n);
+        arb_poly_clear(t);
+        return;
+    }
+
+    a[i = 0] = n;
+    while (n >= 2)
+        a[++i] = (n = (n + 1) / 2);
+
+    arb_poly_init(t, arb_poly_prec(z));
+
+    /* first coefficient (assuming input zero, todo: general case) */
+    arb_poly_set_si(z, 1);
+    fmpz_mul_2exp(arb_poly_coeffs(z), arb_poly_coeffs(z), arb_poly_prec(z));
+    fmpz_set_si(arb_poly_expref(z), -arb_poly_prec(z));
+    fmpz_zero(arb_poly_radref(z));
+
+    /* Newton iteration */
+    for (i--; i >= 0; i--)
+    {
+        m = n;
+        n = a[i];
+
+        arb_poly_log_series(t, z, n);
+        arb_poly_neg(t, t);
+        arb_poly_add(t, t, x);
+        arb_poly_mullow(t, z, t, n);
+        arb_poly_add(z, z, t);
+    }
+
+    arb_poly_clear(t);
+}

arb_poly/integral.c

@@ -0,0 +1,72 @@
+/*=============================================================================
+
+    This file is part of FLINT.
+
+    FLINT 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.
+
+    FLINT 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 FLINT; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+void
+arb_poly_integral(arb_poly_t z, const arb_poly_t x)
+{
+    long i, len, bits, shift;
+
+    len = x->length;
+
+    if (len == 0)
+    {
+        arb_poly_zero(z);
+        return;
+    }
+
+    _arb_poly_fit_length(z, len + 1);
+
+    bits = _fmpz_vec_max_bits(arb_poly_coeffs(x), len);
+    bits = FLINT_ABS(bits);
+    shift = arb_poly_prec(z) - bits + FLINT_BIT_COUNT(len);
+
+    if (shift >= 0)
+    {
+        for (i = len; i > 0; i--)
+        {
+            fmpz_mul_2exp(arb_poly_coeffs(z) + i, arb_poly_coeffs(x) + i - 1, shift);
+            fmpz_tdiv_q_ui(arb_poly_coeffs(z) + i, arb_poly_coeffs(z) + i, i);
+        }
+
+        fmpz_mul_2exp(arb_poly_radref(z), arb_poly_radref(x), shift);
+        fmpz_add_ui(arb_poly_radref(z), arb_poly_radref(z), 1);
+        fmpz_sub_ui(arb_poly_expref(z), arb_poly_expref(x), shift);
+    }
+    else
+    {
+        /* todo: shift down to match target precision? */
+        for (i = len; i > 0; i--)
+            fmpz_tdiv_q_ui(arb_poly_coeffs(z) + i, arb_poly_coeffs(x) + i - 1, i);
+
+        fmpz_add_ui(arb_poly_radref(z), arb_poly_radref(x), 1);
+        fmpz_set(arb_poly_expref(z), arb_poly_expref(x));
+    }
+
+    fmpz_zero(arb_poly_coeffs(z));
+
+    z->length = len + 1;
+}

arb_poly/log_series.c

@@ -0,0 +1,43 @@
+/*=============================================================================
+
+    This file is part of FLINT.
+
+    FLINT 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.
+
+    FLINT 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 FLINT; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+void
+arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n)
+{
+    arb_poly_t t, u;
+
+    arb_poly_init(t, arb_poly_prec(z));
+    arb_poly_init(u, arb_poly_prec(z));
+
+    arb_poly_inv_series(t, x, n);
+    arb_poly_derivative(u, x);
+    arb_poly_mullow(t, t, u, n - 1);
+    arb_poly_integral(z, t);
+
+    arb_poly_clear(t);
+    arb_poly_clear(u);
+}

arb_poly/test/t-derivative.c

@@ -0,0 +1,105 @@
+/*=============================================================================
+
+    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) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("derivative....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arb_poly_t a, b;
+        fmpq_poly_t x, y;
+
+        arb_poly_init(a, 1 + n_randint(state, 200));
+        arb_poly_init(b, 1 + n_randint(state, 200));
+
+        fmpq_poly_init(x);
+        fmpq_poly_init(y);
+
+        arb_poly_randtest(a, state, n_randint(state, 20), 10);
+        arb_poly_randtest(b, state, n_randint(state, 20), 10);
+
+        arb_poly_get_rand_fmpq_poly(x, state, a);
+
+        arb_poly_derivative(b, a);
+        fmpq_poly_derivative(y, x);
+
+        if (!arb_poly_contains_fmpq_poly(b, y))
+        {
+            printf("FAIL: containment\n\n");
+            printf("a = "); arb_poly_debug(a); printf("\n\n");
+            printf("x = "); fmpq_poly_print(x); printf("\n\n");
+            printf("b = "); arb_poly_debug(b); printf("\n\n");
+            printf("y = "); fmpq_poly_print(y); printf("\n\n");
+            abort();
+        }
+
+        arb_poly_clear(a);
+        arb_poly_clear(b);
+
+        fmpq_poly_clear(x);
+        fmpq_poly_clear(y);
+    }
+
+    /* aliasing */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arb_poly_t a;
+        fmpq_poly_t x;
+
+        arb_poly_init(a, 1 + n_randint(state, 200));
+        fmpq_poly_init(x);
+
+        arb_poly_randtest(a, state, n_randint(state, 20), 10);
+        arb_poly_get_rand_fmpq_poly(x, state, a);
+
+        arb_poly_derivative(a, a);
+        fmpq_poly_derivative(x, x);
+
+        if (!arb_poly_contains_fmpq_poly(a, x))
+        {
+            printf("FAIL: aliasing\n\n");
+            printf("a = "); arb_poly_debug(a); printf("\n\n");
+            printf("x = "); fmpq_poly_print(x); printf("\n\n");
+            abort();
+        }
+
+        arb_poly_clear(a);
+        fmpq_poly_clear(x);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb_poly/test/t-integral.c

@@ -0,0 +1,105 @@
+/*=============================================================================
+
+    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) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb_poly.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("integral....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arb_poly_t a, b;
+        fmpq_poly_t x, y;
+
+        arb_poly_init(a, 1 + n_randint(state, 200));
+        arb_poly_init(b, 1 + n_randint(state, 200));
+
+        fmpq_poly_init(x);
+        fmpq_poly_init(y);
+
+        arb_poly_randtest(a, state, n_randint(state, 20), 10);
+        arb_poly_randtest(b, state, n_randint(state, 20), 10);
+
+        arb_poly_get_rand_fmpq_poly(x, state, a);
+
+        arb_poly_integral(b, a);
+        fmpq_poly_integral(y, x);
+
+        if (!arb_poly_contains_fmpq_poly(b, y))
+        {
+            printf("FAIL: containment\n\n");
+            printf("a = "); arb_poly_debug(a); printf("\n\n");
+            printf("x = "); fmpq_poly_print(x); printf("\n\n");
+            printf("b = "); arb_poly_debug(b); printf("\n\n");
+            printf("y = "); fmpq_poly_print(y); printf("\n\n");
+            abort();
+        }
+
+        arb_poly_clear(a);
+        arb_poly_clear(b);
+
+        fmpq_poly_clear(x);
+        fmpq_poly_clear(y);
+    }
+
+    /* aliasing */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        arb_poly_t a;
+        fmpq_poly_t x;
+
+        arb_poly_init(a, 1 + n_randint(state, 200));
+        fmpq_poly_init(x);
+
+        arb_poly_randtest(a, state, n_randint(state, 20), 10);
+        arb_poly_get_rand_fmpq_poly(x, state, a);
+
+        arb_poly_integral(a, a);
+        fmpq_poly_integral(x, x);
+
+        if (!arb_poly_contains_fmpq_poly(a, x))
+        {
+            printf("FAIL: aliasing\n\n");
+            printf("a = "); arb_poly_debug(a); printf("\n\n");
+            printf("x = "); fmpq_poly_print(x); printf("\n\n");
+            abort();
+        }
+
+        arb_poly_clear(a);
+        fmpq_poly_clear(x);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}