Compile a matrix argument as an array, with its inverse and whole powers

BREAKING-CHANGES.md

@@ -1767,17 +1767,76 @@ are now also found wherever the `1` and the square sit in a sum, under the condi
 | `"1 + tan(x)^2 - sec(x)^2".Simplify()` | `1 + tan(x) ^ 2 - sec(x) ^ 2` | `0 provided not cos(x) = 0` |
 | `"a + 1 + tan(y)^2 - sec(y)^2".Simplify()` | `a + 1 + tan(y) ^ 2 - sec(y) ^ 2` | `a provided not cos(y) = 0` |
 
-### A scalar over a matrix argument is not compiled
+### A scalar over a matrix argument is the scalar times its inverse
 
-**Not silent: what compiled to a wrong value now throws `UncompilableNodeException`.** A scalar over a matrix that is an
-argument of the compiled function was compiled as the scalar over every element, where evaluation takes it as the scalar
-times the inverse. The inverse is not compiled, as a negative matrix power is not, so the quotient is refused. A matrix
-written into the expression is inverted before it is compiled, as before
-([#2287](https://github.com/asc-community/AngouriMath/issues/2287)).
+**Not silent: what compiled to a wrong value now gives the right one.** A scalar over a matrix that is an argument of the
+compiled function was compiled as the scalar over every element, where evaluation takes it as the scalar times the
+inverse. Since 2.6.0 it was refused while the inverse was not compiled
+([#2287](https://github.com/asc-community/AngouriMath/issues/2287)), and the inverse compiles now: a matrix over a matrix is the first times the inverse of the second, and a negative power
+is that power of the inverse, as evaluation has them. A singular matrix's inverse is NaN in every entry, as a compiled
+division by zero has no value. A matrix of `int`, `long` or `BigInteger` is not inverted: its inverse is not a matrix of
+whole numbers. A matrix written into the expression is inverted before it is compiled, as before.
 
 | Input | Was (2.6.0) | Now |
 |---|---|---|
-| `"x / A"`, compiled over `GenTensor<double>` and called at `A = [[1, 2], [3, 4]]`, `x = 2` | `[[2, 1], [0.667, 0.5]]` | throws `UncompilableNodeException` |
+| `"x / A"`, compiled over `GenTensor<double>` and called at `A = [[1, 2], [3, 4]]`, `x = 2` | `[[2, 1], [0.667, 0.5]]` | `[[-4, 2], [3, -1]]` |
+| `"A ^ (-1)"`, compiled over `GenTensor<double>` and called at `A = [[1, 2], [3, 4]]` | throws `NotSupportedException` when called | `[[-2, 1], [1.5, -0.5]]` |
+| `"A / B"`, compiled over `GenTensor<double>` | throws `UncompilableNodeException` | `A` times the inverse of `B` |
+
+### A compiled matrix power is the matrix power of a whole exponent
+
+**Silent: a matrix to a power that is not whole has no compiled value, where 2.6.0 compiled it entry by entry.** In 2.6.0
+simplification raised each entry to the power before compiling, so `[[x, 1], [0, y]] ^ x` compiled as `[[x ^ x, 1], [0,
+y ^ x]]`. Evaluation takes a power of a matrix as a matrix power since "A matrix is divided, raised to a power and compared
+as a matrix" above, and leaves one that is not whole unevaluated. The compiled delegate now takes the matrix power where the
+exponent is whole, and throws `InvalidMatrixOperationException` when called with one that is not; it no longer casts the
+exponent to `int`, which truncates. An exponent of type `Complex` compiles, and is whole where its imaginary part is zero and
+its real part is whole.
+
+| Input | Was (2.6.0) | Now |
+|---|---|---|
+| `"[[x, 1], [0, y]] ^ x"`, compiled over `double` to `GenTensor<double>` and called at `(0.7, -1.3)` | `[[0.779, 1], [0, NaN]]`, each entry to the power | throws `InvalidMatrixOperationException` |
+| the same compiled over `Complex` to `GenTensor<Complex>` | throws `UncompilableNodeException` | a delegate; the matrix power at a whole `x`, and `InvalidMatrixOperationException` at any other |
+
+### A power of a matrix that is not square is refused when compiling
+
+**Not silent: what compiled to a wrong value now throws `InvalidMatrixOperationException`.** In 2.6.0 simplification
+squared `[[x, 1, y]]` entry by entry, and the delegate returned that. A matrix that is not square has no power, and
+evaluation leaves it written. Simplification also leaves written any operation on written matrices whose shapes do not
+fit, and `Compile` now refuses each of them, naming the shapes, rather than building a delegate that throws whenever it is
+called. A matrix argument's shape is known only when the delegate is called, and is checked then.
+
+| Input | Was (2.6.0) | Now |
+|---|---|---|
+| `"[[x, 1, y]] ^ 2".ToEntity().Compile<Complex, Complex, GenTensor<Complex>>("x", "y")`, called at `(1, 2)` | `[[1, 1, 4]]` | throws `InvalidMatrixOperationException` when compiling |
+
+### A compiled matrix is taken and given as a two-dimensional array, and FastExpression compiles a number made of matrices
+
+**Addition, not silent.** Under the LINQ compiler `double[,]`, `Complex[,]`, `float[,]`, `int[,]`, `long[,]` and
+`BigInteger[,]` are matrices, as arguments and as the result, read into and out of the `GenTensor` the arithmetic runs
+on, so a caller need not know GenericTensor. FastExpression simplifies an expression with a matrix in it before
+compiling it, as the LINQ compiler does since [#425](https://github.com/asc-community/AngouriMath/issues/425), so one whose
+value is a number compiles. One whose value is a matrix is refused as before, and the message names the LINQ form.
+
+| Input | Was (2.6.0) | Now |
+|---|---|---|
+| `"A * B".ToEntity().Compile<double[,], double[,], double[,]>("A", "B")` | throws `UncompilableNodeException` | the matrix product |
+| `"[1, 0]T * [[a, b], [c, d]] * [0, 1]".ToEntity().Compile("a", "b", "c", "d")` | throws `UncompilableNodeException` | `b` |
+
+### A matrix argument is not simplified as a number before compiling
+
+**Silent: what compiled to a wrong value now gives the right one.** The LINQ compiler simplifies an expression with a
+written matrix in it before compiling ([#425](https://github.com/asc-community/AngouriMath/issues/425)), and
+simplification reads a variable as a number. With a matrix argument beside the written matrices, `[[1, 0]] * A * [[0], [1]]`
+was multiplied out to `0`, and the delegate returned that whatever `A` was. Only the parts that do not mention a matrix
+argument are simplified now, a written matrix beside a matrix argument is over the argument's element type, and a matrix
+known when compiling to be 1 x 1 compiles as a number, its entry, as evaluation has it.
+
+| Input | Was (2.6.0) | Now |
+|---|---|---|
+| `"[[1, 0]] * A * [[0], [1]]"`, compiled from `GenTensor<double>` to `double` and called at `A = [[1, 2], [3, -4]]` | `0` | `2` |
+| `"A + [[1, 2], [3, 4]]"`, compiled over `GenTensor<double>`, the same | throws `UncompilableNodeException` | `[[2, 4], [6, 0]]` |
+| `"[[1, 0]] * A"`, compiled over `GenTensor<double>`, the same | throws `KeyNotFoundException` | `[[1, 2]]` |
 
 ### A square system is not answered with rows that are no solution
 

Sources/AngouriMath/Docs/Usage/Exceptions.md

@@ -169,8 +169,11 @@ the wrappers living in `AngouriMath.Core.Compilation.IntoLinq` — and that is t
 to and to pass a matrix in as: `"A * B".ToEntity().Compile<GenTensor<double, DoubleOperations>,
 GenTensor<double, DoubleOperations>, GenTensor<double, DoubleOperations>>("A", "B")` is the matrix
 product, `2 A + B - 1` scales, adds and shifts element by element, `A ^ 3` is the matrix power
-([#526](https://github.com/asc-community/AngouriMath/issues/526)). Compiling a matrix as a
-number, or a number as a matrix, or `sin(A)`, is still this exception.
+([#526](https://github.com/asc-community/AngouriMath/issues/526)). A two-dimensional array of the
+element type -- `double[,]`, `Complex[,]` -- is a matrix too, as an argument and as the result.
+Compiling a matrix as a number, or a number as a matrix, or `sin(A)`, is still this exception, and
+so is dividing by a matrix of `int`, `long` or `BigInteger`, whose inverse is not one of whole
+numbers.
 
 ### `CannotEvalException`
 
@@ -232,6 +235,10 @@ the matrix is not square, `AsScalar` on something that is not 1x1, a tensor of t
 dimensions, a row of the wrong width added to a `MatrixBuilder`. The messages name the shapes
 involved.
 
+A compiled matrix raises it as well: `Compile` where the shapes of matrices written into the
+expression do not fit, and the compiled delegate where the shapes of matrix arguments do not, or
+where a matrix power's exponent is not a whole number.
+
 ### `BadMatrixShapeException`
 
 The same subject, narrower: matrices being concatenated whose extents disagree along the joining

Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs

@@ -97,7 +97,9 @@ private protected override void CompileNode(Compiler compiler) =>
         partial record Matrix
         {
             private protected override void CompileNode(Compiler compiler) =>
-                throw new UncompilableNodeException($"A {nameof(Matrix)} cannot be compiled: `{Stringize()}`");
+                throw new UncompilableNodeException(
+                    $"A {nameof(Matrix)} cannot be compiled to a FastExpression, whose value is a number: `{Stringize()}`. "
+                    + "Compile it with Compile<..., double[,]> or Compile<..., Complex[,]> to have the matrix");
         }
 
         // Each function and operator processing
@@ -1020,6 +1022,11 @@ internal static FastExpression Compile(Entity func, IEnumerable<Variable> variab
                         varNamespace[varName] = id++;
                 foreach (var constant in Variable.NamedConstants.Values)
                     func = func.Substitute(constant, constant.Value);
+                // An expression built out of matrices often has a number for its value -- [1, 0]T * [[a, b], [c, d]] * [0, 1]
+                // is b -- and the LINQ compiler simplifies such an expression first, for that reason. This one did not, and
+                // refused it. https://github.com/asc-community/AngouriMath/issues/425
+                if (func.Nodes.Any(node => node is Matrix))
+                    func = func.InnerSimplified;
                 func = Compilation.CompiledMembership.Rewritten(func);
                 func = Functions.ExponentialQuotients.Rewritten(func);
                 Compilation.StatementPositions.Check(func);

Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs

@@ -93,9 +93,12 @@ public virtual Expression ConvertType(Expression expr, Type type)
                     throw new UncompilableNodeException($"A statement that may have no truth value cannot be compiled as {type}, which has no NaN");
                 return ConvertType(Expression.Call(ThreeValuedLogic.IversonMethod, ConvertType(expr, typeof(bool?))), type);
             }
+            // A tensor is given back as the two-dimensional array of its element type where that is asked for.
+            if (TensorCompilation.AsArray(expr, type) is { } array)
+                return array;
             // A tensor converts to nothing but itself: a matrix over one element type is not
             // a matrix over another, and it is not a number.
-            if (TensorCompilation.IsTensor(expr.Type) || TensorCompilation.IsTensor(type))
+            if (TensorCompilation.IsMatrixType(expr.Type) || TensorCompilation.IsMatrixType(type))
                 throw new UncompilableNodeException($"A value of type {expr.Type} cannot be compiled as {type}");
 
             bool exprNullable = (Nullable.GetUnderlyingType(expr.Type) is not null) || (expr.Type == typeof(object));

Sources/AngouriMath/Functions/Compilation/IntoLinq/IntoLinqCompiler.cs

@@ -22,17 +22,21 @@ internal static TDelegate Compile<TDelegate>(
             IEnumerable<(Type type, Variable variable)> typesAndNames
             ) where TDelegate : Delegate
         {
+            bool IsMatrixArgument(Variable variable) => typesAndNames.Any(c => c.variable == variable && TensorCompilation.IsMatrixType(c.type));
             // A matrix has no compiled form, but an expression built out of matrices often
             // has a value that is an ordinary number -- [0, 1]T * [[a, b], [c, d]] * [1, 0]
             // is c -- and there is nothing to stop that being compiled. Nothing simplified
             // before compiling, so those failed along with the ones that genuinely cannot
             // be compiled (https://github.com/asc-community/AngouriMath/issues/425).
             // Only expressions that mention a matrix pay for this.
             if (expr.Nodes.Any(node => node is Entity.Matrix))
-                expr = expr.InnerSimplified;
+                expr = TensorCompilation.SimplifiedAroundMatrixArguments(expr, IsMatrixArgument);
             expr = CompiledMembership.Rewritten(expr);
             expr = Functions.ExponentialQuotients.Rewritten(expr);
             StatementPositions.Check(expr);
+            // Two matrices whose shapes are written into the expression and do not fit are refused here, rather than by
+            // a delegate that throws at every point it is called at. https://github.com/asc-community/AngouriMath/issues/526
+            TensorCompilation.CheckShapes(expr, IsMatrixArgument);
             if (MathS.Settings.RefuseNonRealCompilation && returnType is { } real && IsRealType(real)
                 && typesAndNames.All(c => IsRealType(c.type)) && !IsRealWhereverDefined(expr, typesAndNames))
                 throw new UncompilableNodeException(
@@ -42,6 +46,16 @@ internal static TDelegate Compile<TDelegate>(
             var functionArguments = subexpressionsCache.Select(c => c.Value).ToArray(); // copying
             var localVars = new List<ParameterExpression>();
             var variableAssignments = new List<Expression>();
+            // A matrix passed in as a two-dimensional array is read into the tensor the arithmetic runs on, once.
+            foreach (var (_, variable) in typesAndNames)
+                if (subexpressionsCache[variable] is var argument && TensorCompilation.IsMatrixArray(argument.Type))
+                {
+                    var asTensor = TensorCompilation.FromArrayArgument(argument);
+                    var local = Expression.Variable(asTensor.Type);
+                    localVars.Add(local);
+                    variableAssignments.Add(Expression.Assign(local, asTensor));
+                    subexpressionsCache[variable] = local;
+                }
 
             // Linq.Expression refuses a mismatch by throwing its own exceptions -- an
             // InvalidOperationException for `not x` where x is a double, an ArgumentException
@@ -54,10 +68,21 @@ internal static TDelegate Compile<TDelegate>(
             {
                 // Over Complex a constant is a Complex too, so that a function of it is taken over Complex: as a long,
                 // sqrt(-1) widened only to double and was NaN. https://github.com/asc-community/AngouriMath/issues/2235
-                var constantType = typesAndNames.Any(c => c.type == typeof(System.Numerics.Complex)) || returnType == typeof(System.Numerics.Complex)
+                // A matrix over Complex, as an argument or as the result, is a Complex the same way.
+                static bool IsComplex(Type? type)
+                    => type == typeof(System.Numerics.Complex) || (type is not null && TensorCompilation.ElementOfMatrixType(type) == typeof(System.Numerics.Complex));
+                // Over a matrix argument of one element type, a written matrix beside it is over that type too, so that
+                // the two combine: [[1, 0]] * A over double is not a matrix of long times a matrix of double.
+                var matrixElements = typesAndNames.Select(c => TensorCompilation.ElementOfMatrixType(c.type)).OfType<Type>().Distinct().ToList();
+                var constantType = typesAndNames.Any(c => IsComplex(c.type)) || IsComplex(returnType)
                     ? typeof(System.Numerics.Complex)
-                    : null;
+                    : matrixElements.Count == 1 ? matrixElements[0] : null;
                 var tree = BuildTree(expr, subexpressionsCache, variableAssignments, localVars, protocol, constantType);
+                // A matrix known to be 1 x 1 asked for as a number is its entry, as evaluation has it. Any other matrix
+                // is refused as a number by ConvertType.
+                if (TensorCompilation.IsTensor(tree.Type) && returnType is not null && !TensorCompilation.IsMatrixType(returnType)
+                    && TensorCompilation.IsOneByOne(expr, IsMatrixArgument))
+                    tree = TensorCompilation.ToScalar(tree);
                 var treeWithLocals = Expression.Block(localVars, variableAssignments.Append(tree));
                 Expression entireExpression = returnType is not null ? protocol.ConvertType(treeWithLocals, returnType) : treeWithLocals;
                 var finalLambda = Expression.Lambda<TDelegate>(entireExpression, functionArguments);