Added Vectors and Matrics operations, slicing, docs.

lynglib/src/commonMain/kotlin/net/sergeych/lyng/Compiler.kt

@@ -2915,9 +2915,9 @@ class Compiler(
                     operand?.let { left ->
                         // array access via ObjRef
                         val isOptional = t.type == Token.Type.NULL_COALESCE_INDEX
-                        val index = parseStatement() ?: throw ScriptError(t.pos, "Expecting index expression")
+                        val index = parseIndexExpression() ?: throw ScriptError(t.pos, "Expecting index expression")
                         cc.skipTokenOfType(Token.Type.RBRACKET, "missing ']' at the end of the list literal")
-                        operand = IndexRef(left, StatementRef(index), isOptional)
+                        operand = IndexRef(left, index, isOptional)
                     } ?: run {
                         // array literal
                         val entries = parseArrayLiteral()
@@ -3412,6 +3412,20 @@ class Compiler(
         }
     }
 
+    private suspend fun parseIndexExpression(): ObjRef? {
+        val first = parseExpressionLevel() ?: return null
+        if (!cc.skipTokenOfType(Token.Type.COMMA, isOptional = true)) {
+            return first
+        }
+
+        val entries = mutableListOf<ListEntry>(ListEntry.Element(first))
+        do {
+            val next = parseExpressionLevel() ?: throw ScriptError(cc.currentPos(), "Expecting index expression")
+            entries += ListEntry.Element(next)
+        } while (cc.skipTokenOfType(Token.Type.COMMA, isOptional = true))
+        return ListLiteralRef(entries)
+    }
+
     private suspend fun parseDestructuringPattern(): List<ListEntry> {
         // it should be called after Token.Type.LBRACKET is consumed
         val entries = mutableListOf<ListEntry>()

lynglib/src/commonMain/kotlin/net/sergeych/lyng/Script.kt

@@ -29,6 +29,7 @@ import net.sergeych.lyng.obj.*
 import net.sergeych.lyng.pacman.ImportManager
 import net.sergeych.lyng.stdlib_included.complexLyng
 import net.sergeych.lyng.stdlib_included.decimalLyng
+import net.sergeych.lyng.stdlib_included.matrixLyng
 import net.sergeych.lyng.stdlib_included.observableLyng
 import net.sergeych.lyng.stdlib_included.operatorsLyng
 import net.sergeych.lyng.stdlib_included.rootLyng
@@ -839,6 +840,10 @@ class Script(
                     module.eval(Source("lyng.decimal", decimalLyng))
                     ObjBigDecimalSupport.bindTo(module)
                 }
+                addPackage("lyng.matrix") { module ->
+                    module.eval(Source("lyng.matrix", matrixLyng))
+                    ObjMatrixSupport.bindTo(module)
+                }
                 addPackage("lyng.complex") { module ->
                     module.eval(Source("lyng.complex", complexLyng))
                 }

lynglib/src/commonMain/kotlin/net/sergeych/lyng/matrix/MatrixData.kt

@@ -0,0 +1,152 @@
+/*
+ * Copyright 2026 Sergey S. Chernov real.sergeych@gmail.com
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+
+package net.sergeych.lyng.matrix
+
+internal data class MatrixData(
+    val rows: Int,
+    val cols: Int,
+    val values: DoubleArray
+) {
+    init {
+        require(rows > 0) { "matrix must have at least one row" }
+        require(cols > 0) { "matrix must have at least one column" }
+        require(values.size == rows * cols) {
+            "matrix data size ${values.size} does not match shape ${rows}x${cols}"
+        }
+    }
+
+    val isSquare: Boolean get() = rows == cols
+
+    fun at(row: Int, col: Int): Double {
+        require(row in 0 until rows) { "row index $row out of bounds for $rows rows" }
+        require(col in 0 until cols) { "column index $col out of bounds for $cols columns" }
+        return values[row * cols + col]
+    }
+
+    fun rowValues(row: Int): DoubleArray {
+        require(row in 0 until rows) { "row index $row out of bounds for $rows rows" }
+        val start = row * cols
+        return values.copyOfRange(start, start + cols)
+    }
+
+    fun columnValues(col: Int): DoubleArray {
+        require(col in 0 until cols) { "column index $col out of bounds for $cols columns" }
+        return DoubleArray(rows) { row -> values[row * cols + col] }
+    }
+
+    fun slice(rowIndices: IntArray, colIndices: IntArray): MatrixData {
+        require(rowIndices.isNotEmpty()) { "matrix slice must include at least one row" }
+        require(colIndices.isNotEmpty()) { "matrix slice must include at least one column" }
+        val out = DoubleArray(rowIndices.size * colIndices.size)
+        var offset = 0
+        for (row in rowIndices) {
+            require(row in 0 until rows) { "row index $row out of bounds for $rows rows" }
+            val rowOffset = row * cols
+            for (col in colIndices) {
+                require(col in 0 until cols) { "column index $col out of bounds for $cols columns" }
+                out[offset++] = values[rowOffset + col]
+            }
+        }
+        return MatrixData(rowIndices.size, colIndices.size, out)
+    }
+
+    fun plus(other: MatrixData): MatrixData {
+        requireSameShape(other)
+        return MatrixData(rows, cols, DoubleArray(values.size) { index -> values[index] + other.values[index] })
+    }
+
+    fun minus(other: MatrixData): MatrixData {
+        requireSameShape(other)
+        return MatrixData(rows, cols, DoubleArray(values.size) { index -> values[index] - other.values[index] })
+    }
+
+    fun scale(factor: Double): MatrixData =
+        MatrixData(rows, cols, DoubleArray(values.size) { index -> values[index] * factor })
+
+    fun divide(divisor: Double): MatrixData {
+        require(divisor != 0.0) { "matrix division by zero" }
+        return MatrixData(rows, cols, DoubleArray(values.size) { index -> values[index] / divisor })
+    }
+
+    fun transpose(): MatrixData {
+        val out = DoubleArray(values.size)
+        for (row in 0 until rows) {
+            for (col in 0 until cols) {
+                out[col * rows + row] = values[row * cols + col]
+            }
+        }
+        return MatrixData(cols, rows, out)
+    }
+
+    fun multiply(other: MatrixData): MatrixData = PlatformMatrixBackend.multiply(this, other)
+
+    fun multiply(other: VectorData): VectorData = PureMatrixAlgorithms.multiply(this, other)
+
+    fun solve(other: VectorData): VectorData = PureMatrixAlgorithms.solve(this, other)
+
+    fun solve(other: MatrixData): MatrixData = PureMatrixAlgorithms.solve(this, other)
+
+    fun determinant(): Double = PureMatrixAlgorithms.determinant(this)
+
+    fun inverse(): MatrixData = PlatformMatrixBackend.inverse(this)
+
+    fun trace(): Double = PureMatrixAlgorithms.trace(this)
+
+    fun rank(): Int = PureMatrixAlgorithms.rank(this)
+
+    fun toNestedLists(): List<List<Double>> =
+        List(rows) { row ->
+            List(cols) { col -> values[row * cols + col] }
+        }
+
+    fun render(): String = buildString {
+        append("Matrix(")
+        append(rows)
+        append("x")
+        append(cols)
+        append(", [")
+        for (row in 0 until rows) {
+            if (row > 0) append(", ")
+            append("[")
+            for (col in 0 until cols) {
+                if (col > 0) append(", ")
+                append(formatMatrixValue(values[row * cols + col]))
+            }
+            append("]")
+        }
+        append("])")
+    }
+
+    fun compareTo(other: MatrixData): Int {
+        val rowCmp = rows.compareTo(other.rows)
+        if (rowCmp != 0) return rowCmp
+        val colCmp = cols.compareTo(other.cols)
+        if (colCmp != 0) return colCmp
+        for (index in values.indices) {
+            val cmp = values[index].compareTo(other.values[index])
+            if (cmp != 0) return cmp
+        }
+        return 0
+    }
+
+    private fun requireSameShape(other: MatrixData) {
+        require(rows == other.rows && cols == other.cols) {
+            "matrix shape mismatch: ${rows}x${cols} vs ${other.rows}x${other.cols}"
+        }
+    }
+}

lynglib/src/commonMain/kotlin/net/sergeych/lyng/matrix/MatrixFormatting.kt

@@ -0,0 +1,30 @@
+/*
+ * Copyright 2026 Sergey S. Chernov real.sergeych@gmail.com
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+
+package net.sergeych.lyng.matrix
+
+import net.sergeych.lyng.obj.ObjReal
+
+internal fun formatMatrixValue(value: Double): String {
+    if (value.isFinite()) {
+        val asLong = value.toLong()
+        if (asLong.toDouble() == value) {
+            return asLong.toString()
+        }
+    }
+    return ObjReal.of(value).toString()
+}

lynglib/src/commonMain/kotlin/net/sergeych/lyng/matrix/PlatformMatrixBackend.kt

@@ -0,0 +1,23 @@
+/*
+ * Copyright 2026 Sergey S. Chernov real.sergeych@gmail.com
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+
+package net.sergeych.lyng.matrix
+
+internal expect object PlatformMatrixBackend {
+    fun multiply(left: MatrixData, right: MatrixData): MatrixData
+    fun inverse(matrix: MatrixData): MatrixData
+}

lynglib/src/commonMain/kotlin/net/sergeych/lyng/matrix/PureMatrixAlgorithms.kt

@@ -0,0 +1,286 @@
+/*
+ * Copyright 2026 Sergey S. Chernov real.sergeych@gmail.com
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+
+package net.sergeych.lyng.matrix
+
+internal object PureMatrixAlgorithms {
+    private const val singularEpsilon = 1e-12
+
+    fun multiply(left: MatrixData, right: MatrixData): MatrixData {
+        require(left.cols == right.rows) {
+            "matrix multiplication shape mismatch: ${left.rows}x${left.cols} cannot multiply ${right.rows}x${right.cols}"
+        }
+        val out = DoubleArray(left.rows * right.cols)
+        for (row in 0 until left.rows) {
+            val leftRowOffset = row * left.cols
+            val outRowOffset = row * right.cols
+            for (pivot in 0 until left.cols) {
+                val leftValue = left.values[leftRowOffset + pivot]
+                val rightRowOffset = pivot * right.cols
+                for (col in 0 until right.cols) {
+                    out[outRowOffset + col] += leftValue * right.values[rightRowOffset + col]
+                }
+            }
+        }
+        return MatrixData(left.rows, right.cols, out)
+    }
+
+    fun multiply(left: MatrixData, right: VectorData): VectorData {
+        require(left.cols == right.size) {
+            "matrix-vector multiplication shape mismatch: ${left.rows}x${left.cols} cannot multiply length ${right.size}"
+        }
+        val out = DoubleArray(left.rows)
+        for (row in 0 until left.rows) {
+            var sum = 0.0
+            val rowOffset = row * left.cols
+            for (col in 0 until left.cols) {
+                sum += left.values[rowOffset + col] * right.values[col]
+            }
+            out[row] = sum
+        }
+        return VectorData(out)
+    }
+
+    fun inverse(matrix: MatrixData): MatrixData {
+        require(matrix.isSquare) { "matrix inverse requires a square matrix, got ${matrix.rows}x${matrix.cols}" }
+        val n = matrix.rows
+        val width = n * 2
+        val augmented = DoubleArray(n * width)
+
+        for (row in 0 until n) {
+            for (col in 0 until n) {
+                augmented[row * width + col] = matrix.values[row * n + col]
+            }
+            augmented[row * width + (n + row)] = 1.0
+        }
+
+        for (pivotCol in 0 until n) {
+            var pivotRow = pivotCol
+            var pivotAbs = kotlin.math.abs(augmented[pivotRow * width + pivotCol])
+            for (candidate in pivotCol + 1 until n) {
+                val candidateAbs = kotlin.math.abs(augmented[candidate * width + pivotCol])
+                if (candidateAbs > pivotAbs) {
+                    pivotAbs = candidateAbs
+                    pivotRow = candidate
+                }
+            }
+            require(pivotAbs > singularEpsilon) { "matrix is singular and cannot be inverted" }
+            if (pivotRow != pivotCol) {
+                swapRows(augmented, width, pivotRow, pivotCol)
+            }
+
+            val pivotValue = augmented[pivotCol * width + pivotCol]
+            val pivotOffset = pivotCol * width
+            for (col in 0 until width) {
+                augmented[pivotOffset + col] /= pivotValue
+            }
+
+            for (row in 0 until n) {
+                if (row == pivotCol) continue
+                val factor = augmented[row * width + pivotCol]
+                if (factor == 0.0) continue
+                val rowOffset = row * width
+                for (col in 0 until width) {
+                    augmented[rowOffset + col] -= factor * augmented[pivotOffset + col]
+                }
+            }
+        }
+
+        val inverse = DoubleArray(n * n)
+        for (row in 0 until n) {
+            for (col in 0 until n) {
+                inverse[row * n + col] = augmented[row * width + n + col]
+            }
+        }
+        return MatrixData(n, n, inverse)
+    }
+
+    fun determinant(matrix: MatrixData): Double {
+        require(matrix.isSquare) { "matrix determinant requires a square matrix, got ${matrix.rows}x${matrix.cols}" }
+        val n = matrix.rows
+        val work = matrix.values.copyOf()
+        var sign = 1.0
+        var determinant = 1.0
+
+        for (pivotCol in 0 until n) {
+            var pivotRow = pivotCol
+            var pivotAbs = kotlin.math.abs(work[pivotRow * n + pivotCol])
+            for (candidate in pivotCol + 1 until n) {
+                val candidateAbs = kotlin.math.abs(work[candidate * n + pivotCol])
+                if (candidateAbs > pivotAbs) {
+                    pivotAbs = candidateAbs
+                    pivotRow = candidate
+                }
+            }
+            if (pivotAbs <= singularEpsilon) {
+                return 0.0
+            }
+            if (pivotRow != pivotCol) {
+                swapRows(work, n, pivotRow, pivotCol)
+                sign = -sign
+            }
+
+            val pivotValue = work[pivotCol * n + pivotCol]
+            determinant *= pivotValue
+
+            for (row in pivotCol + 1 until n) {
+                val rowOffset = row * n
+                val factor = work[rowOffset + pivotCol] / pivotValue
+                if (factor == 0.0) continue
+                for (col in pivotCol + 1 until n) {
+                    work[rowOffset + col] -= factor * work[pivotCol * n + col]
+                }
+            }
+        }
+
+        return determinant * sign
+    }
+
+    fun trace(matrix: MatrixData): Double {
+        require(matrix.isSquare) { "matrix trace requires a square matrix, got ${matrix.rows}x${matrix.cols}" }
+        var sum = 0.0
+        for (index in 0 until matrix.rows) {
+            sum += matrix.values[index * matrix.cols + index]
+        }
+        return sum
+    }
+
+    fun rank(matrix: MatrixData): Int {
+        val work = matrix.values.copyOf()
+        var rank = 0
+        var pivotRow = 0
+        val rows = matrix.rows
+        val cols = matrix.cols
+
+        for (pivotCol in 0 until cols) {
+            var bestRow = -1
+            var bestAbs = singularEpsilon
+            for (candidate in pivotRow until rows) {
+                val absValue = kotlin.math.abs(work[candidate * cols + pivotCol])
+                if (absValue > bestAbs) {
+                    bestAbs = absValue
+                    bestRow = candidate
+                }
+            }
+            if (bestRow == -1) continue
+            if (bestRow != pivotRow) {
+                swapRows(work, cols, bestRow, pivotRow)
+            }
+            val pivotValue = work[pivotRow * cols + pivotCol]
+            for (row in pivotRow + 1 until rows) {
+                val factor = work[row * cols + pivotCol] / pivotValue
+                if (kotlin.math.abs(factor) <= singularEpsilon) continue
+                val rowOffset = row * cols
+                val pivotOffset = pivotRow * cols
+                for (col in pivotCol until cols) {
+                    work[rowOffset + col] -= factor * work[pivotOffset + col]
+                }
+            }
+            rank += 1
+            pivotRow += 1
+            if (pivotRow == rows) break
+        }
+        return rank
+    }
+
+    fun solve(matrix: MatrixData, rhs: VectorData): VectorData {
+        require(matrix.isSquare) { "matrix solve requires a square matrix, got ${matrix.rows}x${matrix.cols}" }
+        require(matrix.rows == rhs.size) {
+            "matrix solve shape mismatch: ${matrix.rows}x${matrix.cols} cannot solve length ${rhs.size}"
+        }
+        val solution = solveAugmented(matrix.rows, matrix.cols, 1, matrix.values, rhs.values)
+        return VectorData(solution)
+    }
+
+    fun solve(matrix: MatrixData, rhs: MatrixData): MatrixData {
+        require(matrix.isSquare) { "matrix solve requires a square matrix, got ${matrix.rows}x${matrix.cols}" }
+        require(matrix.rows == rhs.rows) {
+            "matrix solve shape mismatch: ${matrix.rows}x${matrix.cols} cannot solve ${rhs.rows}x${rhs.cols}"
+        }
+        val solution = solveAugmented(matrix.rows, matrix.cols, rhs.cols, matrix.values, rhs.values)
+        return MatrixData(matrix.rows, rhs.cols, solution)
+    }
+
+    private fun solveAugmented(rows: Int, cols: Int, rhsCols: Int, matrixValues: DoubleArray, rhsValues: DoubleArray): DoubleArray {
+        val width = cols + rhsCols
+        val augmented = DoubleArray(rows * width)
+
+        for (row in 0 until rows) {
+            val matrixOffset = row * cols
+            val augmentedOffset = row * width
+            for (col in 0 until cols) {
+                augmented[augmentedOffset + col] = matrixValues[matrixOffset + col]
+            }
+            val rhsOffset = row * rhsCols
+            for (col in 0 until rhsCols) {
+                augmented[augmentedOffset + cols + col] = rhsValues[rhsOffset + col]
+            }
+        }
+
+        for (pivotCol in 0 until cols) {
+            var pivotRow = pivotCol
+            var pivotAbs = kotlin.math.abs(augmented[pivotRow * width + pivotCol])
+            for (candidate in pivotCol + 1 until rows) {
+                val candidateAbs = kotlin.math.abs(augmented[candidate * width + pivotCol])
+                if (candidateAbs > pivotAbs) {
+                    pivotAbs = candidateAbs
+                    pivotRow = candidate
+                }
+            }
+            require(pivotAbs > singularEpsilon) { "matrix is singular and cannot be solved" }
+            if (pivotRow != pivotCol) {
+                swapRows(augmented, width, pivotRow, pivotCol)
+            }
+
+            val pivotOffset = pivotCol * width
+            val pivotValue = augmented[pivotOffset + pivotCol]
+            for (col in pivotCol until width) {
+                augmented[pivotOffset + col] /= pivotValue
+            }
+
+            for (row in 0 until rows) {
+                if (row == pivotCol) continue
+                val factor = augmented[row * width + pivotCol]
+                if (kotlin.math.abs(factor) <= singularEpsilon) continue
+                val rowOffset = row * width
+                for (col in pivotCol until width) {
+                    augmented[rowOffset + col] -= factor * augmented[pivotOffset + col]
+                }
+            }
+        }
+
+        val solution = DoubleArray(rows * rhsCols)
+        for (row in 0 until rows) {
+            val augmentedOffset = row * width + cols
+            val solutionOffset = row * rhsCols
+            for (col in 0 until rhsCols) {
+                solution[solutionOffset + col] = augmented[augmentedOffset + col]
+            }
+        }
+        return solution
+    }
+
+    private fun swapRows(data: DoubleArray, stride: Int, firstRow: Int, secondRow: Int) {
+        val firstOffset = firstRow * stride
+        val secondOffset = secondRow * stride
+        for (col in 0 until stride) {
+            val tmp = data[firstOffset + col]
+            data[firstOffset + col] = data[secondOffset + col]
+            data[secondOffset + col] = tmp
+        }
+    }
+}

lynglib/src/commonMain/kotlin/net/sergeych/lyng/matrix/VectorData.kt

@@ -0,0 +1,113 @@
+/*
+ * Copyright 2026 Sergey S. Chernov real.sergeych@gmail.com
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+
+package net.sergeych.lyng.matrix
+
+import kotlin.math.sqrt
+
+internal data class VectorData(
+    val values: DoubleArray
+) {
+    init {
+        require(values.isNotEmpty()) { "vector must have at least one element" }
+    }
+
+    val size: Int get() = values.size
+
+    fun at(index: Int): Double {
+        require(index in values.indices) { "vector index $index out of bounds for length $size" }
+        return values[index]
+    }
+
+    fun plus(other: VectorData): VectorData {
+        require(size == other.size) { "vector size mismatch: $size vs ${other.size}" }
+        return VectorData(DoubleArray(size) { index -> values[index] + other.values[index] })
+    }
+
+    fun minus(other: VectorData): VectorData {
+        require(size == other.size) { "vector size mismatch: $size vs ${other.size}" }
+        return VectorData(DoubleArray(size) { index -> values[index] - other.values[index] })
+    }
+
+    fun scale(factor: Double): VectorData =
+        VectorData(DoubleArray(size) { index -> values[index] * factor })
+
+    fun divide(divisor: Double): VectorData {
+        require(divisor != 0.0) { "vector division by zero" }
+        return VectorData(DoubleArray(size) { index -> values[index] / divisor })
+    }
+
+    fun dot(other: VectorData): Double {
+        require(size == other.size) { "vector size mismatch: $size vs ${other.size}" }
+        var sum = 0.0
+        for (index in values.indices) {
+            sum += values[index] * other.values[index]
+        }
+        return sum
+    }
+
+    fun norm(): Double = sqrt(dot(this))
+
+    fun normalize(): VectorData {
+        val length = norm()
+        require(length != 0.0) { "cannot normalize a zero vector" }
+        return divide(length)
+    }
+
+    fun cross(other: VectorData): VectorData {
+        require(size == 3 && other.size == 3) { "cross product requires two 3D vectors" }
+        return VectorData(
+            doubleArrayOf(
+                values[1] * other.values[2] - values[2] * other.values[1],
+                values[2] * other.values[0] - values[0] * other.values[2],
+                values[0] * other.values[1] - values[1] * other.values[0]
+            )
+        )
+    }
+
+    fun outer(other: VectorData): MatrixData {
+        val out = DoubleArray(size * other.size)
+        for (row in values.indices) {
+            val rowOffset = row * other.size
+            for (col in other.values.indices) {
+                out[rowOffset + col] = values[row] * other.values[col]
+            }
+        }
+        return MatrixData(size, other.size, out)
+    }
+
+    fun toList(): List<Double> = values.asList()
+
+    fun render(): String = buildString {
+        append("Vector([")
+        for (index in values.indices) {
+            if (index > 0) append(", ")
+            append(formatMatrixValue(values[index]))
+        }
+        append("])")
+    }
+
+    fun compareTo(other: VectorData): Int {
+        val sizeCmp = size.compareTo(other.size)
+        if (sizeCmp != 0) return sizeCmp
+        for (index in values.indices) {
+            val cmp = values[index].compareTo(other.values[index])
+            if (cmp != 0) return cmp
+        }
+        return 0
+    }
+}

lynglib/src/commonMain/kotlin/net/sergeych/lyng/obj/ObjBigDecimalSupport.kt

@@ -17,22 +17,14 @@
 
 package net.sergeych.lyng.obj
 
-import com.ionspin.kotlin.bignum.decimal.BigDecimal as IonBigDecimal
 import com.ionspin.kotlin.bignum.decimal.DecimalMode
 import com.ionspin.kotlin.bignum.decimal.RoundingMode
 import com.ionspin.kotlin.bignum.integer.BigInteger
-import net.sergeych.lyng.Arguments
-import net.sergeych.lyng.FrameSlotRef
-import net.sergeych.lyng.InteropOperator
-import net.sergeych.lyng.ModuleScope
-import net.sergeych.lyng.OperatorInteropRegistry
-import net.sergeych.lyng.RecordSlotRef
-import net.sergeych.lyng.Scope
-import net.sergeych.lyng.ScopeFacade
-import net.sergeych.lyng.TypeDecl
+import net.sergeych.lyng.*
 import net.sergeych.lyng.miniast.addPropertyDoc
 import net.sergeych.lyng.miniast.type
 import net.sergeych.lyng.requiredArg
+import com.ionspin.kotlin.bignum.decimal.BigDecimal as IonBigDecimal
 
 object ObjBigDecimalSupport {
     private const val decimalContextVar = "__lyng_decimal_context__"
@@ -88,6 +80,9 @@ object ObjBigDecimalSupport {
         decimalClass.addFn("toReal") {
             ObjReal.of(valueOf(thisObj).doubleValue(false))
         }
+        decimalClass.addFn("toString") {
+            ObjString(valueOf(thisObj).toStringExpanded())
+        }
         decimalClass.addFn("toStringExpanded") {
             ObjString(valueOf(thisObj).toStringExpanded())
         }

lynglib/src/commonMain/kotlin/net/sergeych/lyng/obj/ObjMatrixSupport.kt

@@ -0,0 +1,385 @@
+/*
+ * Copyright 2026 Sergey S. Chernov real.sergeych@gmail.com
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+
+package net.sergeych.lyng.obj
+
+import net.sergeych.lyng.ModuleScope
+import net.sergeych.lyng.Scope
+import net.sergeych.lyng.ScopeFacade
+import net.sergeych.lyng.matrix.MatrixData
+import net.sergeych.lyng.matrix.VectorData
+import net.sergeych.lyng.requiredArg
+
+object ObjMatrixSupport {
+    private sealed interface MatrixAxisIndex {
+        data class Single(val value: Int) : MatrixAxisIndex
+        data class Slice(val values: IntArray) : MatrixAxisIndex
+    }
+
+    private object BoundMarker
+    private val defaultMatrix = MatrixData(1, 1, doubleArrayOf(0.0))
+    private val defaultVector = VectorData(doubleArrayOf(0.0))
+
+    suspend fun bindTo(module: ModuleScope) {
+        val matrixClass = module.requireClass("Matrix")
+        val vectorClass = module.requireClass("Vector")
+        if (matrixClass.kotlinClassData === BoundMarker && vectorClass.kotlinClassData === BoundMarker) return
+
+        bindVectorClass(vectorClass, matrixClass)
+        bindMatrixClass(matrixClass, vectorClass)
+    }
+
+    private fun bindVectorClass(vectorClass: ObjClass, matrixClass: ObjClass) {
+        if (vectorClass.kotlinClassData === BoundMarker) return
+        vectorClass.kotlinClassData = BoundMarker
+        vectorClass.isAbstract = false
+
+        val hooks = vectorClass.bridgeInitHooks ?: mutableListOf<suspend (ScopeFacade, ObjInstance) -> Unit>().also {
+            vectorClass.bridgeInitHooks = it
+        }
+        hooks += { _, instance -> instance.kotlinInstanceData = defaultVector }
+
+        vectorClass.addProperty("size", getter = {
+            ObjInt.of(vectorOf(thisObj).size.toLong())
+        })
+        vectorClass.addProperty("length", getter = {
+            ObjInt.of(vectorOf(thisObj).size.toLong())
+        })
+        vectorClass.addFn("toList") {
+            Obj.from(vectorOf(thisObj).toList())
+        }
+        vectorClass.addFn("get") {
+            ObjReal.of(vectorOf(thisObj).at(requiredArg<ObjInt>(0).value.toInt()))
+        }
+        vectorClass.addFn("plus") {
+            newVector(vectorClass, vectorOf(thisObj).plus(coerceVectorArg(requireScope(), args.firstAndOnly())))
+        }
+        vectorClass.addFn("minus") {
+            newVector(vectorClass, vectorOf(thisObj).minus(coerceVectorArg(requireScope(), args.firstAndOnly())))
+        }
+        vectorClass.addFn("mul") {
+            newVector(vectorClass, vectorOf(thisObj).scale(coerceScalarArg(requireScope(), args.firstAndOnly())))
+        }
+        vectorClass.addFn("div") {
+            newVector(vectorClass, vectorOf(thisObj).divide(coerceScalarArg(requireScope(), args.firstAndOnly())))
+        }
+        vectorClass.addFn("dot") {
+            ObjReal.of(vectorOf(thisObj).dot(coerceVectorArg(requireScope(), args.firstAndOnly())))
+        }
+        vectorClass.addFn("norm") {
+            ObjReal.of(vectorOf(thisObj).norm())
+        }
+        vectorClass.addFn("normalize") {
+            newVector(vectorClass, vectorOf(thisObj).normalize())
+        }
+        vectorClass.addFn("cross") {
+            newVector(vectorClass, vectorOf(thisObj).cross(coerceVectorArg(requireScope(), args.firstAndOnly())))
+        }
+        vectorClass.addFn("outer") {
+            newMatrix(matrixClass, vectorOf(thisObj).outer(coerceVectorArg(requireScope(), args.firstAndOnly())))
+        }
+        vectorClass.addFn("toString") {
+            ObjString(vectorOf(thisObj).render())
+        }
+        vectorClass.addFn("compareTo") {
+            ObjInt.of(vectorOf(thisObj).compareTo(coerceVectorArg(requireScope(), args.firstAndOnly())).toLong())
+        }
+
+        vectorClass.addClassFn("fromList") {
+            newVector(vectorClass, parseVector(requireScope(), requiredArg(0)))
+        }
+        vectorClass.addClassFn("zeros") {
+            val size = requiredArg<ObjInt>(0).value.toInt()
+            if (size <= 0) requireScope().raiseIllegalArgument("vector size must be positive")
+            newVector(vectorClass, VectorData(DoubleArray(size)))
+        }
+    }
+
+    private fun bindMatrixClass(matrixClass: ObjClass, vectorClass: ObjClass) {
+        if (matrixClass.kotlinClassData === BoundMarker) return
+        matrixClass.kotlinClassData = BoundMarker
+        matrixClass.isAbstract = false
+
+        val hooks = matrixClass.bridgeInitHooks ?: mutableListOf<suspend (ScopeFacade, ObjInstance) -> Unit>().also {
+            matrixClass.bridgeInitHooks = it
+        }
+        hooks += { _, instance -> instance.kotlinInstanceData = defaultMatrix }
+
+        matrixClass.addProperty("rows", getter = {
+            ObjInt.of(matrixOf(thisObj).rows.toLong())
+        })
+        matrixClass.addProperty("cols", getter = {
+            ObjInt.of(matrixOf(thisObj).cols.toLong())
+        })
+        matrixClass.addProperty("shape", getter = {
+            ObjList(
+                mutableListOf(
+                    ObjInt.of(matrixOf(thisObj).rows.toLong()),
+                    ObjInt.of(matrixOf(thisObj).cols.toLong())
+                )
+            )
+        })
+        matrixClass.addProperty("isSquare", getter = {
+            matrixOf(thisObj).isSquare.toObj()
+        })
+
+        matrixClass.addFn("plus") {
+            newMatrix(matrixClass, matrixOf(thisObj).plus(coerceMatrixArg(requireScope(), args.firstAndOnly())))
+        }
+        matrixClass.addFn("minus") {
+            newMatrix(matrixClass, matrixOf(thisObj).minus(coerceMatrixArg(requireScope(), args.firstAndOnly())))
+        }
+        matrixClass.addFn("mul") {
+            when (val other = args.firstAndOnly()) {
+                is ObjInstance -> when (other.objClass.className) {
+                    "Matrix" -> newMatrix(matrixClass, matrixOf(thisObj).multiply(matrixOf(other)))
+                    "Vector" -> newVector(vectorClass, matrixOf(thisObj).multiply(vectorOf(other)))
+                    else -> newMatrix(matrixClass, matrixOf(thisObj).scale(coerceScalarArg(requireScope(), other)))
+                }
+                else -> newMatrix(matrixClass, matrixOf(thisObj).scale(coerceScalarArg(requireScope(), other)))
+            }
+        }
+        matrixClass.addFn("div") {
+            newMatrix(matrixClass, matrixOf(thisObj).divide(coerceScalarArg(requireScope(), args.firstAndOnly())))
+        }
+        matrixClass.addFn("transpose") {
+            newMatrix(matrixClass, matrixOf(thisObj).transpose())
+        }
+        matrixClass.addFn("trace") {
+            ObjReal.of(matrixOf(thisObj).trace())
+        }
+        matrixClass.addFn("rank") {
+            ObjInt.of(matrixOf(thisObj).rank().toLong())
+        }
+        matrixClass.addFn("determinant") {
+            ObjReal.of(matrixOf(thisObj).determinant())
+        }
+        matrixClass.addFn("inverse") {
+            newMatrix(matrixClass, matrixOf(thisObj).inverse())
+        }
+        matrixClass.addFn("solve") {
+            when (val rhs = args.firstAndOnly()) {
+                is ObjInstance -> when (rhs.objClass.className) {
+                    "Vector" -> newVector(vectorClass, matrixOf(thisObj).solve(vectorOf(rhs)))
+                    "Matrix" -> newMatrix(matrixClass, matrixOf(thisObj).solve(matrixOf(rhs)))
+                    else -> requireScope().raiseClassCastError("Matrix.solve expects Vector or Matrix")
+                }
+                else -> requireScope().raiseClassCastError("Matrix.solve expects Vector or Matrix")
+            }
+        }
+        matrixClass.addFn("get") {
+            ObjReal.of(
+                matrixOf(thisObj).at(
+                    requiredArg<ObjInt>(0).value.toInt(),
+                    requiredArg<ObjInt>(1).value.toInt()
+                )
+            )
+        }
+        matrixClass.addFn("getAt") {
+            resolveMatrixIndex(matrixClass, matrixOf(thisObj), args.firstAndOnly(), thisObj)
+        }
+        matrixClass.addFn("row") {
+            doubleArrayToObjList(matrixOf(thisObj).rowValues(requiredArg<ObjInt>(0).value.toInt()))
+        }
+        matrixClass.addFn("column") {
+            doubleArrayToObjList(matrixOf(thisObj).columnValues(requiredArg<ObjInt>(0).value.toInt()))
+        }
+        matrixClass.addFn("toList") {
+            Obj.from(matrixOf(thisObj).toNestedLists())
+        }
+        matrixClass.addFn("toString") {
+            ObjString(matrixOf(thisObj).render())
+        }
+        matrixClass.addFn("compareTo") {
+            ObjInt.of(matrixOf(thisObj).compareTo(coerceMatrixArg(requireScope(), args.firstAndOnly())).toLong())
+        }
+
+        matrixClass.addClassFn("fromRows") {
+            newMatrix(matrixClass, parseRows(requireScope(), requiredArg(0)))
+        }
+        matrixClass.addClassFn("zeros") {
+            val rows = requiredArg<ObjInt>(0).value.toInt()
+            val cols = requiredArg<ObjInt>(1).value.toInt()
+            if (rows <= 0) requireScope().raiseIllegalArgument("matrix must have at least one row")
+            if (cols <= 0) requireScope().raiseIllegalArgument("matrix must have at least one column")
+            newMatrix(matrixClass, MatrixData(rows, cols, DoubleArray(rows * cols)))
+        }
+        matrixClass.addClassFn("identity") {
+            val size = requiredArg<ObjInt>(0).value.toInt()
+            if (size <= 0) requireScope().raiseIllegalArgument("identity matrix size must be positive")
+            val values = DoubleArray(size * size)
+            for (index in 0 until size) {
+                values[index * size + index] = 1.0
+            }
+            newMatrix(matrixClass, MatrixData(size, size, values))
+        }
+    }
+
+    private fun matrixOf(obj: Obj): MatrixData {
+        val instance = obj as? ObjInstance ?: error("Matrix receiver must be an object instance")
+        return instance.kotlinInstanceData as? MatrixData ?: defaultMatrix
+    }
+
+    private fun vectorOf(obj: Obj): VectorData {
+        val instance = obj as? ObjInstance ?: error("Vector receiver must be an object instance")
+        return instance.kotlinInstanceData as? VectorData ?: defaultVector
+    }
+
+    private suspend fun ScopeFacade.newMatrix(matrixClass: ObjClass, value: MatrixData): ObjInstance {
+        val instance = call(matrixClass) as? ObjInstance
+            ?: raiseIllegalState("Matrix() did not return an object instance")
+        instance.kotlinInstanceData = value
+        return instance
+    }
+
+    private suspend fun ScopeFacade.newVector(vectorClass: ObjClass, value: VectorData): ObjInstance {
+        val instance = call(vectorClass) as? ObjInstance
+            ?: raiseIllegalState("Vector() did not return an object instance")
+        instance.kotlinInstanceData = value
+        return instance
+    }
+
+    private fun coerceMatrixArg(scope: Scope, value: Obj): MatrixData {
+        val instance = value as? ObjInstance
+            ?: scope.raiseClassCastError("expected Matrix, got ${value.objClass.className}")
+        if (instance.objClass.className != "Matrix") {
+            scope.raiseClassCastError("expected Matrix, got ${instance.objClass.className}")
+        }
+        return instance.kotlinInstanceData as? MatrixData ?: defaultMatrix
+    }
+
+    private fun coerceVectorArg(scope: Scope, value: Obj): VectorData {
+        val instance = value as? ObjInstance
+            ?: scope.raiseClassCastError("expected Vector, got ${value.objClass.className}")
+        if (instance.objClass.className != "Vector") {
+            scope.raiseClassCastError("expected Vector, got ${instance.objClass.className}")
+        }
+        return instance.kotlinInstanceData as? VectorData ?: defaultVector
+    }
+
+    private fun coerceScalarArg(scope: Scope, value: Obj): Double = try {
+        value.toDouble()
+    } catch (_: IllegalArgumentException) {
+        scope.raiseClassCastError("expected matrix scalar (Int or Real), got ${value.objClass.className}")
+    }
+
+    private suspend fun parseRows(scope: Scope, rowsObj: Obj): MatrixData {
+        val outer = asObjList(scope, rowsObj, "Matrix.fromRows expects a list of rows")
+        if (outer.list.isEmpty()) scope.raiseIllegalArgument("matrix must have at least one row")
+
+        val rows = outer.list.size
+        var cols = -1
+        val values = mutableListOf<Double>()
+
+        for (rowObj in outer.list) {
+            val row = asObjList(scope, rowObj, "Matrix rows must be lists")
+            if (cols == -1) {
+                cols = row.list.size
+                if (cols <= 0) scope.raiseIllegalArgument("matrix must have at least one column")
+            } else if (row.list.size != cols) {
+                scope.raiseIllegalArgument("matrix rows must all have the same length")
+            }
+            for (cell in row.list) {
+                values += coerceScalarArg(scope, cell)
+            }
+        }
+        return MatrixData(rows, cols, values.toDoubleArray())
+    }
+
+    private suspend fun parseVector(scope: Scope, valuesObj: Obj): VectorData {
+        val list = asObjList(scope, valuesObj, "Vector.fromList expects a list")
+        if (list.list.isEmpty()) scope.raiseIllegalArgument("vector must have at least one element")
+        return VectorData(list.list.map { coerceScalarArg(scope, it) }.toDoubleArray())
+    }
+
+    private suspend fun ScopeFacade.resolveMatrixIndex(
+        matrixClass: ObjClass,
+        matrix: MatrixData,
+        index: Obj,
+        receiver: Obj
+    ): Obj {
+        val tuple = asObjList(requireScope(), index, "Matrix index must be [row, col]")
+        if (tuple.list.size != 2) {
+            raiseIllegalArgument("Matrix index must contain exactly two selectors: [row, col]")
+        }
+
+        val rowIndex = decodeAxisIndex(requireScope(), tuple.list[0], matrix.rows, "row")
+        val colIndex = decodeAxisIndex(requireScope(), tuple.list[1], matrix.cols, "column")
+
+        return when {
+            rowIndex is MatrixAxisIndex.Single && colIndex is MatrixAxisIndex.Single ->
+                ObjReal.of(matrix.at(rowIndex.value, colIndex.value))
+
+            rowIndex is MatrixAxisIndex.Single && colIndex is MatrixAxisIndex.Slice ->
+                newMatrix(matrixClass, matrix.slice(intArrayOf(rowIndex.value), colIndex.values))
+
+            rowIndex is MatrixAxisIndex.Slice && colIndex is MatrixAxisIndex.Single ->
+                newMatrix(matrixClass, matrix.slice(rowIndex.values, intArrayOf(colIndex.value)))
+
+            rowIndex is MatrixAxisIndex.Slice && colIndex is MatrixAxisIndex.Slice ->
+                newMatrix(matrixClass, matrix.slice(rowIndex.values, colIndex.values))
+
+            else -> requireScope().raiseIllegalState("unreachable matrix index state for ${receiver.objClass.className}")
+        }
+    }
+
+    private suspend fun decodeAxisIndex(scope: Scope, index: Obj, size: Int, axisName: String): MatrixAxisIndex =
+        when (index) {
+            is ObjInt -> {
+                val value = index.value.toInt()
+                if (value !in 0 until size) {
+                    scope.raiseIllegalArgument("$axisName index $value out of bounds for length $size")
+                }
+                MatrixAxisIndex.Single(value)
+            }
+
+            is ObjRange -> {
+                if (index.hasExplicitStep) {
+                    scope.raiseIllegalArgument("Matrix slicing does not support stepped $axisName ranges")
+                }
+                val start = index.startInt(scope)
+                val endExclusive = index.exclusiveIntEnd(scope) ?: size
+                if (start !in 0..size) {
+                    scope.raiseIllegalArgument("$axisName slice start $start out of bounds for length $size")
+                }
+                if (endExclusive !in 0..size) {
+                    scope.raiseIllegalArgument("$axisName slice end $endExclusive out of bounds for length $size")
+                }
+                if (start > endExclusive) {
+                    scope.raiseIllegalArgument("$axisName slice start $start is after end $endExclusive")
+                }
+                if (start == endExclusive) {
+                    scope.raiseIllegalArgument("Matrix slice must include at least one $axisName")
+                }
+                MatrixAxisIndex.Slice(IntArray(endExclusive - start) { start + it })
+            }
+
+            else -> scope.raiseClassCastError("Matrix $axisName selector must be Int or Range, got ${index.objClass.className}")
+        }
+
+    private suspend fun asObjList(scope: Scope, value: Obj, message: String): ObjList = when (value) {
+        is ObjList -> value
+        else -> if (value.isInstanceOf(ObjIterable)) {
+            value.callMethod<ObjList>(scope, "toList")
+        } else {
+            scope.raiseClassCastError(message)
+        }
+    }
+
+    private fun doubleArrayToObjList(values: DoubleArray): ObjList =
+        ObjList(values.map { ObjReal.of(it) }.toMutableList())
+}

lynglib/src/commonTest/kotlin/net/sergeych/lyng/BigDecimalModuleTest.kt

@@ -135,4 +135,15 @@ class BigDecimalModuleTest {
             """.trimIndent()
         )
     }
+
+    @Test
+    fun testDefaultToString() = runTest {
+        eval("""
+            import lyng.decimal
+            
+            var s0 = "0.1".d + "0.1".d
+            assertEquals("0.2", s0.toStringExpanded())
+            assertEquals("0.2", s0.toString())
+        """.trimIndent())
+    }
 }

lynglib/src/commonTest/kotlin/net/sergeych/lyng/MatrixModuleTest.kt

@@ -0,0 +1,155 @@
+/*
+ * Copyright 2026 Sergey S. Chernov real.sergeych@gmail.com
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+
+package net.sergeych.lyng
+
+import kotlinx.coroutines.test.runTest
+import kotlin.test.Test
+
+class MatrixModuleTest {
+    @Test
+    fun testMatrixConstructionAndShape() = runTest {
+        val scope = Script.newScope()
+        scope.eval(
+            """
+            import lyng.matrix
+
+            val a: Matrix = matrix([[1, 2, 3], [4, 5, 6]])
+            val v: Vector = vector([1, 2, 3])
+            assertEquals(2, a.rows)
+            assertEquals(3, a.cols)
+            assertEquals([2, 3], a.shape)
+            assertEquals(false, a.isSquare)
+            assertEquals([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], a.toList())
+            assertEquals("Matrix(2x3, [[1, 2, 3], [4, 5, 6]])", a.toString())
+            assertEquals(3, v.size)
+            assertEquals([1.0, 2.0, 3.0], v.toList())
+            assertEquals("Vector([1, 2, 3])", v.toString())
+            assertEquals([0.2672612419124244, 0.5345224838248488, 0.8017837257372732], v.normalize().toList())
+            """.trimIndent()
+        )
+    }
+
+    @Test
+    fun testMatrixArithmeticAndTranspose() = runTest {
+        val scope = Script.newScope()
+        scope.eval(
+            """
+            import lyng.matrix
+
+            val a: Matrix = matrix([[1, 2, 3], [4, 5, 6]])
+            val b: Matrix = matrix([[7, 8], [9, 10], [11, 12]])
+            val product: Matrix = a * b
+            assertEquals([[58.0, 64.0], [139.0, 154.0]], product.toList())
+
+            val scaled: Matrix = product * 0.5
+            assertEquals([[29.0, 32.0], [69.5, 77.0]], scaled.toList())
+
+            val sum: Matrix = matrix([[1, 2], [3, 4]]) + Matrix.identity(2)
+            assertEquals([[2.0, 2.0], [3.0, 5.0]], sum.toList())
+
+            val transposed: Matrix = a.transpose()
+            assertEquals([[1.0, 4.0], [2.0, 5.0], [3.0, 6.0]], transposed.toList())
+            assertEquals([4.0, 5.0, 6.0], a.row(1))
+            assertEquals([2.0, 5.0], a.column(1))
+
+            val x: Vector = vector([1, 0.5, -1])
+            val y: Vector = a * x
+            assertEquals([-1.0, 0.5], y.toList())
+
+            val shifted: Vector = x + vector([2, 2, 2])
+            assertEquals([3.0, 2.5, 1.0], shifted.toList())
+            val d0: Vector = vector([1, 2, 3])
+            val d1: Vector = vector([2, 0, 0])
+            assertEquals(2.0, d0.dot(d1))
+            val cx0: Vector = vector([1, 0, 0])
+            val cx1: Vector = vector([0, 1, 0])
+            assertEquals([0.0, 0.0, 1.0], cx0.cross(cx1).toList())
+            val o0: Vector = vector([1.5, 3.0])
+            val o1: Vector = vector([2, 2.6666666666666665])
+            assertEquals([[3.0, 4.0], [6.0, 8.0]], o0.outer(o1).toList())
+            """.trimIndent()
+        )
+    }
+
+    @Test
+    fun testMatrixDeterminantAndInverse() = runTest {
+        val scope = Script.newScope()
+        scope.eval(
+            """
+            import lyng.matrix
+
+            val eps = 1e-9
+            val a: Matrix = matrix([[4, 7], [2, 6]])
+            assert(abs(a.determinant() - 10.0) < eps)
+            assertEquals(10.0, a.trace())
+            assertEquals(2, a.rank())
+
+            val inv: Matrix = a.inverse()
+            assert(abs(inv.get(0, 0) - 0.6) < eps)
+            assert(abs(inv.get(0, 1) + 0.7) < eps)
+            assert(abs(inv.get(1, 0) + 0.2) < eps)
+            assert(abs(inv.get(1, 1) - 0.4) < eps)
+
+            val identity: Matrix = a * inv
+            assert(abs(identity.get(0, 0) - 1.0) < eps)
+            assert(abs(identity.get(0, 1)) < eps)
+            assert(abs(identity.get(1, 0)) < eps)
+            assert(abs(identity.get(1, 1) - 1.0) < eps)
+
+            val rhs: Vector = vector([1, 0])
+            val solution: Vector = a.solve(rhs)
+            assert(abs(solution.get(0) - 0.6) < eps)
+            assert(abs(solution.get(1) + 0.2) < eps)
+
+            val rhsMatrix: Matrix = Matrix.identity(2)
+            val solvedMatrix: Matrix = a.solve(rhsMatrix)
+            assert(abs(solvedMatrix.get(0, 0) - 0.6) < eps)
+            assert(abs(solvedMatrix.get(1, 1) - 0.4) < eps)
+
+            val lowRank: Matrix = matrix([[1, 2], [2, 4]])
+            assertEquals(1, lowRank.rank())
+            """.trimIndent()
+        )
+    }
+
+    @Test
+    fun testMatrixBracketIndexingAndSlices() = runTest {
+        val scope = Script.newScope()
+        scope.eval(
+            """
+            import lyng.matrix
+
+            val m: Matrix = matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
+
+            assertEquals(7.0, m[1, 2])
+
+            val columnSlice: Matrix = m[0..2, 2]
+            assertEquals([[3.0], [7.0], [11.0]], columnSlice.toList())
+
+            val topLeft: Matrix = m[0..1, 0..1]
+            assertEquals([[1.0, 2.0], [5.0, 6.0]], topLeft.toList())
+
+            val tail: Matrix = m[1.., 1..]
+            assertEquals([[6.0, 7.0, 8.0], [10.0, 11.0, 12.0]], tail.toList())
+
+            val rowSlice: Matrix = m[1, 1..2]
+            assertEquals([[6.0, 7.0]], rowSlice.toList())
+            """.trimIndent()
+        )
+    }
+}

lynglib/src/matrixMultikMain/kotlin/net/sergeych/lyng/matrix/PlatformMatrixBackend.multik.kt

@@ -0,0 +1,53 @@
+/*
+ * Copyright 2026 Sergey S. Chernov real.sergeych@gmail.com
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+
+package net.sergeych.lyng.matrix
+
+import org.jetbrains.kotlinx.multik.api.mk
+import org.jetbrains.kotlinx.multik.api.ndarray
+import org.jetbrains.kotlinx.multik.api.linalg.dot
+import org.jetbrains.kotlinx.multik.api.linalg.inv
+import org.jetbrains.kotlinx.multik.ndarray.data.get
+
+internal actual object PlatformMatrixBackend {
+    actual fun multiply(left: MatrixData, right: MatrixData): MatrixData {
+        require(left.cols == right.rows) {
+            "matrix multiplication shape mismatch: ${left.rows}x${left.cols} cannot multiply ${right.rows}x${right.cols}"
+        }
+        val leftArray = mk.ndarray(left.values, left.rows, left.cols)
+        val rightArray = mk.ndarray(right.values, right.rows, right.cols)
+        val product = leftArray dot rightArray
+        return MatrixData(left.rows, right.cols, extract(product, left.rows, right.cols))
+    }
+
+    actual fun inverse(matrix: MatrixData): MatrixData {
+        require(matrix.isSquare) { "matrix inverse requires a square matrix, got ${matrix.rows}x${matrix.cols}" }
+        val source = mk.ndarray(matrix.values, matrix.rows, matrix.cols)
+        val inverse = mk.linalg.inv(source)
+        return MatrixData(matrix.rows, matrix.cols, extract(inverse, matrix.rows, matrix.cols))
+    }
+
+    private fun extract(array: org.jetbrains.kotlinx.multik.ndarray.data.D2Array<Double>, rows: Int, cols: Int): DoubleArray {
+        val out = DoubleArray(rows * cols)
+        for (row in 0 until rows) {
+            for (col in 0 until cols) {
+                out[row * cols + col] = array[row, col]
+            }
+        }
+        return out
+    }
+}

lynglib/src/matrixPureMain/kotlin/net/sergeych/lyng/matrix/PlatformMatrixBackend.pure.kt

@@ -0,0 +1,26 @@
+/*
+ * Copyright 2026 Sergey S. Chernov real.sergeych@gmail.com
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ *
+ */
+
+package net.sergeych.lyng.matrix
+
+internal actual object PlatformMatrixBackend {
+    actual fun multiply(left: MatrixData, right: MatrixData): MatrixData =
+        PureMatrixAlgorithms.multiply(left, right)
+
+    actual fun inverse(matrix: MatrixData): MatrixData =
+        PureMatrixAlgorithms.inverse(matrix)
+}