MAINT: is_square and is_empty for arb_mat and acb_mat

acb_mat.h

@@ -126,6 +126,18 @@ int acb_mat_contains_fmpz_mat(const acb_mat_t mat1, const fmpz_mat_t mat2);
 
 int acb_mat_is_real(const acb_mat_t mat);
 
+ACB_MAT_INLINE int
+acb_mat_is_empty(const acb_mat_t mat)
+{
+    return (mat->r == 0) || (mat->c == 0);
+}
+
+ACB_MAT_INLINE int
+acb_mat_is_square(const acb_mat_t mat)
+{
+    return (mat->r == mat->c);
+}
+
 /* Special matrices */
 
 void acb_mat_zero(acb_mat_t mat);

acb_mat/det.c

@@ -153,7 +153,15 @@ acb_mat_det_inplace(acb_t det, acb_mat_t A, slong prec)
 void
 acb_mat_det(acb_t det, const acb_mat_t A, slong prec)
 {
-    slong n = acb_mat_nrows(A);
+    slong n;
+
+    if (!acb_mat_is_square(A))
+    {
+        flint_printf("acb_mat_det: a square matrix is required!\n");
+        abort();
+    }
+
+    n = acb_mat_nrows(A);
 
     if (n == 0)
     {

acb_mat/exp.c

@@ -193,19 +193,18 @@ acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
     acb_mat_t T;
     int is_real;
 
-    dim = acb_mat_nrows(A);
+    if (!acb_mat_is_square(A))
-
-    if (dim != acb_mat_ncols(A))
     {
         flint_printf("acb_mat_exp: a square matrix is required!\n");
         abort();
     }
 
-    if (dim == 0)
+    if (acb_mat_is_empty(A))
-    {
         return;
-    }
+
-    else if (dim == 1)
+    dim = acb_mat_nrows(A);
+
+    if (dim == 1)
     {
         acb_exp(acb_mat_entry(B, 0, 0), acb_mat_entry(A, 0, 0), prec);
         return;

acb_mat/lu.c

@@ -33,14 +33,12 @@ acb_mat_lu(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec)
     slong i, j, m, n, r, row, col;
     int result;
 
+    if (acb_mat_is_empty(A))
+        return 1;
+
     m = acb_mat_nrows(A);
     n = acb_mat_ncols(A);
 
-    result = 1;
-
-    if (m == 0 || n == 0)
-        return result;
-
     acb_mat_set(LU, A);
 
     a = LU->rows;
@@ -52,6 +50,8 @@ acb_mat_lu(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec)
     acb_init(d);
     acb_init(e);
 
+    result = 1;
+
     while (row < m && col < n)
     {
         r = acb_mat_find_pivot_partial(LU, row, m, col);

acb_mat/sqr.c

@@ -40,10 +40,7 @@ acb_mat_sqr(acb_mat_t B, const acb_mat_t A, slong prec)
     }
 
     if (n == 0)
-    {
-        acb_mat_zero(B);
         return;
-    }
 
     if (n == 1)
     {

acb_mat/trace.c

@@ -28,18 +28,23 @@
 void
 acb_mat_trace(acb_t trace, const acb_mat_t mat, slong prec)
 {
-    slong i, n = acb_mat_nrows(mat);
+    slong i;
 
-    if (n == 0)
+    if (!acb_mat_is_square(mat))
+    {
+        flint_printf("acb_mat_trace: a square matrix is required!\n");
+        abort();
+    }
+
+    if (acb_mat_is_empty(mat))
     {
         acb_zero(trace);
+        return;
     }
-    else
+
+    acb_set(trace, acb_mat_entry(mat, 0, 0));
+    for (i = 1; i < acb_mat_nrows(mat); i++)
     {
-        acb_set(trace, acb_mat_entry(mat, 0, 0));
+        acb_add(trace, trace, acb_mat_entry(mat, i, i), prec);
-        for (i = 1; i < n; i++)
-        {
-            acb_add(trace, trace, acb_mat_entry(mat, i, i), prec);
-        }
     }
 }

acb_mat/transpose.c

@@ -28,7 +28,6 @@
 void
 acb_mat_transpose(acb_mat_t B, const acb_mat_t A)
 {
-    acb_struct tmp;
     slong i, j;
 
     if (acb_mat_nrows(B) != acb_mat_ncols(A) || acb_mat_ncols(B) != acb_mat_nrows(A))
@@ -37,15 +36,16 @@ acb_mat_transpose(acb_mat_t B, const acb_mat_t A)
         abort();
     }
 
+    if (acb_mat_is_empty(A))
+        return;
+
     if (A == B)  /* In-place, guaranteed to be square */
     {
         for (i = 0; i < acb_mat_nrows(A) - 1; i++)
         {
             for (j = i + 1; j < acb_mat_ncols(A); j++)
             {
-                tmp = *acb_mat_entry(A, i, j);
+                acb_swap(acb_mat_entry(A, i, j), acb_mat_entry(A, j, i));
-                *acb_mat_entry(A, i, j) = *acb_mat_entry(A, j, i);
-                *acb_mat_entry(A, j, i) = tmp;
             }
         }
     }

arb_mat.h

@@ -118,6 +118,18 @@ int arb_mat_contains_fmpq_mat(const arb_mat_t mat1, const fmpq_mat_t mat2);
 
 int arb_mat_contains_fmpz_mat(const arb_mat_t mat1, const fmpz_mat_t mat2);
 
+ARB_MAT_INLINE int
+arb_mat_is_empty(const arb_mat_t mat)
+{
+    return (mat->r == 0) || (mat->c == 0);
+}
+
+ARB_MAT_INLINE int
+arb_mat_is_square(const arb_mat_t mat)
+{
+    return (mat->r == mat->c);
+}
+
 /* Special matrices */
 
 void arb_mat_zero(arb_mat_t mat);

arb_mat/det.c

@@ -137,7 +137,15 @@ arb_mat_det_inplace(arb_t det, arb_mat_t A, slong prec)
 void
 arb_mat_det(arb_t det, const arb_mat_t A, slong prec)
 {
-    slong n = arb_mat_nrows(A);
+    slong n;
+
+    if (!arb_mat_is_square(A))
+    {
+        flint_printf("arb_mat_det: a square matrix is required!\n");
+        abort();
+    }
+
+    n = arb_mat_nrows(A);
 
     if (n == 0)
     {

arb_mat/exp.c

@@ -186,19 +186,18 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
     mag_t norm, err;
     arb_mat_t T;
 
-    dim = arb_mat_nrows(A);
+    if (!arb_mat_is_square(A))
-
-    if (dim != arb_mat_ncols(A))
     {
         flint_printf("arb_mat_exp: a square matrix is required!\n");
         abort();
     }
 
-    if (dim == 0)
+    if (arb_mat_is_empty(A))
-    {
         return;
-    }
+
-    else if (dim == 1)
+    dim = arb_mat_nrows(A);
+
+    if (dim == 1)
     {
         arb_exp(arb_mat_entry(B, 0, 0), arb_mat_entry(A, 0, 0), prec);
         return;

arb_mat/lu.c

@@ -33,14 +33,12 @@ arb_mat_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec)
     slong i, j, m, n, r, row, col;
     int result;
 
+    if (arb_mat_is_empty(A))
+        return 1;
+
     m = arb_mat_nrows(A);
     n = arb_mat_ncols(A);
 
-    result = 1;
-
-    if (m == 0 || n == 0)
-        return result;
-
     arb_mat_set(LU, A);
 
     a = LU->rows;
@@ -52,6 +50,8 @@ arb_mat_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec)
     arb_init(d);
     arb_init(e);
 
+    result = 1;
+
     while (row < m && col < n)
     {
         r = arb_mat_find_pivot_partial(LU, row, m, col);

arb_mat/sqr_classical.c

@@ -40,10 +40,7 @@ arb_mat_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec)
     }
 
     if (n == 0)
-    {
-        arb_mat_zero(B);
         return;
-    }
 
     if (n == 1)
     {

arb_mat/trace.c

@@ -28,18 +28,23 @@
 void
 arb_mat_trace(arb_t trace, const arb_mat_t mat, slong prec)
 {
-    slong i, n = arb_mat_nrows(mat);
+    slong i;
 
-    if (n == 0)
+    if (!arb_mat_is_square(mat))
+    {
+        flint_printf("arb_mat_trace: a square matrix is required!\n");
+        abort();
+    }
+
+    if (arb_mat_is_empty(mat))
     {
         arb_zero(trace);
+        return;
     }
-    else
+
+    arb_set(trace, arb_mat_entry(mat, 0, 0));
+    for (i = 1; i < arb_mat_nrows(mat); i++)
     {
-        arb_set(trace, arb_mat_entry(mat, 0, 0));
+        arb_add(trace, trace, arb_mat_entry(mat, i, i), prec);
-        for (i = 1; i < n; i++)
-        {
-            arb_add(trace, trace, arb_mat_entry(mat, i, i), prec);
-        }
     }
 }

arb_mat/transpose.c

@@ -28,7 +28,6 @@
 void
 arb_mat_transpose(arb_mat_t B, const arb_mat_t A)
 {
-    arb_struct tmp;
     slong i, j;
 
     if (arb_mat_nrows(B) != arb_mat_ncols(A) || arb_mat_ncols(B) != arb_mat_nrows(A))
@@ -37,15 +36,16 @@ arb_mat_transpose(arb_mat_t B, const arb_mat_t A)
         abort();
     }
 
+    if (arb_mat_is_empty(A))
+        return;
+
     if (A == B)  /* In-place, guaranteed to be square */
     {
         for (i = 0; i < arb_mat_nrows(A) - 1; i++)
         {
             for (j = i + 1; j < arb_mat_ncols(A); j++)
             {
-                tmp = *arb_mat_entry(A, i, j);
+                arb_swap(arb_mat_entry(A, i, j), arb_mat_entry(A, j, i));
-                *arb_mat_entry(A, i, j) = *arb_mat_entry(A, j, i);
-                *arb_mat_entry(A, j, i) = tmp;
             }
         }
     }
@@ -56,4 +56,3 @@ arb_mat_transpose(arb_mat_t B, const arb_mat_t A)
                 arb_set(arb_mat_entry(B, i, j), arb_mat_entry(A, j, i));
     }
 }
-

doc/source/acb_mat.rst

@@ -124,6 +124,15 @@ Comparisons
 
     Returns nonzero iff all entries in *mat* have zero imaginary part.
 
+.. function:: int acb_mat_is_empty(const acb_mat_t mat)
+
+    Returns nonzero iff the number of rows or the number of columns in *mat* is zero.
+
+.. function:: int acb_mat_is_square(const acb_mat_t mat)
+
+    Returns nonzero iff the number of rows is equal to the number of columns in *mat*.
+
+
 Special matrices
 -------------------------------------------------------------------------------
 

doc/source/arb_mat.rst

@@ -116,6 +116,14 @@ Comparisons
 
     Returns nonzero iff *mat1* and *mat2* certainly do not represent the same matrix.
 
+.. function:: int arb_mat_is_empty(const arb_mat_t mat)
+
+    Returns nonzero iff the number of rows or the number of columns in *mat* is zero.
+
+.. function:: int arb_mat_is_square(const arb_mat_t mat)
+
+    Returns nonzero iff the number of rows is equal to the number of columns in *mat*.
+
 Special matrices
 -------------------------------------------------------------------------------