diff --git a/HISTORY.md b/HISTORY.md
index d656f33a97..a53bf5b556 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -1687,7 +1687,7 @@ Breaking changes:
 - Dropped support for (non-primitive) instances of `Number`, `Boolean`, and
   `String` from functions `clone` and `typeof`.
 - Dropped official support for IE9 (probably still works, but it's not tested).
-- Fixed #851: More consistent behavior of sqrt, nthRoot, and pow.
+- Fixed #851: More consistent behavior of `sqrt`, `nthRoot`, and `pow`.
   Thanks @dakotablair.
 - Fixed #1103: Calling `toTex` on node that contains `derivative` causing
   an exception. Thanks @joelhoover.
@@ -2520,7 +2520,7 @@ Non breaking changes:
 
 ## 2015-07-12, version 1.7.1
 
-- Fixed #397: Inaccuracies in nthRoot for very large values, and wrong results
+- Fixed #397: Inaccuracies in `nthRoot` for very large values, and wrong results
   for very small values. (backported from v2)
 - Fixed #405: Parser throws error when defining a function in a multiline
   expression.
@@ -3061,7 +3061,7 @@ Non breaking changes:
 
 - Extended the import function and some other minor improvements.
 - Fixed a bug in merging one dimensional vectors into a matrix.
-- Fixed a bug in function subtract, when subtracting a complex number from a
+- Fixed a bug in function `subtract`, when subtracting a complex number from a
   real number.
 
 ## 2013-05-10, version 0.8.1
diff --git a/docs/datatypes/matrices.md b/docs/datatypes/matrices.md
index 01178f7413..a1d40d37ee 100644
--- a/docs/datatypes/matrices.md
+++ b/docs/datatypes/matrices.md
@@ -12,10 +12,13 @@ dense and sparse matrices.
 Math.js supports two types of matrices:
 
 - `Array`, a regular JavaScript array. A multidimensional array can be created
-  by nesting arrays.
+  by nesting arrays. The following terminology will be used to define different kinds of nested arrays.
+  - **Recangular array**: All elements of an array are arrays of the same size. Like `[[1, 2], [3, 4]]`
+  - **Jagged arrays**: All elements of an array are arrays but not of the same size. Like `[[1, 2], [3, 4, 5]]`
+  - **Heterogeneous arrays**: If not all the elements inside an array are arrays. Like `[[1, 2], 3]`.
 - `Matrix`, a matrix implementation by math.js. A `Matrix` is an object wrapped
   around a regular JavaScript `Array`, providing utility functions for easy
-  matrix manipulation such as `subset`, `size`, `resize`, `clone`, and more.
+  matrix manipulation such as `subset`, `size`, `resize`, `clone`, and more. Nested arrays must be rectangular to be converted to a `Matrix`.
 
 In most cases, the type of matrix output from functions is determined by the
 function input: An `Array` as input will return an `Array`, a `Matrix` as input
@@ -224,6 +227,27 @@ If you have a matrix where the first dimension means `x` and the second
 means `y`, this will look confusing since `x` is printed as _column_ 
 (vertically) and `y` as _row_ (horizontally).
 
+## Non-rectangular arrays
+
+By nesting arrays it is possible to have arrays that are not rectangular, for example.
+```js
+[[1, 2], [3, 4]]    // rectangular of size [2, 2]
+[[1, 2], [3, 4, 5]] // jagged
+[[1, 2], 3]         // heterogeneous
+```
+
+Jagged and heterogeneous arrays can't be converted to a matrix, but many operations are available for them.
+```js
+const A = [[1, 2], 3]
+math.add(A, 1)
+math.map(A, a => a+1)
+math.forEach(A, a => console.log(a))
+```
+Some matrix functions expect a rectangular array and might provide unexpected results with non rectangular arrays, for example.
+```js
+math.size([[1, 2], [3]]) // [2, 2]
+```
+The process of validation for rectangularity is expensive and is mandatory to create a `Matrix`, thus there might be a performance benefit of not converting an `Array` to a `Matrix`.
 
 ## Resizing
 
diff --git a/package-lock.json b/package-lock.json
index 2637b244f6..39af161ed7 100644
--- a/package-lock.json
+++ b/package-lock.json
@@ -117,6 +117,7 @@
       "integrity": "sha512-CGOfOJqWjg2qW/Mb6zNsDm+u5vFQ8DxXfbM09z69p5Z6+mE1ikP2jUXw+j42Pf1XTYED2Rni5f95npYeuwMDQA==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "@babel/code-frame": "^7.29.0",
         "@babel/generator": "^7.29.0",
@@ -2302,6 +2303,7 @@
       "integrity": "sha512-FXx2pKgId/WyYo2jXw63kk7/+TY7u7AziEJxJAnSFzHlqTAS3Ync6SvgYAN/k4/PQpnnVuzoMuVnByKK2qp0ag==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "@types/estree": "*",
         "@types/json-schema": "*"
@@ -2366,6 +2368,7 @@
       "integrity": "sha512-ixiWrCSRi33uqBMRuICcKECW7rtgY43TbsHDpM2XK7lXispd48opW+0IXrBVxv9NMhaz/Ue9kyj6r3NTVyXm8A==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "undici-types": "~6.20.0"
       }
@@ -2415,6 +2418,7 @@
       "integrity": "sha512-4z2nCSBfVIMnbuu8uinj+f0o4qOeggYJLbjpPHka3KH1om7e+H9yLKTYgksTaHcGco+NClhhY2vyO3HsMH1RGw==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "@typescript-eslint/scope-manager": "8.55.0",
         "@typescript-eslint/types": "8.55.0",
@@ -2828,6 +2832,7 @@
       "integrity": "sha512-NZyJarBfL7nWwIq+FDL6Zp/yHEhePMNnnJ0y3qfieCrmNvYct8uvtiV41UvlSe6apAfk0fY1FbWx+NwfmpvtTg==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "bin": {
         "acorn": "bin/acorn"
       },
@@ -3812,6 +3817,7 @@
         }
       ],
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "baseline-browser-mapping": "^2.9.0",
         "caniuse-lite": "^1.0.30001759",
@@ -5262,6 +5268,7 @@
       "deprecated": "This version is no longer supported. Please see https://eslint.org/version-support for other options.",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "@eslint-community/eslint-utils": "^4.2.0",
         "@eslint-community/regexpp": "^4.6.1",
@@ -5347,6 +5354,7 @@
       "integrity": "sha512-NSWl5BFQWEPi1j4TjVNItzYV7dZXZ+wP6I6ZhrBGpChQhZRUaElihE9uRRkcbRnNb76UMKDF3r+WTmNcGPKsqw==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "bin": {
         "eslint-config-prettier": "bin/cli.js"
       },
@@ -5462,6 +5470,7 @@
       "integrity": "sha512-whOE1HFo/qJDyX4SnXzP4N6zOWn79WhnCUY/iDR0mPfQZO8wcYE4JClzI2oZrhBnnMUCBCHZhO6VQyoBU95mZA==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "@rtsao/scc": "^1.1.0",
         "array-includes": "^3.1.9",
@@ -5577,6 +5586,7 @@
       "integrity": "sha512-6TyDmZ1HXoFQXnhCTUjVFULReoBPOAjpuiKELMkeP40yffI/1ZRO+d9ug/VC6fqISo2WkuIBk3cvuRPALaWlOQ==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "@eslint-community/eslint-utils": "^4.4.0",
         "builtins": "^5.0.1",
@@ -5690,6 +5700,7 @@
       "integrity": "sha512-57Zzfw8G6+Gq7axm2Pdo3gW/Rx3h9Yywgn61uE/3elTCOePEHVrn2i5CdfBwA1BLK0Q0WqctICIUSqXZW/VprQ==",
       "dev": true,
       "license": "ISC",
+      "peer": true,
       "engines": {
         "node": "^12.22.0 || ^14.17.0 || >=16.0.0"
       },
@@ -8366,6 +8377,7 @@
       "integrity": "sha512-LrtUxbdvt1gOpo3gxG+VAJlJAEMhbWlM4YrFQgql98FwF7+K8K12LYO4hnDdUkNjeztYrOXEMqgTajSWgmtI/w==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "@colors/colors": "1.5.0",
         "body-parser": "^1.19.0",
@@ -10324,6 +10336,7 @@
       "integrity": "sha512-UOnG6LftzbdaHZcKoPFtOcCKztrQ57WkHDeRD9t/PTQtmT0NHSeWWepj6pS0z/N7+08BHFDQVUrfmfMRcZwbMg==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "bin": {
         "prettier": "bin/prettier.cjs"
       },
@@ -11018,6 +11031,7 @@
       "integrity": "sha512-B/gBuNg5SiMTrPkC+A2+cW0RszwxYmn6VYxB/inlBStS5nx6xHIt/ehKRhIMhqusl7a8LjQoZnjCs5vhwxOQ1g==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "fast-deep-equal": "^3.1.3",
         "fast-uri": "^3.0.1",
@@ -12079,6 +12093,7 @@
       "integrity": "sha512-5gTmgEY/sqK6gFXLIsQNH19lWb4ebPDLA4SdLP7dsWkIXHWlG66oPuVvXSGFPppYZz8ZDZq0dYYrbHfBCVUb1Q==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "engines": {
         "node": ">=12"
       },
@@ -12420,6 +12435,7 @@
       "integrity": "sha512-jl1vZzPDinLr9eUt3J/t7V6FgNEw9QjvBPdysz9KfQDD41fQrC2Y4vKQdiaUpFT4bXlb1RHhLpp8wtm6M5TgSw==",
       "dev": true,
       "license": "Apache-2.0",
