more arb_mat special matrices: pascal, stirling

arb_mat.h

@@ -148,6 +148,10 @@ void arb_mat_ones(arb_mat_t mat);
 
 void arb_mat_hilbert(arb_mat_t mat, slong prec);
 
+void arb_mat_pascal(arb_mat_t mat, int triangular, slong prec);
+
+void arb_mat_stirling(arb_mat_t mat, int kind, slong prec);
+
 void arb_mat_dct(arb_mat_t mat, int type, slong prec);
 
 void arb_mat_transpose(arb_mat_t mat1, const arb_mat_t mat2);

arb_mat/pascal.c

@@ -0,0 +1,76 @@
+/*
+    Copyright (C) 2018 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_mat.h"
+
+void
+arb_mat_pascal(arb_mat_t mat, int triangular, slong prec)
+{
+    slong R, C, i, j;
+
+    R = arb_mat_nrows(mat);
+    C = arb_mat_ncols(mat);
+
+    if (R == 0 || C == 0)
+        return;
+
+    if (triangular == 1)
+    {
+        for (i = 1; i < R; i++)
+            for (j = 0; j < i && j < C; j++)
+                arb_zero(arb_mat_entry(mat, i, j));
+
+        for (j = 0; j < C; j++)
+            arb_one(arb_mat_entry(mat, 0, j));
+
+        for (i = 1; i < R && i < C; i++)
+            arb_one(arb_mat_entry(mat, i, i));
+
+        for (i = 1; i < R; i++)
+            for (j = i + 1; j < C; j++)
+                arb_add(arb_mat_entry(mat, i, j),
+                    arb_mat_entry(mat, i, j - 1),
+                    arb_mat_entry(mat, i - 1, j - 1), prec);
+    }
+    else if (triangular == -1)
+    {
+        for (i = 0; i < R; i++)
+            for (j = i + 1; j < C; j++)
+                arb_zero(arb_mat_entry(mat, i, j));
+
+        for (i = 0; i < R; i++)
+            arb_one(arb_mat_entry(mat, i, 0));
+
+        for (i = 1; i < R && i < C; i++)
+            arb_one(arb_mat_entry(mat, i, i));
+
+        for (i = 2; i < R; i++)
+            for (j = 1; j < i && j < C; j++)
+                arb_add(arb_mat_entry(mat, i, j),
+                    arb_mat_entry(mat, i - 1, j - 1),
+                    arb_mat_entry(mat, i - 1, j), prec);
+    }
+    else
+    {
+        for (j = 0; j < C; j++)
+            arb_one(arb_mat_entry(mat, 0, j));
+
+        for (i = 1; i < R; i++)
+            arb_one(arb_mat_entry(mat, i, 0));
+
+        for (i = 1; i < R; i++)
+            for (j = 1; j < C; j++)
+                arb_add(arb_mat_entry(mat, i, j),
+                    arb_mat_entry(mat, i, j - 1),
+                    arb_mat_entry(mat, i - 1, j), prec);
+    }
+}
+

arb_mat/stirling.c

@@ -0,0 +1,117 @@
+/*
+    Copyright (C) 2018 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_mat.h"
+
+static void _stirling_number_1u_vec_next(arb_ptr row,
+    arb_srcptr prev, slong n, slong klen, slong prec)
+{
+    slong k;
+
+    if (klen > n) arb_one(row + n);
+    if (n != 0 && klen != 0) arb_zero(row);
+
+    for (k = FLINT_MIN(n, klen) - 1; k >= 1; k--)
+    {
+        arb_mul_ui(row + k, prev + k, n - 1, prec);
+        arb_add(row + k, prev + k - 1, row + k, prec);
+    }
+
+    for (k = n + 1; k < klen; k++)
+        arb_zero(row + k);
+}
+
+static void _stirling_number_1_vec_next(arb_ptr row,
+    arb_srcptr prev, slong n, slong klen, slong prec)
+{
+    slong k;
+
+    if (klen > n) arb_one(row + n);
+    if (n != 0 && klen != 0) arb_zero(row);
+
+    for (k = FLINT_MIN(n, klen) - 1; k >= 1; k--)
+    {
+        arb_mul_ui(row + k, prev + k, n - 1, prec);
+        arb_sub(row + k, prev + k - 1, row + k, prec);
+    }
+
+    for (k = n + 1; k < klen; k++)
+        arb_zero(row + k);
+}
+
+static void _stirling_number_2_vec_next(arb_ptr row,
+    arb_srcptr prev, slong n, slong klen, slong prec)
+{
+    slong k;
+
+    if (klen > n) arb_one(row + n);
+    if (n != 0 && klen != 0) arb_zero(row);
+
+    for (k = FLINT_MIN(n, klen) - 1; k >= 1; k--)
+    {
+        arb_mul_ui(row + k, prev + k, k, prec);
+        arb_add(row + k, prev + k - 1, row + k, prec);
+    }
+
+    for (k = n + 1; k < klen; k++)
+        arb_zero(row + k);
+}
+
+static void
+_stirling_matrix_1u(arb_mat_t mat, slong prec)
+{
+    slong n;
+
+    if (arb_mat_is_empty(mat))
+        return;
+
+    for (n = 0; n < mat->r; n++)
+        _stirling_number_1u_vec_next(mat->rows[n],
+            mat->rows[n - (n != 0)], n, mat->c, prec);
+}
+
+static void
+_stirling_matrix_1(arb_mat_t mat, slong prec)
+{
+    slong n;
+
+    if (arb_mat_is_empty(mat))
+        return;
+
+    for (n = 0; n < mat->r; n++)
+        _stirling_number_1_vec_next(mat->rows[n],
+            mat->rows[n - (n != 0)], n, mat->c, prec);
+}
+
+static void
+_stirling_matrix_2(arb_mat_t mat, slong prec)
+{
+    slong n;
+
+    if (arb_mat_is_empty(mat))
+        return;
+
+    for (n = 0; n < mat->r; n++)
+        _stirling_number_2_vec_next(mat->rows[n],
+            mat->rows[n - (n != 0)], n, mat->c, prec);
+}
+
+void
+arb_mat_stirling(arb_mat_t mat, int kind, slong prec)
+{
+    if (kind == 0)
+        _stirling_matrix_1u(mat, prec);
+    else if (kind == 1)
+        _stirling_matrix_1(mat, prec);
+    else
+        _stirling_matrix_2(mat, prec);
+}
+

arb_mat/test/t-pascal.c

@@ -0,0 +1,88 @@
+/*
+    Copyright (C) 2018 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_mat.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("pascal....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100 * arb_test_multiplier(); iter++)
+    {
+        arb_mat_t A;
+        fmpz_t t;
+        slong n, m, i, j, prec;
+
+        n = n_randint(state, 10);
+        m = n_randint(state, 10);
+        prec = 2 + n_randint(state, 200);
+
+        fmpz_init(t);
+        arb_mat_init(A, n, m);
+        arb_mat_randtest(A, state, 100, 10);
+
+        arb_mat_pascal(A, 0, prec);
+        for (i = 0; i < n; i++)
+        {
+            for (j = 0; j < m; j++)
+            {
+                fmpz_bin_uiui(t, i + j, i);
+                if (!arb_contains_fmpz(arb_mat_entry(A, i, j), t))
+                {
+                    flint_printf("FAIL: containment (0)\n");
+                    flint_abort();
+                }
+            }
+        }
+
+        arb_mat_pascal(A, 1, prec);
+        for (i = 0; i < n; i++)
+        {
+            for (j = 0; j < m; j++)
+            {
+                fmpz_bin_uiui(t, j, i);
+                if (!arb_contains_fmpz(arb_mat_entry(A, i, j), t))
+                {
+                    flint_printf("FAIL: containment (1)\n");
+                    flint_abort();
+                }
+            }
+        }
+
+        arb_mat_pascal(A, -1, prec);
+        for (i = 0; i < n; i++)
+        {
+            for (j = 0; j < m; j++)
+            {
+                fmpz_bin_uiui(t, i, j);
+                if (!arb_contains_fmpz(arb_mat_entry(A, i, j), t))
+                {
+                    flint_printf("FAIL: containment (-1)\n");
+                    flint_abort();
+                }
+            }
+        }
+
+        arb_mat_clear(A);
+        fmpz_clear(t);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb_mat/test/t-stirling.c

@@ -0,0 +1,71 @@
+/*
+    Copyright (C) 2018 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_mat.h"
+#include "flint/arith.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("stirling....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100 * arb_test_multiplier(); iter++)
+    {
+        arb_mat_t A;
+        fmpz_mat_t B;
+        slong n, m, prec;
+
+        n = n_randint(state, 10);
+        m = n_randint(state, 10);
+        prec = 2 + n_randint(state, 200);
+
+        arb_mat_init(A, n, m);
+        fmpz_mat_init(B, n, m);
+        arb_mat_randtest(A, state, 100, 10);
+
+        arb_mat_stirling(A, 0, prec);
+        arith_stirling_matrix_1u(B);
+        if (!arb_mat_contains_fmpz_mat(A, B))
+        {
+            flint_printf("FAIL: containment (0)\n");
+            flint_abort();
+        }
+
+        arb_mat_stirling(A, 1, prec);
+        arith_stirling_matrix_1(B);
+        if (!arb_mat_contains_fmpz_mat(A, B))
+        {
+            flint_printf("FAIL: containment (1)\n");
+            flint_abort();
+        }
+
+        arb_mat_stirling(A, 2, prec);
+        arith_stirling_matrix_2(B);
+        if (!arb_mat_contains_fmpz_mat(A, B))
+        {
+            flint_printf("FAIL: containment (2)\n");
+            flint_abort();
+        }
+
+        arb_mat_clear(A);
+        fmpz_mat_clear(B);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/arb_mat.rst

@@ -149,6 +149,34 @@ Special matrices
 
     Sets *mat* to the Hilbert matrix, which has entries `A_{j,k} = 1/(j+k+1)`.
 
+.. function:: void arb_mat_pascal(arb_mat_t mat, int triangular, slong prec)
+
+    Sets *mat* to a Pascal matrix, whose entries are binomial coefficients.
+    If *triangular* is 0, constructs a full symmetric matrix
+    with the rows of Pascal's triangle as successive antidiagonals.
+    If *triangular* is 1, constructs the upper triangular matrix with
+    the rows of Pascal's triangle as columns, and if *triangular* is -1,
+    constructs the lower triangular matrix with the rows of Pascal's
+    triangle as rows.
+
+    The entries are computed using recurrence relations.
+    When the dimensions get large, some precision loss is possible; in that
+    case, the user may wish to create the matrix at slightly higher precision
+    and then round it to the final precision.
+
+.. function:: void arb_mat_stirling(arb_mat_t mat, int kind, slong prec)
+
+    Sets *mat* to a Stirling matrix, whose entries are Stirling numbers.
+    If *kind* is 0, the entries are set to the unsigned Stirling numbers
+    of the first kind. If *kind* is 1, the entries are set to the signed
+    Stirling numbers of the first kind. If *kind* is 2, the entries are
+    set to the Stirling numbers of the second kind.
+
+    The entries are computed using recurrence relations.
+    When the dimensions get large, some precision loss is possible; in that
+    case, the user may wish to create the matrix at slightly higher precision
+    and then round it to the final precision.
+
 .. function:: void arb_mat_dct(arb_mat_t mat, int type, slong prec)
 
     Sets *mat* to the DCT (discrete cosine transform) matrix of order *n*