From fc4c28f72cab0df2377959e896bcfca2380037c6 Mon Sep 17 00:00:00 2001 From: fredrik Date: Sun, 24 Jun 2018 17:29:24 +0200 Subject: [PATCH] more arb_mat special matrices: pascal, stirling --- arb_mat.h | 4 ++ arb_mat/pascal.c | 76 +++++++++++++++++++++++++ arb_mat/stirling.c | 117 ++++++++++++++++++++++++++++++++++++++ arb_mat/test/t-pascal.c | 88 ++++++++++++++++++++++++++++ arb_mat/test/t-stirling.c | 71 +++++++++++++++++++++++ doc/source/arb_mat.rst | 28 +++++++++ 6 files changed, 384 insertions(+) create mode 100644 arb_mat/pascal.c create mode 100644 arb_mat/stirling.c create mode 100644 arb_mat/test/t-pascal.c create mode 100644 arb_mat/test/t-stirling.c diff --git a/arb_mat.h b/arb_mat.h index 71093676..5616f285 100644 --- a/arb_mat.h +++ b/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); diff --git a/arb_mat/pascal.c b/arb_mat/pascal.c new file mode 100644 index 00000000..23422309 --- /dev/null +++ b/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 . +*/ + +#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); + } +} + diff --git a/arb_mat/stirling.c b/arb_mat/stirling.c new file mode 100644 index 00000000..9c5af176 --- /dev/null +++ b/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 . +*/ + +#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); +} + diff --git a/arb_mat/test/t-pascal.c b/arb_mat/test/t-pascal.c new file mode 100644 index 00000000..e7374bfc --- /dev/null +++ b/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 . +*/ + +#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; +} diff --git a/arb_mat/test/t-stirling.c b/arb_mat/test/t-stirling.c new file mode 100644 index 00000000..b2c7a556 --- /dev/null +++ b/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 . +*/ + +#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; +} diff --git a/doc/source/arb_mat.rst b/doc/source/arb_mat.rst index 70d5ebcb..6c988c1c 100644 --- a/doc/source/arb_mat.rst +++ b/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*