+      "peer": true,
       "bin": {
         "tsc": "bin/tsc",
         "tsserver": "bin/tsserver"
@@ -12994,6 +13010,7 @@
       "integrity": "sha512-gX/dMkRQc7QOMzgTe6KsYFM7DxeIONQSui1s0n/0xht36HvrgbxtM1xBlgx596NbpHuQU8P7QpKwrZYwUX48nw==",
       "dev": true,
       "license": "MIT",
+      "peer": true,
       "dependencies": {
         "@types/eslint-scope": "^3.7.7",
         "@types/estree": "^1.0.8",
diff --git a/src/expression/parse.js b/src/expression/parse.js
index c300b8cf98..65de3eebba 100644
--- a/src/expression/parse.js
+++ b/src/expression/parse.js
@@ -82,20 +82,20 @@ export const createParse = /* #__PURE__ */ factory(name, dependencies, ({
    *
    * History:
    *
-   *     v0.9   Created
-   *     v0.13  Switched to one-based indices
-   *     v0.14  Added `[1,2;3,4]` notation for matrices
-   *     v0.18  Dropped the `function` keyword
-   *     v0.20  Added ternary conditional
-   *     v0.27  Allow multi-line expressions; allow functions that receive
-   *            unevaluated parameters (`rawArgs`)
+   *     v15.1  Add optional chaining operator
+   *     v14.8  Add nullish coalescing operator
+   *     v12.4  Allow trailing commas in matrices
+   *     v9.5   Support for calculations with percentages
+   *     v7.3   Supported binary, octal, and hexadecimal notation
    *     v3     Add object notation; allow assignments internal to other
    *            expressions
-   *     v7.3   Supported binary, octal, and hexadecimal notation
-   *     v9.5   Support for calculations with percentages
-   *     v12.4  Allow trailing commas in matrices
-   *     v14.8  Add nullish coalescing operator
-   *     v15.1  Add optional chaining operator
+   *     v0.27  Allow multi-line expressions; allow functions that receive
+   *            unevaluated parameters (`rawArgs`)
+   *     v0.20  Added ternary conditional
+   *     v0.18  Dropped the `function` keyword
+   *     v0.14  Added `[1,2;3,4]` notation for matrices
+   *     v0.13  Switched to one-based indices
+   *     v0.9   Created
    *
    * @param {string | string[] | Matrix} expr          Expression to be parsed
    * @param {{nodes: Object<string, Node>}} [options]  Available options:
diff --git a/src/function/arithmetic/add.js b/src/function/arithmetic/add.js
index 1acfdb30b9..c2c5fdbc53 100644
--- a/src/function/arithmetic/add.js
+++ b/src/function/arithmetic/add.js
@@ -10,19 +10,17 @@ const dependencies = [
   'matrix',
   'addScalar',
   'equalScalar',
-  'DenseMatrix',
-  'SparseMatrix',
-  'concat'
+  'DenseMatrix'
 ]
 
 export const createAdd = /* #__PURE__ */ factory(
   name,
   dependencies,
-  ({ typed, matrix, addScalar, equalScalar, DenseMatrix, SparseMatrix, concat }) => {
+  ({ typed, matrix, addScalar, equalScalar, DenseMatrix }) => {
     const matAlgo01xDSid = createMatAlgo01xDSid({ typed })
     const matAlgo04xSidSid = createMatAlgo04xSidSid({ typed, equalScalar })
     const matAlgo10xSids = createMatAlgo10xSids({ typed, DenseMatrix })
-    const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+    const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
     /**
      * Add two or more values, `x + y`.
      * For matrices, the function is evaluated element wise.
@@ -55,10 +53,9 @@ export const createAdd = /* #__PURE__ */ factory(
      *
      * History:
      *
-     *     v13     Handle bigint arguments
      *     v11.6   Support matrix broadcasting
+     *     v5.0.3  Support matrices with a `datatype` defined
      *     v3.8    Allow more than two arguments
-     *     v0.0.2  Created
      *
      * @param  {number | BigNumber | bigint | Fraction | Complex | Unit | Array | Matrix} x First value to add
      * @param  {number | BigNumber | bigint | Fraction | Complex | Unit | Array | Matrix} y Second value to add
diff --git a/src/function/arithmetic/dotDivide.js b/src/function/arithmetic/dotDivide.js
index 3587faf762..f071de70d7 100644
--- a/src/function/arithmetic/dotDivide.js
+++ b/src/function/arithmetic/dotDivide.js
@@ -13,17 +13,16 @@ const dependencies = [
   'equalScalar',
   'divideScalar',
   'DenseMatrix',
-  'concat',
   'SparseMatrix'
 ]
 
-export const createDotDivide = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, divideScalar, DenseMatrix, concat, SparseMatrix }) => {
+export const createDotDivide = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, divideScalar, DenseMatrix, SparseMatrix }) => {
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo07xSSf = createMatAlgo07xSSf({ typed, SparseMatrix })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Divide two matrices element wise. The function accepts both matrices and
@@ -47,6 +46,12 @@ export const createDotDivide = /* #__PURE__ */ factory(name, dependencies, ({ ty
    *
    *    divide, multiply, dotMultiply
    *
+   * History:
+   *
+   *    v14.0.0  Return a sparse matrix for sparse inputs
+   *    v11.6    Support matrix broadcasting
+   *    v0.23    Renamed from `edivide`
+   *
    * @param  {number | BigNumber | Fraction | Complex | Unit | Array | Matrix} x Numerator
    * @param  {number | BigNumber | Fraction | Complex | Unit | Array | Matrix} y Denominator
    * @return {number | BigNumber | Fraction | Complex | Unit | Array | Matrix}                    Quotient, `x ./ y`
diff --git a/src/function/arithmetic/dotMultiply.js b/src/function/arithmetic/dotMultiply.js
index d49287bc6b..6d6faa109d 100644
--- a/src/function/arithmetic/dotMultiply.js
+++ b/src/function/arithmetic/dotMultiply.js
@@ -9,15 +9,14 @@ const dependencies = [
   'typed',
   'matrix',
   'equalScalar',
-  'multiplyScalar',
-  'concat'
+  'multiplyScalar'
 ]
 
-export const createDotMultiply = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, multiplyScalar, concat }) => {
+export const createDotMultiply = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, multiplyScalar }) => {
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo09xS0Sf = createMatAlgo09xS0Sf({ typed, equalScalar })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Multiply two matrices element wise. The function accepts both matrices and
@@ -41,6 +40,11 @@ export const createDotMultiply = /* #__PURE__ */ factory(name, dependencies, ({
    *
    *    multiply, divide, dotDivide
    *
+   * History:
+   *
+   *    v11.6   Support matrix broadcasting
+   *    v0.23   Renamed from `emultiply`
+   *
    * @param  {number | BigNumber | Fraction | Complex | Unit | Array | Matrix} x Left hand value
    * @param  {number | BigNumber | Fraction | Complex | Unit | Array | Matrix} y Right hand value
    * @return {number | BigNumber | Fraction | Complex | Unit | Array | Matrix}                    Multiplication of `x` and `y`
diff --git a/src/function/arithmetic/dotPow.js b/src/function/arithmetic/dotPow.js
index a1c25f1522..91c7227ce8 100644
--- a/src/function/arithmetic/dotPow.js
+++ b/src/function/arithmetic/dotPow.js
@@ -12,16 +12,15 @@ const dependencies = [
   'matrix',
   'pow',
   'DenseMatrix',
-  'concat',
   'SparseMatrix'
 ]
 
-export const createDotPow = /* #__PURE__ */ factory(name, dependencies, ({ typed, equalScalar, matrix, pow, DenseMatrix, concat, SparseMatrix }) => {
+export const createDotPow = /* #__PURE__ */ factory(name, dependencies, ({ typed, equalScalar, matrix, pow, DenseMatrix, SparseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo07xSSf = createMatAlgo07xSSf({ typed, SparseMatrix })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   const powScalarSignatures = {}
   for (const signature in pow.signatures) {
@@ -52,6 +51,12 @@ export const createDotPow = /* #__PURE__ */ factory(name, dependencies, ({ typed
    *
    *    pow, sqrt, multiply
    *
+   * History:
+   *
+   *    v14.0.0  Return a sparse matrix for sparse inputs
+   *    v11.6    Support matrix broadcasting
+   *    v0.23    Renamed from `epow`
+   *
    * @param  {number | BigNumber | Complex | Unit | Array | Matrix} x  The base
    * @param  {number | BigNumber | Complex | Unit | Array | Matrix} y  The exponent
    * @return {number | BigNumber | Complex | Unit | Array | Matrix}                     The value of `x` to the power `y`
diff --git a/src/function/arithmetic/gcd.js b/src/function/arithmetic/gcd.js
index d370c498d6..296ab6dc67 100644
--- a/src/function/arithmetic/gcd.js
+++ b/src/function/arithmetic/gcd.js
@@ -32,7 +32,7 @@ export const createGcd = /* #__PURE__ */ factory(name, dependencies, ({ typed, m
   const matAlgo01xDSid = createMatAlgo01xDSid({ typed })
   const matAlgo04xSidSid = createMatAlgo04xSidSid({ typed, equalScalar })
   const matAlgo10xSids = createMatAlgo10xSids({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Calculate the greatest common divisor for two or more values or arrays.
@@ -56,6 +56,13 @@ export const createGcd = /* #__PURE__ */ factory(name, dependencies, ({ typed, m
    *
    *    lcm, xgcd
    *
+   * History:
+   *
+   *    v11.7  Accept arrays as input
+   *    v11.6  Support matrix broadcasting
+   *    v0.26  Implement BigNumber support
+   *    v0.6   Created
+   *
    * @param {... number | BigNumber | Fraction | Array | Matrix} args  Two or more integer numbers
    * @return {number | BigNumber | Fraction | Array | Matrix}                           The greatest common divisor
    */
diff --git a/src/function/arithmetic/lcm.js b/src/function/arithmetic/lcm.js
index d53dd3e799..9588dcdd98 100644
--- a/src/function/arithmetic/lcm.js
+++ b/src/function/arithmetic/lcm.js
@@ -9,15 +9,14 @@ const name = 'lcm'
 const dependencies = [
   'typed',
   'matrix',
-  'equalScalar',
-  'concat'
+  'equalScalar'
 ]
 
-export const createLcm = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, concat }) => {
+export const createLcm = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar }) => {
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo06xS0S0 = createMatAlgo06xS0S0({ typed, equalScalar })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   const lcmTypes = 'number | BigNumber | Fraction | Matrix | Array'
   const lcmManySignature = {}
@@ -56,6 +55,13 @@ export const createLcm = /* #__PURE__ */ factory(name, dependencies, ({ typed, m
    *
    *    gcd, xgcd
    *
+   * History:
+   *
+   *    v11.6  Support matrix broadcasting
+   *    v2.3   Support Fractions
+   *    v0.26  Implement BigNumber support
+   *    v0.6   Created
+   *
    * @param {... number | BigNumber | Array | Matrix} args  Two or more integer numbers
    * @return {number | BigNumber | Array | Matrix}                           The least common multiple
    */
diff --git a/src/function/arithmetic/mod.js b/src/function/arithmetic/mod.js
index 6c00537d38..c1801cfe7d 100644
--- a/src/function/arithmetic/mod.js
+++ b/src/function/arithmetic/mod.js
@@ -15,18 +15,17 @@ const dependencies = [
   'matrix',
   'equalScalar',
   'zeros',
-  'DenseMatrix',
-  'concat'
+  'DenseMatrix'
 ]
 
-export const createMod = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, round, matrix, equalScalar, zeros, DenseMatrix, concat }) => {
+export const createMod = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, round, matrix, equalScalar, zeros, DenseMatrix }) => {
   const floor = createFloor({ typed, config, round, matrix, equalScalar, zeros, DenseMatrix })
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo05xSfSf = createMatAlgo05xSfSf({ typed, equalScalar })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Calculates the modulus, the remainder of an integer division.
@@ -60,10 +59,11 @@ export const createMod = /* #__PURE__ */ factory(name, dependencies, ({ typed, c
    *
    * History:
    *
-   *    v13    Handle bigints
-   *    v11.6  Support matrix broadcasting
-   *    v2     Handle Fractions
-   *    v0.2   Created
+   *    v13     Handle bigints
+   *    v12.1   Support negative divisors for BigNumber and Fraction
+   *    v11.6   Support matrix broadcasting
+   *    v2      Handle Fractions
+   *    v0.2    Created
    *
    * @param  {number | BigNumber | bigint | Fraction | Array | Matrix} x Dividend
    * @param  {number | BigNumber | bigint | Fraction | Array | Matrix} y Divisor
diff --git a/src/function/arithmetic/nthRoot.js b/src/function/arithmetic/nthRoot.js
index aedbf140b4..69426a0889 100644
--- a/src/function/arithmetic/nthRoot.js
+++ b/src/function/arithmetic/nthRoot.js
@@ -11,16 +11,15 @@ const dependencies = [
   'typed',
   'matrix',
   'equalScalar',
-  'BigNumber',
-  'concat'
+  'BigNumber'
 ]
 
-export const createNthRoot = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, BigNumber, concat }) => {
+export const createNthRoot = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, BigNumber }) => {
   const matAlgo01xDSid = createMatAlgo01xDSid({ typed })
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo06xS0S0 = createMatAlgo06xS0S0({ typed, equalScalar })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Calculate the nth root of a value.
@@ -49,7 +48,10 @@ export const createNthRoot = /* #__PURE__ */ factory(name, dependencies, ({ type
    * History:
    *
    *     v11.6  Support matrix broadcasting
-   *     v2     Handle Complex
+   *     v5.0.0 Improve consistency with sqrt and pow
+   *     v3.1.1 Support negative roots of zero
+   *     v2.4.1 Support negative values like `nthRoot(-2, 3)`
+   *     v2     Support roots of Complex
    *     v1.1   Created
    *
    * @param {number | BigNumber | Array | Matrix | Complex} a
diff --git a/src/function/arithmetic/subtract.js b/src/function/arithmetic/subtract.js
index aa6732484b..0ee9303bc1 100644
--- a/src/function/arithmetic/subtract.js
+++ b/src/function/arithmetic/subtract.js
@@ -12,12 +12,10 @@ const dependencies = [
   'matrix',
   'equalScalar',
   'subtractScalar',
-  'unaryMinus',
-  'DenseMatrix',
-  'concat'
+  'DenseMatrix'
 ]
 
-export const createSubtract = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, subtractScalar, unaryMinus, DenseMatrix, concat }) => {
+export const createSubtract = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, subtractScalar, DenseMatrix }) => {
   // TODO: split function subtract in two: subtract and subtractScalar
 
   const matAlgo01xDSid = createMatAlgo01xDSid({ typed })
@@ -25,7 +23,7 @@ export const createSubtract = /* #__PURE__ */ factory(name, dependencies, ({ typ
   const matAlgo05xSfSf = createMatAlgo05xSfSf({ typed, equalScalar })
   const matAlgo10xSids = createMatAlgo10xSids({ typed, DenseMatrix })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Subtract two values, `x - y`.
@@ -53,6 +51,12 @@ export const createSubtract = /* #__PURE__ */ factory(name, dependencies, ({ typ
    *
    *    add
    *
+   * History:
+   *
+   *    v11.6  Support matrix broadcasting
+   *    v5.0.3 Support matrices with a `datatype` defined
+   *    v0.8.2 Support complex - real mixed subtraction
+   *
    * @param  {number | BigNumber | bigint | Fraction | Complex | Unit | Array | Matrix} x Initial value
    * @param  {number | BigNumber | bigint | Fraction | Complex | Unit | Array | Matrix} y Value to subtract from `x`
    * @return {number | BigNumber | bigint | Fraction | Complex | Unit | Array | Matrix} Subtraction of `x` and `y`
diff --git a/src/function/bitwise/bitAnd.js b/src/function/bitwise/bitAnd.js
index d6e169b09e..6a07dfeef5 100644
--- a/src/function/bitwise/bitAnd.js
+++ b/src/function/bitwise/bitAnd.js
@@ -10,15 +10,14 @@ const name = 'bitAnd'
 const dependencies = [
   'typed',
   'matrix',
-  'equalScalar',
-  'concat'
+  'equalScalar'
 ]
 
-export const createBitAnd = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, concat }) => {
+export const createBitAnd = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar }) => {
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo06xS0S0 = createMatAlgo06xS0S0({ typed, equalScalar })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Bitwise AND two values, `x & y`.
@@ -38,6 +37,11 @@ export const createBitAnd = /* #__PURE__ */ factory(name, dependencies, ({ typed
    *
    *    bitNot, bitOr, bitXor, leftShift, rightArithShift, rightLogShift
    *
+   * History:
+   *
+   *    v11.6  Support matrix broadcasting
+   *    v1.2   Created
+   *
    * @param  {number | BigNumber | bigint | Array | Matrix} x First value to and
    * @param  {number | BigNumber | bigint | Array | Matrix} y Second value to and
    * @return {number | BigNumber | bigint | Array | Matrix} AND of `x` and `y`
diff --git a/src/function/bitwise/bitOr.js b/src/function/bitwise/bitOr.js
index c0f7dea215..130abb207b 100644
--- a/src/function/bitwise/bitOr.js
+++ b/src/function/bitwise/bitOr.js
@@ -11,15 +11,14 @@ const dependencies = [
   'typed',
   'matrix',
   'equalScalar',
-  'DenseMatrix',
-  'concat'
+  'DenseMatrix'
 ]
 
-export const createBitOr = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, DenseMatrix, concat }) => {
+export const createBitOr = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, DenseMatrix }) => {
   const matAlgo01xDSid = createMatAlgo01xDSid({ typed })
   const matAlgo04xSidSid = createMatAlgo04xSidSid({ typed, equalScalar })
   const matAlgo10xSids = createMatAlgo10xSids({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Bitwise OR two values, `x | y`.
@@ -40,6 +39,11 @@ export const createBitOr = /* #__PURE__ */ factory(name, dependencies, ({ typed,
    *
    *    bitAnd, bitNot, bitXor, leftShift, rightArithShift, rightLogShift
    *
+   * History:
+   *
+   *    v11.6  Support matrix broadcasting
+   *    v1.2   Created
+   *
    * @param  {number | BigNumber | bigint | Array | Matrix} x First value to or
    * @param  {number | BigNumber | bigint | Array | Matrix} y Second value to or
    * @return {number | BigNumber | bigint | Array | Matrix} OR of `x` and `y`
diff --git a/src/function/bitwise/bitXor.js b/src/function/bitwise/bitXor.js
index a46acacc8a..29aeb59c1a 100644
--- a/src/function/bitwise/bitXor.js
+++ b/src/function/bitwise/bitXor.js
@@ -11,15 +11,14 @@ const dependencies = [
   'typed',
   'matrix',
   'DenseMatrix',
-  'concat',
   'SparseMatrix'
 ]
 
-export const createBitXor = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, DenseMatrix, concat, SparseMatrix }) => {
+export const createBitXor = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, DenseMatrix, SparseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo07xSSf = createMatAlgo07xSSf({ typed, SparseMatrix })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Bitwise XOR two values, `x ^ y`.
@@ -39,6 +38,12 @@ export const createBitXor = /* #__PURE__ */ factory(name, dependencies, ({ typed
    *
    *    bitAnd, bitNot, bitOr, leftShift, rightArithShift, rightLogShift
    *
+   * History:
+   *
+   *    v14.0.0 Return a sparse matrix for sparse inputs
+   *    v11.6   Support matrix broadcasting
+   *    v1.2    Created
+   *
    * @param  {number | BigNumber | bigint | Array | Matrix} x First value to xor
    * @param  {number | BigNumber | bigint | Array | Matrix} y Second value to xor
    * @return {number | BigNumber | bigint | Array | Matrix} XOR of `x` and `y`
diff --git a/src/function/bitwise/leftShift.js b/src/function/bitwise/leftShift.js
index 18b406fc61..6552d8ea8e 100644
--- a/src/function/bitwise/leftShift.js
+++ b/src/function/bitwise/leftShift.js
@@ -1,14 +1,15 @@
 import { createMatAlgo02xDS0 } from '../../type/matrix/utils/matAlgo02xDS0.js'
 import { createMatAlgo11xS0s } from '../../type/matrix/utils/matAlgo11xS0s.js'
 import { createMatAlgo14xDs } from '../../type/matrix/utils/matAlgo14xDs.js'
+import { createMatAlgo15xAs } from '../../type/matrix/utils/matAlgo15xAs.js'
 import { createMatAlgo01xDSid } from '../../type/matrix/utils/matAlgo01xDSid.js'
 import { createMatAlgo10xSids } from '../../type/matrix/utils/matAlgo10xSids.js'
 import { createMatAlgo08xS0Sid } from '../../type/matrix/utils/matAlgo08xS0Sid.js'
 import { factory } from '../../utils/factory.js'
 import { createMatrixAlgorithmSuite } from '../../type/matrix/utils/matrixAlgorithmSuite.js'
-import { createUseMatrixForArrayScalar } from './useMatrixForArrayScalar.js'
 import { leftShiftNumber } from '../../plain/number/index.js'
 import { leftShiftBigNumber } from '../../utils/bignumber/bitwise.js'
+import { deepMap, clone } from '../../utils/array.js'
 
 const name = 'leftShift'
 const dependencies = [
@@ -16,19 +17,18 @@ const dependencies = [
   'matrix',
   'equalScalar',
   'zeros',
-  'DenseMatrix',
-  'concat'
+  'DenseMatrix'
 ]
 
-export const createLeftShift = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, zeros, DenseMatrix, concat }) => {
+export const createLeftShift = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, zeros, DenseMatrix }) => {
   const matAlgo01xDSid = createMatAlgo01xDSid({ typed })
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo08xS0Sid = createMatAlgo08xS0Sid({ typed, equalScalar })
   const matAlgo10xSids = createMatAlgo10xSids({ typed, DenseMatrix })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
   const matAlgo14xDs = createMatAlgo14xDs({ typed })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
-  const useMatrixForArrayScalar = createUseMatrixForArrayScalar({ typed, matrix })
+  const matAlgo15xAs = createMatAlgo15xAs()
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Bitwise left logical shift of a value x by y number of bits, `x << y`.
@@ -49,6 +49,11 @@ export const createLeftShift = /* #__PURE__ */ factory(name, dependencies, ({ ty
    *
    *    leftShift, bitNot, bitOr, bitXor, rightArithShift, rightLogShift
    *
+   * History:
+   *
+   *    v11.6  Support matrix broadcasting
+   *    v1.2   Created
+   *
    * @param  {number | BigNumber | bigint | Array | Matrix} x Value to be shifted
    * @param  {number | BigNumber | bigint} y Amount of shifts
    * @return {number | BigNumber | bigint | Array | Matrix} `x` shifted left `y` times
@@ -78,6 +83,14 @@ export const createLeftShift = /* #__PURE__ */ factory(name, dependencies, ({ ty
         return matAlgo14xDs(x, y, self, false)
       }),
 
+      'Array, number | BigNumber': typed.referToSelf(self => (x, y) => {
+        // check scalar
+        if (equalScalar(y, 0)) {
+          return clone(x)
+        }
+        return matAlgo15xAs(x, y, self, false)
+      }),
+
       'number | BigNumber, SparseMatrix': typed.referToSelf(self => (x, y) => {
         // check scalar
         if (equalScalar(x, 0)) {
@@ -92,9 +105,15 @@ export const createLeftShift = /* #__PURE__ */ factory(name, dependencies, ({ ty
           return zeros(y.size(), y.storage())
         }
         return matAlgo14xDs(y, x, self, true)
+      }),
+      'number | BigNumber, Array': typed.referToSelf(self => (x, y) => {
+        // check scalar
+        if (equalScalar(x, 0)) {
+          return deepMap(y, () => x)
+        }
+        return matAlgo15xAs(y, x, self, true)
       })
     },
-    useMatrixForArrayScalar,
     matrixAlgorithmSuite({
       SS: matAlgo08xS0Sid,
       DS: matAlgo01xDSid,
diff --git a/src/function/bitwise/rightArithShift.js b/src/function/bitwise/rightArithShift.js
index c84a411d06..9242707fed 100644
--- a/src/function/bitwise/rightArithShift.js
+++ b/src/function/bitwise/rightArithShift.js
@@ -2,13 +2,14 @@ import { rightArithShiftBigNumber } from '../../utils/bignumber/bitwise.js'
 import { createMatAlgo02xDS0 } from '../../type/matrix/utils/matAlgo02xDS0.js'
 import { createMatAlgo11xS0s } from '../../type/matrix/utils/matAlgo11xS0s.js'
 import { createMatAlgo14xDs } from '../../type/matrix/utils/matAlgo14xDs.js'
+import { createMatAlgo15xAs } from '../../type/matrix/utils/matAlgo15xAs.js'
 import { createMatAlgo01xDSid } from '../../type/matrix/utils/matAlgo01xDSid.js'
 import { createMatAlgo10xSids } from '../../type/matrix/utils/matAlgo10xSids.js'
 import { createMatAlgo08xS0Sid } from '../../type/matrix/utils/matAlgo08xS0Sid.js'
 import { factory } from '../../utils/factory.js'
 import { createMatrixAlgorithmSuite } from '../../type/matrix/utils/matrixAlgorithmSuite.js'
-import { createUseMatrixForArrayScalar } from './useMatrixForArrayScalar.js'
 import { rightArithShiftNumber } from '../../plain/number/index.js'
+import { clone, deepMap } from '../../utils/array.js'
 
 const name = 'rightArithShift'
 const dependencies = [
@@ -16,19 +17,18 @@ const dependencies = [
   'matrix',
   'equalScalar',
   'zeros',
-  'DenseMatrix',
-  'concat'
+  'DenseMatrix'
 ]
 
-export const createRightArithShift = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, zeros, DenseMatrix, concat }) => {
+export const createRightArithShift = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, zeros, DenseMatrix }) => {
   const matAlgo01xDSid = createMatAlgo01xDSid({ typed })
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo08xS0Sid = createMatAlgo08xS0Sid({ typed, equalScalar })
   const matAlgo10xSids = createMatAlgo10xSids({ typed, DenseMatrix })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
   const matAlgo14xDs = createMatAlgo14xDs({ typed })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
-  const useMatrixForArrayScalar = createUseMatrixForArrayScalar({ typed, matrix })
+  const matAlgo15xAs = createMatAlgo15xAs()
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Bitwise right arithmetic shift of a value x by y number of bits, `x >> y`.
@@ -49,6 +49,11 @@ export const createRightArithShift = /* #__PURE__ */ factory(name, dependencies,
    *
    *    bitAnd, bitNot, bitOr, bitXor, rightArithShift, rightLogShift
    *
+   * History:
+   *
+   *    v11.6  Support matrix broadcasting
+   *    v1.2   Created
+   *
    * @param  {number | BigNumber | bigint | Array | Matrix} x Value to be shifted
    * @param  {number | BigNumber | bigint} y Amount of shifts
    * @return {number | BigNumber | bigint | Array | Matrix} `x` zero-filled shifted right `y` times
@@ -78,6 +83,14 @@ export const createRightArithShift = /* #__PURE__ */ factory(name, dependencies,
         return matAlgo14xDs(x, y, self, false)
       }),
 
+      'Array, number | BigNumber': typed.referToSelf(self => (x, y) => {
+        // check scalar
+        if (equalScalar(y, 0)) {
+          return clone(x)
+        }
+        return matAlgo15xAs(x, y, self, false)
+      }),
+
       'number | BigNumber, SparseMatrix': typed.referToSelf(self => (x, y) => {
         // check scalar
         if (equalScalar(x, 0)) {
@@ -92,9 +105,15 @@ export const createRightArithShift = /* #__PURE__ */ factory(name, dependencies,
           return zeros(y.size(), y.storage())
         }
         return matAlgo14xDs(y, x, self, true)
+      }),
+      'number | BigNumber, Array': typed.referToSelf(self => (x, y) => {
+        // check scalar
+        if (equalScalar(x, 0)) {
+          return deepMap(y, () => x)
+        }
+        return matAlgo15xAs(y, x, self, true)
       })
     },
-    useMatrixForArrayScalar,
     matrixAlgorithmSuite({
       SS: matAlgo08xS0Sid,
       DS: matAlgo01xDSid,
diff --git a/src/function/bitwise/rightLogShift.js b/src/function/bitwise/rightLogShift.js
index 55cef7d355..c90ecfbbfe 100644
--- a/src/function/bitwise/rightLogShift.js
+++ b/src/function/bitwise/rightLogShift.js
@@ -4,10 +4,11 @@ import { createMatAlgo14xDs } from '../../type/matrix/utils/matAlgo14xDs.js'
 import { createMatAlgo01xDSid } from '../../type/matrix/utils/matAlgo01xDSid.js'
 import { createMatAlgo10xSids } from '../../type/matrix/utils/matAlgo10xSids.js'
 import { createMatAlgo08xS0Sid } from '../../type/matrix/utils/matAlgo08xS0Sid.js'
+import { createMatAlgo15xAs } from '../../type/matrix/utils/matAlgo15xAs.js'
 import { factory } from '../../utils/factory.js'
 import { createMatrixAlgorithmSuite } from '../../type/matrix/utils/matrixAlgorithmSuite.js'
 import { rightLogShiftNumber } from '../../plain/number/index.js'
-import { createUseMatrixForArrayScalar } from './useMatrixForArrayScalar.js'
+import { deepMap, clone } from '../../utils/array.js'
 
 const name = 'rightLogShift'
 const dependencies = [
@@ -15,19 +16,18 @@ const dependencies = [
   'matrix',
   'equalScalar',
   'zeros',
-  'DenseMatrix',
-  'concat'
+  'DenseMatrix'
 ]
 
-export const createRightLogShift = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, zeros, DenseMatrix, concat }) => {
+export const createRightLogShift = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, zeros, DenseMatrix }) => {
   const matAlgo01xDSid = createMatAlgo01xDSid({ typed })
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo08xS0Sid = createMatAlgo08xS0Sid({ typed, equalScalar })
   const matAlgo10xSids = createMatAlgo10xSids({ typed, DenseMatrix })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
   const matAlgo14xDs = createMatAlgo14xDs({ typed })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
-  const useMatrixForArrayScalar = createUseMatrixForArrayScalar({ typed, matrix })
+  const matAlgo15xAs = createMatAlgo15xAs()
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Bitwise right logical shift of value x by y number of bits, `x >>> y`.
@@ -48,6 +48,11 @@ export const createRightLogShift = /* #__PURE__ */ factory(name, dependencies, (
    *
    *    bitAnd, bitNot, bitOr, bitXor, leftShift, rightLogShift
    *
+   * History:
+   *
+   *    v11.6  Support matrix broadcasting
+   *    v1.2   Created
+   *
    * @param  {number | Array | Matrix} x Value to be shifted
    * @param  {number} y Amount of shifts
    * @return {number | Array | Matrix} `x` zero-filled shifted right `y` times
@@ -76,6 +81,14 @@ export const createRightLogShift = /* #__PURE__ */ factory(name, dependencies, (
         return matAlgo14xDs(x, y, self, false)
       }),
 
+      'Array, number | BigNumber': typed.referToSelf(self => (x, y) => {
+        // check scalar
+        if (equalScalar(y, 0)) {
+          return clone(x)
+        }
+        return matAlgo15xAs(x, y, self, false)
+      }),
+
       'number | BigNumber, SparseMatrix': typed.referToSelf(self => (x, y) => {
         // check scalar
         if (equalScalar(x, 0)) {
@@ -90,9 +103,16 @@ export const createRightLogShift = /* #__PURE__ */ factory(name, dependencies, (
           return zeros(y.size(), y.storage())
         }
         return matAlgo14xDs(y, x, self, true)
+      }),
+
+      'number | BigNumber, Array': typed.referToSelf(self => (x, y) => {
+        // check scalar
+        if (equalScalar(x, 0)) {
+          return deepMap(y, () => x)
+        }
+        return matAlgo15xAs(y, x, self, true)
       })
     },
-    useMatrixForArrayScalar,
     matrixAlgorithmSuite({
       SS: matAlgo08xS0Sid,
       DS: matAlgo01xDSid,
diff --git a/src/function/bitwise/useMatrixForArrayScalar.js b/src/function/bitwise/useMatrixForArrayScalar.js
deleted file mode 100644
index ff7f231e2f..0000000000
--- a/src/function/bitwise/useMatrixForArrayScalar.js
+++ /dev/null
@@ -1,15 +0,0 @@
-import { factory } from '../../utils/factory.js'
-
-export const createUseMatrixForArrayScalar = /* #__PURE__ */ factory('useMatrixForArrayScalar', ['typed', 'matrix'], ({ typed, matrix }) => ({
-  'Array, number': typed.referTo('DenseMatrix, number',
-    selfDn => (x, y) => selfDn(matrix(x), y).valueOf()),
-
-  'Array, BigNumber': typed.referTo('DenseMatrix, BigNumber',
-    selfDB => (x, y) => selfDB(matrix(x), y).valueOf()),
-
-  'number, Array': typed.referTo('number, DenseMatrix',
-    selfnD => (x, y) => selfnD(x, matrix(y)).valueOf()),
-
-  'BigNumber, Array': typed.referTo('BigNumber, DenseMatrix',
-    selfBD => (x, y) => selfBD(x, matrix(y)).valueOf())
-}))
diff --git a/src/function/logical/and.js b/src/function/logical/and.js
index 827751d042..2d98109502 100644
--- a/src/function/logical/and.js
+++ b/src/function/logical/and.js
@@ -12,16 +12,15 @@ const dependencies = [
   'matrix',
   'equalScalar',
   'zeros',
-  'not',
-  'concat'
+  'not'
 ]
 
-export const createAnd = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, zeros, not, concat }) => {
+export const createAnd = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, zeros, not }) => {
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo06xS0S0 = createMatAlgo06xS0S0({ typed, equalScalar })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
   const matAlgo14xDs = createMatAlgo14xDs({ typed })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Logical `and`. Test whether two values are both defined with a nonzero/nonempty value.
@@ -46,6 +45,11 @@ export const createAnd = /* #__PURE__ */ factory(name, dependencies, ({ typed, m
    *
    *    not, or, xor
    *
+   * History:
+   *    v12.2  Lazy evaluation
+   *    v11.6  Support matrix broadcasting
+   *    v1.2   Created
+   *
    * @param  {number | BigNumber | bigint | Complex | Unit | Array | Matrix} x First value to check
    * @param  {number | BigNumber | bigint | Complex | Unit | Array | Matrix} y Second value to check
    * @return {boolean | Array | Matrix}
diff --git a/src/function/logical/or.js b/src/function/logical/or.js
index 12d6c5b48d..d8f369dfcb 100644
--- a/src/function/logical/or.js
+++ b/src/function/logical/or.js
@@ -10,15 +10,14 @@ const dependencies = [
   'typed',
   'matrix',
   'equalScalar',
-  'DenseMatrix',
-  'concat'
+  'DenseMatrix'
 ]
 
-export const createOr = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, DenseMatrix, concat }) => {
+export const createOr = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, DenseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo05xSfSf = createMatAlgo05xSfSf({ typed, equalScalar })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Logical `or`. Test if at least one value is defined with a nonzero/nonempty value.
@@ -43,6 +42,12 @@ export const createOr = /* #__PURE__ */ factory(name, dependencies, ({ typed, ma
    *
    *    and, not, xor
    *
+   * History:
+   *
+   *    v12.2  Lazy evaluation
+   *    v11.6  Support matrix broadcasting
+   *    v1.2   Created
+   *
    * @param  {number | BigNumber | bigint | Complex | Unit | Array | Matrix} x First value to check
    * @param  {number | BigNumber | bigint | Complex | Unit | Array | Matrix} y Second value to check
    * @return {boolean | Array | Matrix}
diff --git a/src/function/logical/xor.js b/src/function/logical/xor.js
index 1ae3a05f5b..659cbe2a1d 100644
--- a/src/function/logical/xor.js
+++ b/src/function/logical/xor.js
@@ -10,15 +10,14 @@ const dependencies = [
   'typed',
   'matrix',
   'DenseMatrix',
-  'concat',
   'SparseMatrix'
 ]
 
-export const createXor = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, DenseMatrix, concat, SparseMatrix }) => {
+export const createXor = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, DenseMatrix, SparseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo07xSSf = createMatAlgo07xSSf({ typed, SparseMatrix })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Logical `xor`. Test whether one and only one value is defined with a nonzero/nonempty value.
@@ -43,6 +42,12 @@ export const createXor = /* #__PURE__ */ factory(name, dependencies, ({ typed, m
    *
    *    and, not, or
    *
+   * History:
+   *
+   *    v14.0.0 Return a sparse matrix for sparse inputs
+   *    v11.6   Support matrix broadcasting
+   *    v1.2    Created
+   *
    * @param  {number | BigNumber | bigint | Complex | Unit | Array | Matrix} x First value to check
    * @param  {number | BigNumber | bigint | Complex | Unit | Array | Matrix} y Second value to check
    * @return {boolean | Array | Matrix}
diff --git a/src/function/relational/compare.js b/src/function/relational/compare.js
index 41a4c4adf4..a96db151d2 100644
--- a/src/function/relational/compare.js
+++ b/src/function/relational/compare.js
@@ -15,15 +15,14 @@ const dependencies = [
   'equalScalar',
   'BigNumber',
   'Fraction',
-  'DenseMatrix',
-  'concat'
+  'DenseMatrix'
 ]
 
-export const createCompare = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, equalScalar, matrix, BigNumber, Fraction, DenseMatrix, concat }) => {
+export const createCompare = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, equalScalar, matrix, BigNumber, Fraction, DenseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo05xSfSf = createMatAlgo05xSfSf({ typed, equalScalar })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
   const compareUnits = createCompareUnits({ typed })
 
   /**
@@ -60,9 +59,10 @@ export const createCompare = /* #__PURE__ */ factory(name, dependencies, ({ type
    *
    * History:
    *
-   *    v0.19  Created
-   *    v4     Changed to compare strings by numerical value
    *    v13    Change to use separate relative and absolute tolerances
+   *    v11.6  Support matrix broadcasting
+   *    v4     Changed to compare strings by numerical value
+   *    v0.19  Created
    *
    * @param  {number | BigNumber | bigint | Fraction | Unit | string | Array | Matrix} x First value to compare
    * @param  {number | BigNumber | bigint | Fraction | Unit | string | Array | Matrix} y Second value to compare
diff --git a/src/function/relational/compareText.js b/src/function/relational/compareText.js
index 310b55cd0e..754db4ec55 100644
--- a/src/function/relational/compareText.js
+++ b/src/function/relational/compareText.js
@@ -5,14 +5,13 @@ import { createMatrixAlgorithmSuite } from '../../type/matrix/utils/matrixAlgori
 const name = 'compareText'
 const dependencies = [
   'typed',
-  'matrix',
-  'concat'
+  'matrix'
 ]
 
 _compareText.signature = 'any, any'
 
-export const createCompareText = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, concat }) => {
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+export const createCompareText = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix }) => {
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Compare two strings lexically. Comparison is case sensitive.
diff --git a/src/function/relational/equal.js b/src/function/relational/equal.js
index c606a06e98..1e88cfb597 100644
--- a/src/function/relational/equal.js
+++ b/src/function/relational/equal.js
@@ -13,7 +13,7 @@ const dependencies = [
   'SparseMatrix'
 ]
 
-export const createEqual = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, DenseMatrix, concat, SparseMatrix }) => {
+export const createEqual = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, DenseMatrix, SparseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo07xSSf = createMatAlgo07xSSf({ typed, SparseMatrix })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
@@ -64,9 +64,10 @@ export const createEqual = /* #__PURE__ */ factory(name, dependencies, ({ typed,
    *
    * History:
    *
-   *    v13    Handle bigints
+   *    v14    Return a sparse matrix for sparse inputs
+   *    v13    Handle bigints; use separate relative and absolute tolerances
    *    v11.6  Support matrix broadcasting
-   *    v4     Compare strings by their numerical value
+   *    v4     Compare strings by numerical, allow tolerance in comparison
    *    v0.24  Handle `null` and `undefined`
    *    v0.23  Compare collections elementwise
    *    v0.20  Compare floating-point numbers within epsilon, allowing roundoff
diff --git a/src/function/relational/larger.js b/src/function/relational/larger.js
index 8507d21511..bf0c484f56 100644
--- a/src/function/relational/larger.js
+++ b/src/function/relational/larger.js
@@ -14,15 +14,14 @@ const dependencies = [
   'bignumber',
   'matrix',
   'DenseMatrix',
-  'concat',
   'SparseMatrix'
 ]
 
-export const createLarger = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, bignumber, matrix, DenseMatrix, concat, SparseMatrix }) => {
+export const createLarger = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, bignumber, matrix, DenseMatrix, SparseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo07xSSf = createMatAlgo07xSSf({ typed, SparseMatrix })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
   const compareUnits = createCompareUnits({ typed })
 
   /**
@@ -52,6 +51,15 @@ export const createLarger = /* #__PURE__ */ factory(name, dependencies, ({ typed
    *
    *    equal, unequal, smaller, smallerEq, largerEq, compare
    *
+   * History:
+   *
+   *    v14.1  Support bigints
+   *    v14    Return a sparse matrix for sparse inputs
+   *    v11.6  Support matrix broadcasting
+   *    v4     Compare strings by numerical value, allow tolerance in comparison
+   *    v3     Use nearly-equal comparison for BigNumbers
+   *    v0.15  Removed support for complex numbers
+   *
    * @param  {number | BigNumber | bigint | Fraction | boolean | Unit | string | Array | Matrix} x First value to compare
    * @param  {number | BigNumber | bigint | Fraction | boolean | Unit | string | Array | Matrix} y Second value to compare
    * @return {boolean | Array | Matrix} Returns true when the x is larger than y, else returns false
diff --git a/src/function/relational/largerEq.js b/src/function/relational/largerEq.js
index 6760df7f6c..0c24d86cad 100644
--- a/src/function/relational/largerEq.js
+++ b/src/function/relational/largerEq.js
@@ -13,15 +13,14 @@ const dependencies = [
   'config',
   'matrix',
   'DenseMatrix',
-  'concat',
   'SparseMatrix'
 ]
 
-export const createLargerEq = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, matrix, DenseMatrix, concat, SparseMatrix }) => {
+export const createLargerEq = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, matrix, DenseMatrix, SparseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo07xSSf = createMatAlgo07xSSf({ typed, SparseMatrix })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
   const compareUnits = createCompareUnits({ typed })
 
   /**
@@ -47,6 +46,16 @@ export const createLargerEq = /* #__PURE__ */ factory(name, dependencies, ({ typ
    *
    *    equal, unequal, smaller, smallerEq, larger, compare
    *
+   * History:
+   *
+   *    v14.0.0 Return a sparse matrix for sparse inputs
+   *    v11.6  Support matrix broadcasting
+   *    v4     Compare strings by numerical value, allow tolerance in comparison
+   *    v3     Use nearly-equal comparison for BigNumbers
+   *    0.23   Renamed from `largereq`to `largerEq`
+   *    v0.15  Removed support for complex numbers
+   *    v0.2   Created
+   *
    * @param  {number | BigNumber | bigint | Fraction | boolean | Unit | string | Array | Matrix} x First value to compare
    * @param  {number | BigNumber | bigint | Fraction | boolean | Unit | string | Array | Matrix} y Second value to compare
    * @return {boolean | Array | Matrix} Returns true when the x is larger or equal to y, else returns false
diff --git a/src/function/relational/smaller.js b/src/function/relational/smaller.js
index 87807f0240..f53cb16bc8 100644
--- a/src/function/relational/smaller.js
+++ b/src/function/relational/smaller.js
@@ -14,15 +14,14 @@ const dependencies = [
   'bignumber',
   'matrix',
   'DenseMatrix',
-  'concat',
   'SparseMatrix'
 ]
 
-export const createSmaller = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, bignumber, matrix, DenseMatrix, concat, SparseMatrix }) => {
+export const createSmaller = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, bignumber, matrix, DenseMatrix, SparseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo07xSSf = createMatAlgo07xSSf({ typed, SparseMatrix })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
   const compareUnits = createCompareUnits({ typed })
 
   /**
@@ -52,6 +51,14 @@ export const createSmaller = /* #__PURE__ */ factory(name, dependencies, ({ type
    *
    *    equal, unequal, smallerEq, smaller, smallerEq, compare
    *
+   * History:
+   *
+   *    v14.1  Support bigints
+   *    v14    Return a sparse matrix for sparse inputs
+   *    v11.6  Support matrix broadcasting
+   *    v4     Compare strings by numerical value, allow tolerance in comparison
+   *    v3     Use nearly-equal comparison for BigNumbers
+   *
    * @param  {number | BigNumber | bigint | Fraction | boolean | Unit | string | Array | Matrix} x First value to compare
    * @param  {number | BigNumber | bigint | Fraction | boolean | Unit | string | Array | Matrix} y Second value to compare
    * @return {boolean | Array | Matrix} Returns true when the x is smaller than y, else returns false
diff --git a/src/function/relational/smallerEq.js b/src/function/relational/smallerEq.js
index b314418a4c..3a035c9225 100644
--- a/src/function/relational/smallerEq.js
+++ b/src/function/relational/smallerEq.js
@@ -13,15 +13,14 @@ const dependencies = [
   'config',
   'matrix',
   'DenseMatrix',
-  'concat',
   'SparseMatrix'
 ]
 
-export const createSmallerEq = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, matrix, DenseMatrix, concat, SparseMatrix }) => {
+export const createSmallerEq = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, matrix, DenseMatrix, SparseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo07xSSf = createMatAlgo07xSSf({ typed, SparseMatrix })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
   const compareUnits = createCompareUnits({ typed })
 
   /**
@@ -47,6 +46,14 @@ export const createSmallerEq = /* #__PURE__ */ factory(name, dependencies, ({ ty
    *
    *    equal, unequal, smaller, larger, largerEq, compare
    *
+   * History:
+   *
+   *    v14.0.0 Return a sparse matrix for sparse inputs
+   *    v11.6  Support matrix broadcasting
+   *    v4     Compare strings by numerical value, allow tolerance in comparison
+   *    v0.23  Renamed from `smalleq`to `smallerEq`
+   *    v0.15  Removed support for complex numbers
+   *
    * @param  {number | BigNumber | bigint | Fraction | boolean | Unit | string | Array | Matrix} x First value to compare
    * @param  {number | BigNumber | bigint | Fraction | boolean | Unit | string | Array | Matrix} y Second value to compare
    * @return {boolean | Array | Matrix} Returns true when the x is smaller than y, else returns false
diff --git a/src/function/relational/unequal.js b/src/function/relational/unequal.js
index f0a15a2041..696b878c38 100644
--- a/src/function/relational/unequal.js
+++ b/src/function/relational/unequal.js
@@ -7,19 +7,17 @@ import { createMatrixAlgorithmSuite } from '../../type/matrix/utils/matrixAlgori
 const name = 'unequal'
 const dependencies = [
   'typed',
-  'config',
   'equalScalar',
   'matrix',
   'DenseMatrix',
-  'concat',
   'SparseMatrix'
 ]
 
-export const createUnequal = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, equalScalar, matrix, DenseMatrix, concat, SparseMatrix }) => {
+export const createUnequal = /* #__PURE__ */ factory(name, dependencies, ({ typed, equalScalar, matrix, DenseMatrix, SparseMatrix }) => {
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo07xSSf = createMatAlgo07xSSf({ typed, SparseMatrix })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Test whether two values are unequal.
@@ -66,11 +64,12 @@ export const createUnequal = /* #__PURE__ */ factory(name, dependencies, ({ type
    *
    * History:
    *
-   *    v13    Handle bigints
-   *    v11.6  Support matrix broadcasting
-   *    v4     Compare strings by their numeric values
-   *    v0.24  Handle `null` and `undefined`
-   *    v0.2   Created
+   *    v14.0.0 Return a sparse matrix for sparse inputs
+   *    v13     Handle bigints
+   *    v11.6   Support matrix broadcasting
+   *    v4      Compare strings by their numeric values
+   *    v0.24   Handle `null` and `undefined`
+   *    v0.2    Created
    *
    * @param  {number | BigNumber | Fraction | boolean | Complex | Unit | string | Array | Matrix | undefined} x First value to compare
    * @param  {number | BigNumber | Fraction | boolean | Complex | Unit | string | Array | Matrix | undefined} y Second value to compare
diff --git a/src/function/trigonometry/atan2.js b/src/function/trigonometry/atan2.js
index ec0eb544d3..899cc8de73 100644
--- a/src/function/trigonometry/atan2.js
+++ b/src/function/trigonometry/atan2.js
@@ -12,17 +12,16 @@ const dependencies = [
   'matrix',
   'equalScalar',
   'BigNumber',
-  'DenseMatrix',
-  'concat'
+  'DenseMatrix'
 ]
 
-export const createAtan2 = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, BigNumber, DenseMatrix, concat }) => {
+export const createAtan2 = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, equalScalar, BigNumber, DenseMatrix }) => {
   const matAlgo02xDS0 = createMatAlgo02xDS0({ typed, equalScalar })
   const matAlgo03xDSf = createMatAlgo03xDSf({ typed })
   const matAlgo09xS0Sf = createMatAlgo09xS0Sf({ typed, equalScalar })
   const matAlgo11xS0s = createMatAlgo11xS0s({ typed, equalScalar })
   const matAlgo12xSfs = createMatAlgo12xSfs({ typed, DenseMatrix })
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Calculate the inverse tangent function with two arguments, y/x.
@@ -50,6 +49,11 @@ export const createAtan2 = /* #__PURE__ */ factory(name, dependencies, ({ typed,
    *
    *    tan, atan, sin, cos
    *
+  * History:
+  *
+  *    v11.6  Support matrix broadcasting
+  *    v1.5   Support BigNumber arguments
+  *
    * @param {number | Array | Matrix} y  Second dimension
    * @param {number | Array | Matrix} x  First dimension
    * @return {number | Array | Matrix} Four-quadrant inverse tangent
diff --git a/src/function/unit/to.js b/src/function/unit/to.js
index e0705a5bd1..f362691b4c 100644
--- a/src/function/unit/to.js
+++ b/src/function/unit/to.js
@@ -4,12 +4,11 @@ import { createMatrixAlgorithmSuite } from '../../type/matrix/utils/matrixAlgori
 const name = 'to'
 const dependencies = [
   'typed',
-  'matrix',
-  'concat'
+  'matrix'
 ]
 
-export const createTo = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, concat }) => {
-  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix, concat })
+export const createTo = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix }) => {
+  const matrixAlgorithmSuite = createMatrixAlgorithmSuite({ typed, matrix })
 
   /**
    * Change the unit of a value.
@@ -30,6 +29,12 @@ export const createTo = /* #__PURE__ */ factory(name, dependencies, ({ typed, ma
    *
    *    unit
    *
+   * History:
+   *
+   *    v11.6  Support matrix broadcasting
+   *    v2.2   Fix simplified unit conversions
+   *    v0.18  Renamed from `in` to `to`
+   *
    * @param {Unit | Array | Matrix} x     The unit to be converted.
    * @param {Unit | Array | Matrix} unit  New unit. Can be a string like "cm"
    *                                      or a unit without value.
diff --git a/src/type/matrix/utils/matAlgo14xDs.js b/src/type/matrix/utils/matAlgo14xDs.js
index 342b2b19d9..14acc1ba0e 100644
--- a/src/type/matrix/utils/matAlgo14xDs.js
+++ b/src/type/matrix/utils/matAlgo14xDs.js
@@ -1,5 +1,4 @@
 import { factory } from '../../../utils/factory.js'
-import { clone } from '../../../utils/object.js'
 
 const name = 'matAlgo14xDs'
 const dependencies = ['typed']
@@ -22,54 +21,19 @@ export const createMatAlgo14xDs = /* #__PURE__ */ factory(name, dependencies, ({
    */
   return function matAlgo14xDs (a, b, callback, inverse) {
     // a arrays
-    const adata = a._data
-    const asize = a._size
     const adt = a._datatype
 
-    // datatype
-    let dt
     // callback signature to use
     let cf = callback
 
     // process data types
     if (typeof adt === 'string') {
-      // datatype
-      dt = adt
       // convert b to the same datatype
-      b = typed.convert(b, dt)
+      b = typed.convert(b, adt)
       // callback
-      cf = typed.find(callback, [dt, dt])
+      cf = typed.find(callback, [adt, adt])
     }
 
-    // populate cdata, iterate through dimensions
-    const cdata = asize.length > 0 ? _iterate(cf, 0, asize, asize[0], adata, b, inverse) : []
-
-    // c matrix
-    return a.createDenseMatrix({
-      data: cdata,
-      size: clone(asize),
-      datatype: dt
-    })
-  }
-
-  // recursive function
-  function _iterate (f, level, s, n, av, bv, inverse) {
-    // initialize array for this level
-    const cv = []
-    // check we reach the last level
-    if (level === s.length - 1) {
-      // loop arrays in last level
-      for (let i = 0; i < n; i++) {
-        // invoke callback and store value
-        cv[i] = inverse ? f(bv, av[i]) : f(av[i], bv)
-      }
-    } else {
-      // iterate current level
-      for (let j = 0; j < n; j++) {
-        // iterate next level
-        cv[j] = _iterate(f, level + 1, s, s[level + 1], av[j], bv, inverse)
-      }
-    }
-    return cv
+    return inverse ? a.map(v => cf(b, v), false, true) : a.map(v => cf(v, b), false, true)
   }
 })
diff --git a/src/type/matrix/utils/matAlgo15xAs.js b/src/type/matrix/utils/matAlgo15xAs.js
new file mode 100644
index 0000000000..163c1fad76
--- /dev/null
+++ b/src/type/matrix/utils/matAlgo15xAs.js
@@ -0,0 +1,27 @@
+import { factory } from '../../../utils/factory.js'
+import { deepMap as map } from '../../../utils/array.js'
+
+const name = 'matAlgo15xAs'
+const dependencies = []
+
+export const createMatAlgo15xAs = /* #__PURE__ */ factory(name, dependencies, () => {
+  /**
+   * Deeply iterates over Array items and invokes the callback function
+   * f(Aij..z, b) for each nested scalar.
+   * The callback function is invoked once for each item.
+   *
+   * C(i,j,...z) = f(Aij..z, b)
+   *
+   * @param {Array}    a                 The Array instance (A)
+   * @param {Scalar}   b                 The Scalar value
+   * @param {Function} cf                The f(Aij..z,b) operation to invoke
+   * @param {boolean}  inverse           A true value indicates callback should be invoked f(b,Aij..z)
+   *
+   * @return {Array}                    Array (C)
+   *
+   * https://github.com/josdejong/mathjs/pull/346#issuecomment-97659042
+   */
+  return function matAlgo15xAs (a, b, cf, inverse) {
+    return inverse ? map(a, v => cf(b, v), true) : map(a, v => cf(v, b), true)
+  }
+})
diff --git a/src/type/matrix/utils/matrixAlgorithmSuite.js b/src/type/matrix/utils/matrixAlgorithmSuite.js
index 216935b20e..db02f6d38d 100644
--- a/src/type/matrix/utils/matrixAlgorithmSuite.js
+++ b/src/type/matrix/utils/matrixAlgorithmSuite.js
@@ -2,6 +2,7 @@ import { factory } from '../../../utils/factory.js'
 import { extend } from '../../../utils/object.js'
 import { createMatAlgo13xDD } from './matAlgo13xDD.js'
 import { createMatAlgo14xDs } from './matAlgo14xDs.js'
+import { createMatAlgo15xAs } from './matAlgo15xAs.js'
 import { broadcast } from './broadcast.js'
 
 const name = 'matrixAlgorithmSuite'
@@ -11,6 +12,7 @@ export const createMatrixAlgorithmSuite = /* #__PURE__ */ factory(
   name, dependencies, ({ typed, matrix }) => {
     const matAlgo13xDD = createMatAlgo13xDD({ typed })
     const matAlgo14xDs = createMatAlgo14xDs({ typed })
+    const matAlgo15xAs = createMatAlgo15xAs()
 
     /**
      * Return a signatures object with the usual boilerplate of
@@ -115,9 +117,9 @@ export const createMatrixAlgorithmSuite = /* #__PURE__ */ factory(
           matrixSignatures[scalar + ', DenseMatrix'] =
             (x, y) => matAlgo14xDs(y, x, elop, true)
           matrixSignatures['Array,' + scalar] =
-            (x, y) => matAlgo14xDs(matrix(x), y, elop, false).valueOf()
+            (x, y) => matAlgo15xAs(x, y, elop, false)
           matrixSignatures[scalar + ', Array'] =
-            (x, y) => matAlgo14xDs(matrix(y), x, elop, true).valueOf()
+            (x, y) => matAlgo15xAs(y, x, elop, true)
         } else {
           matrixSignatures['DenseMatrix,' + scalar] =
             typed.referToSelf(self => (x, y) => {
@@ -129,11 +131,11 @@ export const createMatrixAlgorithmSuite = /* #__PURE__ */ factory(
             })
           matrixSignatures['Array,' + scalar] =
             typed.referToSelf(self => (x, y) => {
-              return matAlgo14xDs(matrix(x), y, self, false).valueOf()
+              return matAlgo15xAs(x, y, self, false)
             })
           matrixSignatures[scalar + ', Array'] =
             typed.referToSelf(self => (x, y) => {
-              return matAlgo14xDs(matrix(y), x, self, true).valueOf()
+              return matAlgo15xAs(y, x, self, true)
             })
         }
       }
diff --git a/src/utils/optimizeCallback.js b/src/utils/optimizeCallback.js
index 5e6a0b94fb..04a6938977 100644
--- a/src/utils/optimizeCallback.js
+++ b/src/utils/optimizeCallback.js
@@ -1,38 +1,36 @@
 import typed from 'typed-function'
-import { get, arraySize } from './array.js'
 import { typeOf as _typeOf } from './is.js'
 
 /**
- * Simplifies a callback function by reducing its complexity and potentially improving its performance.
+ * Simplifies a callback function to be used in mapping functions like `map`, `forEach`, etc. It determines if the callback is unary and optimizes it accordingly.
  *
  * @param {Function} callback The original callback function to simplify.
  * @param {Array|Matrix} array The array that will be used with the callback function.
  * @param {string} name The name of the function that is using the callback.
  * @param {boolean} isUnary If true, the callback function is unary and will be optimized as such.
- * @returns {Function} Returns a simplified version of the callback function.
+ * @returns {Object} Returns an object with properties `isUnary` and `fn`.
  */
 export function optimizeCallback (callback, array, name, isUnary) {
+  const isMatrix = array && array.isMatrix
   if (typed.isTypedFunction(callback)) {
     let numberOfArguments
     if (isUnary) {
       numberOfArguments = 1
     } else {
-      const size = array.isMatrix ? array.size() : arraySize(array)
-
+      const firstIndex = isMatrix ? array.size().map(() => 0) : [0] // it only needs to be an array to be recognized
+      const firstValue = findFirst(array)
+      const isEmpty = isMatrix ? array.size().at(-1) === 0 : array.length === 0
       // Check the size of the last dimension to see if the array/matrix is empty
-      const isEmpty = size.length ? size[size.length - 1] === 0 : true
-      if (isEmpty) {
+      if (isEmpty || firstValue === undefined) {
         // don't optimize callbacks for empty arrays/matrix, as they will never be called
         // and in fact will throw an exception when we try to access the first element below
         return { isUnary, fn: callback }
       }
 
-      const firstIndex = size.map(() => 0)
-      const firstValue = array.isMatrix ? array.get(firstIndex) : get(array, firstIndex)
       numberOfArguments = _findNumberOfArgumentsTyped(callback, firstValue, firstIndex, array)
     }
     let fastCallback
-    if (array.isMatrix && (array.dataType !== 'mixed' && array.dataType !== undefined)) {
+    if (isMatrix && (array.dataType !== 'mixed' && array.dataType !== undefined)) {
       const singleSignature = _findSingleSignatureWithArity(callback, numberOfArguments)
       fastCallback = (singleSignature !== undefined) ? singleSignature : callback
     } else {
@@ -46,10 +44,35 @@ export function optimizeCallback (callback, array, name, isUnary) {
     }
     return { isUnary: false, fn: (...args) => _tryFunctionWithArgs(fastCallback, args, name, callback.name) }
   }
-  if (isUnary === undefined) {
+  if (isUnary) {
+    return { isUnary, fn: callback }
+  } else {
     return { isUnary: _findIfCallbackIsUnary(callback), fn: callback }
+  }
+}
+
+export function findFirst (array) {
+  if (array && array.isMatrix) {
+    if (array.size().at(-1) === 0) {
+      return undefined
+    }
+    const size = array.size()
+    const firstIndex = size.map(() => 0)
+    return array.get(firstIndex)
   } else {
-    return { isUnary, fn: callback }
+    return traverse(array)
+  }
+  function traverse (value) {
+    if (Array.isArray(value)) {
+      for (let i = 0; i < value.length; i++) {
+        const found = traverse(value[i])
+        if (found !== undefined) {
+          return found
+        }
+      }
+    } else {
+      return value
+    }
   }
 }
 
diff --git a/test/benchmark/broadcast.js b/test/benchmark/broadcast.js
new file mode 100644
index 0000000000..bb509c0b7e
--- /dev/null
+++ b/test/benchmark/broadcast.js
@@ -0,0 +1,33 @@
+import { Bench } from 'tinybench'
+import { matrix, add, subtract, random } from '../../lib/esm/index.js'
+import { formatTaskResult } from './utils/formatTaskResult.js'
+
+const array = random([500, 500], -10, 10)
+const genericMatrix = matrix(array)
+const numberMatrix = matrix(array, 'dense', 'number')
+
+// console.log('data', array)
+// console.log('abs(data)', abs(array))npm run
+
+const bench = new Bench({ time: 100, iterations: 100 })
+  .add('add(array, 1)', () => {
+    add(array, 1)
+  })
+  .add('add(matrix, 1)', () => {
+    add(genericMatrix, 1)
+  })
+  .add('add(numberMatrix, 1)', () => {
+    add(numberMatrix, 1)
+  })
+  .add('subtract(array, 1)', () => {
+    subtract(array, 1)
+  })
+  .add('subtract(matrix, 1)', () => {
+    subtract(genericMatrix, 1)
+  })
+  .add('subtract(numberMatrix, 1)', () => {
+    subtract(numberMatrix, 1)
+  })
+
+bench.addEventListener('cycle', (event) => console.log(formatTaskResult(bench, event.task)))
+await bench.run()
diff --git a/test/unit-tests/function/arithmetic/add.test.js b/test/unit-tests/function/arithmetic/add.test.js
index aa5c385bf3..1aab403f21 100644
--- a/test/unit-tests/function/arithmetic/add.test.js
+++ b/test/unit-tests/function/arithmetic/add.test.js
@@ -31,6 +31,11 @@ describe('add', function () {
       assert.deepStrictEqual(add([3, 4], 2), [5, 6])
     })
 
+    it('should add a scalar and a jagged array correctly', function () {
+      assert.deepStrictEqual(add(2, [[3, 4], 5]), [[5, 6], 7])
+      assert.deepStrictEqual(add([[3, 4], 5], 2), [[5, 6], 7])
+    })
+
     it('should add broadcastable arrays correctly', function () {
       const a2 = [1, 2]
       const a3 = [[3], [4]]
diff --git a/test/unit-tests/function/arithmetic/dotDivide.test.js b/test/unit-tests/function/arithmetic/dotDivide.test.js
index e99fa1cde9..feb4c35c94 100644
--- a/test/unit-tests/function/arithmetic/dotDivide.test.js
+++ b/test/unit-tests/function/arithmetic/dotDivide.test.js
@@ -68,6 +68,9 @@ describe('dotDivide', function () {
       const a = [[1, 2], [3, 4]]
       assert.deepStrictEqual(dotDivide(a, 2), [[0.5, 1], [1.5, 2]])
       assert.deepStrictEqual(dotDivide([], 2), [])
+      // jagged array
+      const ja = [[1, 2, 3], [4], [5, 6]]
+      assert.deepStrictEqual(dotDivide(ja, 2), [[0.5, 1, 1.5], [2], [2.5, 3]])
     })
 
     it('should divide 1 over a array element-wise', function () {
diff --git a/test/unit-tests/function/arithmetic/dotMultiply.test.js b/test/unit-tests/function/arithmetic/dotMultiply.test.js
index 8e7fb57e75..9468737962 100644
--- a/test/unit-tests/function/arithmetic/dotMultiply.test.js
+++ b/test/unit-tests/function/arithmetic/dotMultiply.test.js
@@ -76,6 +76,10 @@ describe('dotMultiply', function () {
       approxDeepEqual(dotMultiply(3, a), [[3, 0], [9, 12]])
       approxDeepEqual(dotMultiply([1, 2, 3, 4], 2), [2, 4, 6, 8])
       approxDeepEqual(dotMultiply(2, [1, 2, 3, 4]), [2, 4, 6, 8])
+      // jagged array
+      const ja = [[1, 2, 3], [4], [5, 6]]
+      approxDeepEqual(dotMultiply(ja, 2), [[2, 4, 6], [8], [10, 12]])
+      approxDeepEqual(dotMultiply(2, ja), [[2, 4, 6], [8], [10, 12]])
     })
 
     it('should multiply broadcastable arrays element-wise', function () {
diff --git a/test/unit-tests/function/arithmetic/dotPow.test.js b/test/unit-tests/function/arithmetic/dotPow.test.js
index 86bcff95c5..a543c14ec4 100644
--- a/test/unit-tests/function/arithmetic/dotPow.test.js
+++ b/test/unit-tests/function/arithmetic/dotPow.test.js
@@ -92,12 +92,16 @@ describe('dotPow', function () {
       approxDeepEqual(dotPow([[1, 2], [0, 4]], 2), [[1, 4], [0, 16]])
       approxDeepEqual(dotPow([[1, 2], [0, 4]], 2.5), [[1, 5.65685424949238], [0, 32]])
       approxDeepEqual(dotPow([[1, 2, 3], [4, 5, 0]], 2), [[1, 4, 9], [16, 25, 0]])
+      // jagged array
+      approxDeepEqual(dotPow([[1, 2, 3], [4], [5, 6]], 2), [[1, 4, 9], [16], [25, 36]])
     })
 
     it('should elevate scalar .^ array', function () {
       approxDeepEqual(dotPow(2, [[1, 2], [0, 4]]), [[2, 4], [1, 16]])
       approxDeepEqual(dotPow(2.5, [[1, 2], [0, 4]]), [[2.5, 6.25], [1, 39.0625]])
       approxDeepEqual(dotPow(2, [[1, 2, 3], [4, 5, 0]]), [[2, 4, 8], [16, 32, 1]])
+      // jagged array
+      approxDeepEqual(dotPow(2, [[1, 2, 3], [4], [5, 6]]), [[2, 4, 8], [16], [32, 64]])
     })
 
     it('should elevate broadcastable arrays element-wise', function () {
diff --git a/test/unit-tests/function/arithmetic/mod.test.js b/test/unit-tests/function/arithmetic/mod.test.js
index 5055a0f0b5..9906e3990e 100644
--- a/test/unit-tests/function/arithmetic/mod.test.js
+++ b/test/unit-tests/function/arithmetic/mod.test.js
@@ -155,6 +155,11 @@ describe('mod', function () {
       approxDeepEqual(mod(3, [[4, 3], [2, 1]]), [[3, 0], [1, 0]])
     })
 
+    it('should perform element-wise modulus on a jagged array and a scalar', function () {
+      approxDeepEqual(mod([[-4, -3, 0, -1], [0, 1, 2]], 3), [[2, 0, 0, 2], [0, 1, 2]])
+      approxDeepEqual(mod(3, [[4, 3, 2], [2, 1]]), [[3, 0, 1], [1, 0]])
+    })
+
     it('should perform element-wise modulus on broadcastable arrays', function () {
       approxDeepEqual(mod([-40, -31], [[3], [1]]), [[2, 2], [0, 0]])
       approxDeepEqual(mod([[-40], [-31]], [3, 1]), [[2, 0], [2, 0]])
diff --git a/test/unit-tests/function/arithmetic/subtract.test.js b/test/unit-tests/function/arithmetic/subtract.test.js
index 1302032171..5a8c1525f7 100644
--- a/test/unit-tests/function/arithmetic/subtract.test.js
+++ b/test/unit-tests/function/arithmetic/subtract.test.js
@@ -148,7 +148,12 @@ describe('subtract', function () {
       assert.deepStrictEqual(subtract([3, 0], 2), [1, -2])
     })
 
-    it('should substract broadcastable arrays correctly', function () {
+    it('should subtract a scalar and a jagged array correctly', function () {
+      assert.deepStrictEqual(subtract(2, [[3, 4], 5]), [[-1, -2], -3])
+      assert.deepStrictEqual(subtract([[3, 4], 5], 2), [[1, 2], 3])
+    })
+
+    it('should subtract broadcastable arrays correctly', function () {
       const a2 = [1, 2]
       const a3 = [[3], [4]]
       const a4 = subtract(a2, a3)
diff --git a/test/unit-tests/function/bitwise/bitAnd.test.js b/test/unit-tests/function/bitwise/bitAnd.test.js
index 9f1bee44e0..4fa95b61ac 100644
--- a/test/unit-tests/function/bitwise/bitAnd.test.js
+++ b/test/unit-tests/function/bitwise/bitAnd.test.js
@@ -199,6 +199,11 @@ describe('bitAnd', function () {
     assert.deepStrictEqual(bitAnd([3, 9], 12), [0, 8])
   })
 
+  it('should bitwise and a scalar and an jagged array correctly', function () {
+    assert.deepStrictEqual(bitAnd(12, [[3, 9], 5]), [[0, 8], 4])
+    assert.deepStrictEqual(bitAnd([[3, 9], 5], 12), [[0, 8], 4])
+  })
+
   it('should bitwise and broadcastable arrays correctly', function () {
     assert.deepStrictEqual(bitAnd([12, 13], [[3], [9]]), [[0, 1], [8, 9]])
     assert.deepStrictEqual(bitAnd([[12], [13]], [3, 9]), [[0, 8], [1, 9]])
diff --git a/test/unit-tests/function/bitwise/bitOr.test.js b/test/unit-tests/function/bitwise/bitOr.test.js
index d0a9aa27f7..c4a9a2bb82 100644
--- a/test/unit-tests/function/bitwise/bitOr.test.js
+++ b/test/unit-tests/function/bitwise/bitOr.test.js
@@ -199,6 +199,11 @@ describe('bitOr', function () {
     assert.deepStrictEqual(bitOr([3, 9], 12), [15, 13])
   })
 
+  it('should bitwise or a scalar and a jagged array correctly', function () {
+    assert.deepStrictEqual(bitOr(12, [[3, 9], 5]), [[15, 13], 13])
+    assert.deepStrictEqual(bitOr([[3, 9], 5], 12), [[15, 13], 13])
+  })
+
   it('should bitwise or broadcastable arrays correctly', function () {
     const a = [6, 4, 28]
     const b = [[13], [92], [101]]
diff --git a/test/unit-tests/function/bitwise/leftShift.test.js b/test/unit-tests/function/bitwise/leftShift.test.js
index 321d701438..f7c8b66794 100644
--- a/test/unit-tests/function/bitwise/leftShift.test.js
+++ b/test/unit-tests/function/bitwise/leftShift.test.js
@@ -111,6 +111,11 @@ describe('leftShift', function () {
       assert.deepStrictEqual(leftShift(2, [[1, 2], [8, 0]]), [[4, 8], [512, 2]])
     })
 
+    it('should left shift a jagged array and scalar', function () {
+      assert.deepStrictEqual(leftShift([[1, 2], [8, 0, 2]], 2), [[4, 8], [32, 0, 8]])
+      assert.deepStrictEqual(leftShift(2, [[1, 2], [8, 0, 8]]), [[4, 8], [512, 2, 512]])
+    })
+
     it('should left shift array - array', function () {
       assert.deepStrictEqual(leftShift([[1, 2], [8, 0]], [[4, 8], [32, 0]]), [[16, 512], [8, 0]])
       assert.deepStrictEqual(leftShift([[4, 8], [32, 0]], [[1, 2], [8, 0]]), [[8, 32], [8192, 0]])
diff --git a/test/unit-tests/function/matrix/map.test.js b/test/unit-tests/function/matrix/map.test.js
index 4a864c23b5..ea4d2b8f21 100644
--- a/test/unit-tests/function/matrix/map.test.js
+++ b/test/unit-tests/function/matrix/map.test.js
@@ -16,6 +16,13 @@ describe('map', function () {
     assert.ok(Array.isArray(arr2))
   })
 
+  it('should map jagged arrays', function () {
+    const map = math.map
+    const sqrt = math.square
+    assert.deepStrictEqual(map([[2], []], sqrt), [[4], []])
+    assert.deepStrictEqual(map([[], [2]], sqrt), [[], [4]])
+  })
+
   it('should map two arrays', function () {
     const arrA = [[1, 2, 3], [4, 5, 6]]
     const arrB = [[10, 20, 30], [40, 50, 60]]
diff --git a/test/unit-tests/utils/optimizeCallback.test.js b/test/unit-tests/utils/optimizeCallback.test.js
new file mode 100644
index 0000000000..4e65679807
--- /dev/null
+++ b/test/unit-tests/utils/optimizeCallback.test.js
@@ -0,0 +1,124 @@
+import assert from 'assert'
+import math from '../../../src/defaultInstance.js'
+import { optimizeCallback, findFirst } from '../../../src/utils/optimizeCallback.js'
+
+describe('utils.optimizeCallback', function () {
+  describe('optimizeCallback', function () {
+    function unaryCallback (a) {
+      return a + 1
+    }
+    const typedUnaryCallback = math.typed('unaryCallback', {
+      number: unaryCallback
+    })
+
+    function ternaryCallback (a, index, c) {
+      return a + index[0] + c[0]
+    }
+    const typedTernaryCallback = math.typed('ternaryCallback', {
+      'number, Array, Array': ternaryCallback
+    })
+
+    const arrayOfNumbers = [1, 2, 3]
+    const matrixOfNumbers = math.matrix(arrayOfNumbers)
+    const sparseMatrixOfNumbers = math.matrix(arrayOfNumbers, 'sparse')
+    const name = 'myFunction'
+
+    it('should find the right number of arguments for a function', function () {
+      const optimizedUnary = optimizeCallback(unaryCallback, arrayOfNumbers, name)
+      assert.strictEqual(optimizedUnary.isUnary, true)
+      const optimizedTernary = optimizeCallback(ternaryCallback, arrayOfNumbers, name)
+      assert.strictEqual(optimizedTernary.isUnary, false)
+      const optimizedUnaryMatrix = optimizeCallback(unaryCallback, matrixOfNumbers, name)
+      assert.strictEqual(optimizedUnaryMatrix.isUnary, true)
+      const optimizedTernaryMatrix = optimizeCallback(ternaryCallback, matrixOfNumbers, name)
+      assert.strictEqual(optimizedTernaryMatrix.isUnary, false)
+      const optimizedUnarySparseMatrix = optimizeCallback(unaryCallback, sparseMatrixOfNumbers, name)
+      assert.strictEqual(optimizedUnarySparseMatrix.isUnary, true)
+      const optimizedTernarySparseMatrix = optimizeCallback(ternaryCallback, sparseMatrixOfNumbers, name)
+      assert.strictEqual(optimizedTernarySparseMatrix.isUnary, false)
+    })
+
+    it('should find the right number of arguments for a typed function', function () {
+      const optimizedTypedUnary = optimizeCallback(typedUnaryCallback, arrayOfNumbers, name)
+      assert.strictEqual(optimizedTypedUnary.isUnary, true)
+      const optimizedTypedTernary = optimizeCallback(typedTernaryCallback, arrayOfNumbers, name)
+      assert.strictEqual(optimizedTypedTernary.isUnary, false)
+      const optimizedTypedUnaryMatrix = optimizeCallback(typedUnaryCallback, matrixOfNumbers, name)
+      assert.strictEqual(optimizedTypedUnaryMatrix.isUnary, true)
+      const optimizedTypedTernaryMatrix = optimizeCallback(typedTernaryCallback, matrixOfNumbers, name)
+      assert.strictEqual(optimizedTypedTernaryMatrix.isUnary, false)
+      const optimizedTypedUnarySparseMatrix = optimizeCallback(typedUnaryCallback, sparseMatrixOfNumbers, name)
+      assert.strictEqual(optimizedTypedUnarySparseMatrix.isUnary, true)
+      const optimizedTypedTernarySparseMatrix = optimizeCallback(typedTernaryCallback, sparseMatrixOfNumbers, name)
+      assert.strictEqual(optimizedTypedTernarySparseMatrix.isUnary, false)
+    })
+
+    it('should return unary functions as unary when isUnary is set to true', function () {
+      const optimizedUnary = optimizeCallback(unaryCallback, arrayOfNumbers, name, true)
+      assert.strictEqual(optimizedUnary.isUnary, true)
+      const optimizedTernary = optimizeCallback(ternaryCallback, arrayOfNumbers, name, true)
+      assert.strictEqual(optimizedTernary.isUnary, true)
+      const optimizedTypedUnary = optimizeCallback(typedUnaryCallback, arrayOfNumbers, name, true)
+      assert.strictEqual(optimizedTypedUnary.isUnary, true)
+      const optimizedTypedTernary = optimizeCallback(typedTernaryCallback, arrayOfNumbers, name, true)
+      assert.strictEqual(optimizedTypedTernary.isUnary, true)
+    })
+
+    it('should run the optimized callback', function () {
+      const optimizedUnary = optimizeCallback(unaryCallback, arrayOfNumbers, name)
+      const args = [1, [2, 2], [3]]
+      assert.strictEqual(unaryCallback(1), 2)
+      assert.strictEqual(optimizedUnary.fn(1), 2)
+      const optimizedTernary = optimizeCallback(ternaryCallback, arrayOfNumbers, name)
+      assert.strictEqual(optimizedTernary.fn(...args), 6)
+      const optimizedTypedUnary = optimizeCallback(typedUnaryCallback, arrayOfNumbers, name)
+      assert.strictEqual(optimizedTypedUnary.fn(1), 2)
+      const optimizedTypedTernary = optimizeCallback(typedTernaryCallback, arrayOfNumbers, name)
+      assert.strictEqual(optimizedTypedTernary.fn(...args), 6)
+      const optimizedUnaryMatrix = optimizeCallback(unaryCallback, matrixOfNumbers, name)
+      assert.strictEqual(optimizedUnaryMatrix.fn(1), 2)
+      const optimizedTernaryMatrix = optimizeCallback(ternaryCallback, matrixOfNumbers, name)
+      assert.strictEqual(optimizedTernaryMatrix.fn(...args), 6)
+      const optimizedTypedUnaryMatrix = optimizeCallback(typedUnaryCallback, matrixOfNumbers, name)
+      assert.strictEqual(optimizedTypedUnaryMatrix.fn(1), 2)
+      const optimizedTypedTernaryMatrix = optimizeCallback(typedTernaryCallback, matrixOfNumbers, name)
+      assert.strictEqual(optimizedTypedTernaryMatrix.fn(...args), 6)
+      const optimizedUnarySparseMatrix = optimizeCallback(unaryCallback, sparseMatrixOfNumbers, name)
+      assert.strictEqual(optimizedUnarySparseMatrix.fn(1), 2)
+      const optimizedTernarySparseMatrix = optimizeCallback(ternaryCallback, sparseMatrixOfNumbers, name)
+      assert.strictEqual(optimizedTernarySparseMatrix.fn(...args), 6)
+    })
+  }
+  )
+
+  describe('findFirst', function () {
+    it('should find the first value in a nested array', function () {
+      assert.strictEqual(findFirst([[[1, 2], [3, 4]], [[5, 6], [7, 8]]]), 1)
+      assert.strictEqual(findFirst([[[[1], [2]], [[3], [4]]], [[[5], [6]], [[7], [8]]]]), 1)
+      assert.strictEqual(findFirst([[[[1, 2], [3, 4]], [[5, 6], [7, 8]]], [[[9, 10], [11, 12]], [[13, 14], [15, 16]]]]), 1)
+    })
+
+    it('should find the first value in a matrix', function () {
+      assert.strictEqual(findFirst(math.matrix([[[1, 2], [3, 4]], [[5, 6], [7, 8]]])), 1)
+      assert.strictEqual(findFirst(math.matrix([[[[1], [2]], [[3], [4]]], [[[5], [6]], [[7], [8]]]])), 1)
+      assert.strictEqual(findFirst(math.matrix([[[[1, 2], [3, 4]], [[5, 6], [7, 8]]], [[[9, 10], [11, 12]], [[13, 14], [15, 16]]]])), 1)
+    })
+
+    it('should find the first value in a sparse matrix', function () {
+      assert.strictEqual(findFirst(math.matrix([[1, 2], [3, 4]], 'sparse')), 1)
+      assert.strictEqual(findFirst(math.matrix([[1], [2]], 'sparse')), 1)
+      assert.strictEqual(findFirst(math.matrix([[1, 2], [3, 4], [5, 6]], 'sparse')), 1)
+    })
+
+    it('should return undefined for an empty array or matrix', function () {
+      assert.strictEqual(findFirst([]), undefined)
+      assert.strictEqual(findFirst(math.matrix([])), undefined)
+    })
+
+    it('should find the first value in a jagged array if the first element is empty', function () {
+      assert.strictEqual(findFirst([[], [1, 2], [3, 4]]), 1)
+      assert.strictEqual(findFirst([[], [[1], [2]], [[3], [4]]]), 1)
+      assert.strictEqual(findFirst([[], [[[1, 2], [3, 4]], [[5, 6], [7, 8]]], [[[9, 10], [11, 12]], [[13, 14], [15, 16]]]]), 1)
+    })
+  })
+})