README.md

@@ -1,9 +1,12 @@
+[![JetBrains Research](https://jb.gg/badges/research.svg)](https://confluence.jetbrains.com/display/ALL/JetBrains+on+GitHub)
+[![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382)
+
+![Gradle build](https://github.com/mipt-npm/kmath/workflows/Gradle%20build/badge.svg)
+
 Bintray:        [ ![Download](https://api.bintray.com/packages/mipt-npm/scientifik/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion)
 
 Bintray-dev:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/scientifik/kmath-core/_latestVersion)
 
-[![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382)
-
 # KMath
 Could be pronounced as `key-math`.
 The Kotlin MATHematics library is intended as a Kotlin-based analog to Python's `numpy` library. In contrast to `numpy` and `scipy` it is modular and has a lightweight core.
@@ -40,6 +43,8 @@ can be used for a wide variety of purposes from high performance calculations to
 
 * **Streaming** Streaming operations on mathematical objects and objects buffers.
 
+* **Type-safe dimensions** Type-safe dimensions for matrix operations.
+
 * **Commons-math wrapper** It is planned to gradually wrap most parts of [Apache commons-math](http://commons.apache.org/proper/commons-math/)
                            library in Kotlin code and maybe rewrite some parts to better suit the Kotlin programming paradigm, however there is no fixed roadmap for that. Feel free
                            to submit a feature request if you want something to be done first.

build.gradle.kts

@@ -1,8 +1,8 @@
 plugins {
-    id("scientifik.publish") version "0.2.6" apply false
+    id("scientifik.publish") version "0.4.2" apply false
 }
 
-val kmathVersion by extra("0.1.4-dev-1")
+val kmathVersion by extra("0.1.4-dev-4")
 
 val bintrayRepo by extra("scientifik")
 val githubProject by extra("kmath")

examples/build.gradle.kts

@@ -4,8 +4,8 @@ import org.jetbrains.kotlin.gradle.tasks.KotlinCompile
 plugins {
     java
     kotlin("jvm")
-    kotlin("plugin.allopen") version "1.3.60"
+    kotlin("plugin.allopen") version "1.3.71"
-    id("kotlinx.benchmark") version "0.2.0-dev-5"
+    id("kotlinx.benchmark") version "0.2.0-dev-7"
 }
 
 configure<AllOpenExtension> {
@@ -13,10 +13,9 @@ configure<AllOpenExtension> {
 }
 
 repositories {
-    maven("https://dl.bintray.com/kotlin/kotlin-eap")
     maven("http://dl.bintray.com/kyonifer/maven")
-    maven ("https://dl.bintray.com/orangy/maven")
+    maven("https://dl.bintray.com/mipt-npm/scientifik")
-
+    maven("https://dl.bintray.com/mipt-npm/dev")
     mavenCentral()
 }
 
@@ -29,12 +28,11 @@ dependencies {
     implementation(project(":kmath-coroutines"))
     implementation(project(":kmath-commons"))
     implementation(project(":kmath-koma"))
+    implementation(project(":kmath-viktor"))
+    implementation(project(":kmath-dimensions"))
     implementation("com.kyonifer:koma-core-ejml:0.12")
-    implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:${Scientifik.ioVersion}")
+    implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:0.2.0-npm-dev-6")
-
+    implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-7")
-    implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-2")
-
-
     "benchmarksCompile"(sourceSets.main.get().compileClasspath)
 }
 

examples/src/benchmarks/kotlin/scientifik/kmath/structures/ViktorBenchmark.kt

@@ -0,0 +1,71 @@
+package scientifik.kmath.structures
+
+import org.jetbrains.bio.viktor.F64Array
+import org.openjdk.jmh.annotations.Benchmark
+import org.openjdk.jmh.annotations.Scope
+import org.openjdk.jmh.annotations.State
+import scientifik.kmath.operations.RealField
+import scientifik.kmath.viktor.ViktorNDField
+
+
+@State(Scope.Benchmark)
+class ViktorBenchmark {
+    final val dim = 1000
+    final val n = 100
+
+    // automatically build context most suited for given type.
+    final val autoField = NDField.auto(RealField, dim, dim)
+    final val realField = NDField.real(dim, dim)
+
+    final val viktorField = ViktorNDField(intArrayOf(dim, dim))
+
+    @Benchmark
+    fun `Automatic field addition`() {
+        autoField.run {
+            var res = one
+            repeat(n) {
+                res += 1.0
+            }
+        }
+    }
+
+    @Benchmark
+    fun `Viktor field addition`() {
+        viktorField.run {
+            var res = one
+            repeat(n) {
+                res += one
+            }
+        }
+    }
+
+    @Benchmark
+    fun `Raw Viktor`() {
+        val one = F64Array.full(init = 1.0, shape = *intArrayOf(dim, dim))
+        var res = one
+        repeat(n) {
+            res = res + one
+        }
+    }
+
+    @Benchmark
+    fun `Real field log`() {
+        realField.run {
+            val fortyTwo = produce { 42.0 }
+            var res = one
+
+            repeat(n) {
+                res = ln(fortyTwo)
+            }
+        }
+    }
+
+    @Benchmark
+    fun `Raw Viktor log`() {
+        val fortyTwo = F64Array.full(dim, dim, init = 42.0)
+        var res: F64Array
+        repeat(n) {
+            res = fortyTwo.log()
+        }
+    }
+}

examples/src/benchmarks/kotlin/scientifik/kmath/utils/utils.kt

@@ -0,0 +1,8 @@
+package scientifik.kmath.utils
+
+import kotlin.system.measureTimeMillis
+
+internal inline fun measureAndPrint(title: String, block: () -> Unit) {
+    val time = measureTimeMillis(block)
+    println("$title completed in $time millis")
+}

examples/src/main/kotlin/scientifik/kmath/structures/NDField.kt

@@ -4,7 +4,13 @@ import kotlinx.coroutines.GlobalScope
 import scientifik.kmath.operations.RealField
 import kotlin.system.measureTimeMillis
 
-fun main(args: Array<String>) {
+internal inline fun measureAndPrint(title: String, block: () -> Unit) {
+    val time = measureTimeMillis(block)
+    println("$title completed in $time millis")
+}
+
+
+fun main() {
     val dim = 1000
     val n = 1000
 
@@ -15,8 +21,7 @@ fun main(args: Array<String>) {
     //A generic boxing field. It should be used for objects, not primitives.
     val genericField = NDField.boxing(RealField, dim, dim)
 
-
+    measureAndPrint("Automatic field addition") {
-    val autoTime = measureTimeMillis {
         autoField.run {
             var res = one
             repeat(n) {
@@ -25,18 +30,14 @@ fun main(args: Array<String>) {
         }
     }
 
-    println("Automatic field addition completed in $autoTime millis")
+    measureAndPrint("Element addition"){
-
-    val elementTime = measureTimeMillis {
         var res = genericField.one
         repeat(n) {
             res += 1.0
         }
     }
 
-    println("Element addition completed in $elementTime millis")
+    measureAndPrint("Specialized addition") {
-
-    val specializedTime = measureTimeMillis {
         specializedField.run {
             var res: NDBuffer<Double> = one
             repeat(n) {
@@ -45,10 +46,7 @@ fun main(args: Array<String>) {
         }
     }
 
-    println("Specialized addition completed in $specializedTime millis")
+    measureAndPrint("Lazy addition") {
-
-
-    val lazyTime = measureTimeMillis {
         val res = specializedField.one.mapAsync(GlobalScope) {
             var c = 0.0
             repeat(n) {
@@ -60,9 +58,7 @@ fun main(args: Array<String>) {
         res.elements().forEach { it.second }
     }
 
-    println("Lazy addition completed in $lazyTime millis")
+    measureAndPrint("Generic addition") {
-
-    val genericTime = measureTimeMillis {
         //genericField.run(action)
         genericField.run {
             var res: NDBuffer<Double> = one
@@ -72,6 +68,4 @@ fun main(args: Array<String>) {
         }
     }
 
-    println("Generic addition completed in $genericTime millis")
-
 }

examples/src/main/kotlin/scientifik/kmath/structures/typeSafeDimensions.kt

@@ -0,0 +1,35 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.dimensions.D2
+import scientifik.kmath.dimensions.D3
+import scientifik.kmath.dimensions.DMatrixContext
+import scientifik.kmath.dimensions.Dimension
+import scientifik.kmath.operations.RealField
+
+fun DMatrixContext<Double, RealField>.simple() {
+    val m1 = produce<D2, D3> { i, j -> (i + j).toDouble() }
+    val m2 = produce<D3, D2> { i, j -> (i + j).toDouble() }
+
+    //Dimension-safe addition
+    m1.transpose() + m2
+}
+
+
+object D5 : Dimension {
+    override val dim: UInt = 5u
+}
+
+fun DMatrixContext<Double, RealField>.custom() {
+    val m1 = produce<D2, D5> { i, j -> (i + j).toDouble() }
+    val m2 = produce<D5, D2> { i, j -> (i - j).toDouble() }
+    val m3 = produce<D2, D2> { i, j -> (i - j).toDouble() }
+
+    (m1 dot m2) + m3
+}
+
+fun main() {
+    DMatrixContext.real.run {
+        simple()
+        custom()
+    }
+}

gradle/wrapper/gradle-wrapper.properties

@@ -1,5 +1,5 @@
 distributionBase=GRADLE_USER_HOME
 distributionPath=wrapper/dists
-distributionUrl=https\://services.gradle.org/distributions/gradle-6.0-bin.zip
+distributionUrl=https\://services.gradle.org/distributions/gradle-6.3-bin.zip
 zipStoreBase=GRADLE_USER_HOME
 zipStorePath=wrapper/dists

gradlew.bat

@@ -29,6 +29,9 @@ if "%DIRNAME%" == "" set DIRNAME=.
 set APP_BASE_NAME=%~n0
 set APP_HOME=%DIRNAME%
 
+@rem Resolve any "." and ".." in APP_HOME to make it shorter.
+for %%i in ("%APP_HOME%") do set APP_HOME=%%~fi
+
 @rem Add default JVM options here. You can also use JAVA_OPTS and GRADLE_OPTS to pass JVM options to this script.
 set DEFAULT_JVM_OPTS="-Xmx64m" "-Xms64m"
 

kmath-commons/build.gradle.kts

@@ -8,7 +8,6 @@ dependencies {
     api(project(":kmath-core"))
     api(project(":kmath-coroutines"))
     api(project(":kmath-prob"))
+    api(project(":kmath-functions"))
     api("org.apache.commons:commons-math3:3.6.1")
-    testImplementation("org.jetbrains.kotlin:kotlin-test")
-    testImplementation("org.jetbrains.kotlin:kotlin-test-junit")
 }

kmath-commons/src/test/kotlin/scientifik/kmath/commons/expressions/AutoDiffTest.kt

@@ -1,7 +1,7 @@
 package scientifik.kmath.commons.expressions
 
-import org.junit.Test
 import scientifik.kmath.expressions.invoke
+import kotlin.test.Test
 import kotlin.test.assertEquals
 
 inline fun <R> diff(order: Int, vararg parameters: Pair<String, Double>, block: DerivativeStructureField.() -> R) =

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt

@@ -41,16 +41,6 @@ fun <T : Any> Structure2D.Companion.square(vararg elements: T): FeaturedMatrix<T
     return BufferMatrix(size, size, buffer)
 }
 
-fun <T : Any> Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder<T> = MatrixBuilder(rows, columns)
-
-class MatrixBuilder<T : Any>(val rows: Int, val columns: Int) {
-    operator fun invoke(vararg elements: T): FeaturedMatrix<T> {
-        if (rows * columns != elements.size) error("The number of elements ${elements.size} is not equal $rows * $columns")
-        val buffer = elements.asBuffer()
-        return BufferMatrix(rows, columns, buffer)
-    }
-}
-
 val Matrix<*>.features get() = (this as? FeaturedMatrix)?.features?: emptySet()
 
 /**

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixBuilder.kt

@@ -0,0 +1,14 @@
+package scientifik.kmath.linear
+
+import scientifik.kmath.structures.Structure2D
+import scientifik.kmath.structures.asBuffer
+
+class MatrixBuilder<T : Any>(val rows: Int, val columns: Int) {
+    operator fun invoke(vararg elements: T): FeaturedMatrix<T> {
+        if (rows * columns != elements.size) error("The number of elements ${elements.size} is not equal $rows * $columns")
+        val buffer = elements.asBuffer()
+        return BufferMatrix(rows, columns, buffer)
+    }
+}
+
+fun <T : Any> Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder<T> = MatrixBuilder(rows, columns)

kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/functions.kt

@@ -1,37 +0,0 @@
-package scientifik.kmath.coroutines
-
-import scientifik.kmath.operations.RealField
-import scientifik.kmath.operations.SpaceOperations
-import kotlin.jvm.JvmName
-
-/**
- * A suspendable univariate function defined in algebraic context
- */
-interface UFunction<T, C : SpaceOperations<T>> {
-    suspend operator fun C.invoke(arg: T): T
-}
-
-suspend fun UFunction<Double, RealField>.invoke(arg: Double) = RealField.invoke(arg)
-
-/**
- * A suspendable multivariate (N->1) function defined on algebraic context
- */
-interface MFunction<T, C : SpaceOperations<T>> {
-    /**
-     * The input dimension of the function
-     */
-    val dimension: UInt
-
-    suspend operator fun C.invoke(vararg args: T): T
-}
-
-suspend fun MFunction<Double, RealField>.invoke(args: DoubleArray) = RealField.invoke(*args.toTypedArray())
-@JvmName("varargInvoke")
-suspend fun MFunction<Double, RealField>.invoke(vararg args: Double) = RealField.invoke(*args.toTypedArray())
-
-/**
- * A suspendable univariate function with parameter
- */
-interface ParametricUFunction<T, P, C : SpaceOperations<T>> {
-    suspend operator fun C.invoke(arg: T, parameter: P): T
-}

kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt

@@ -6,14 +6,14 @@ annotation class KMathContext
 /**
  * Marker interface for any algebra
  */
-interface Algebra
+interface Algebra<T>
 
-inline operator fun <T : Algebra, R> T.invoke(block: T.() -> R): R = run(block)
+inline operator fun <T : Algebra<*>, R> T.invoke(block: T.() -> R): R = run(block)
 
 /**
  * Space-like operations without neutral element
  */
-interface SpaceOperations<T> : Algebra {
+interface SpaceOperations<T> : Algebra<T> {
     /**
      * Addition operation for two context elements
      */

kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/AlgebraExtensions.kt

@@ -2,3 +2,14 @@ package scientifik.kmath.operations
 
 fun <T> Space<T>.sum(data: Iterable<T>): T = data.fold(zero) { left, right -> add(left, right) }
 fun <T> Space<T>.sum(data: Sequence<T>): T = data.fold(zero) { left, right -> add(left, right) }
+
+fun <T : Any, S : Space<T>> Iterable<T>.sumWith(space: S): T = space.sum(this)
+
+//TODO optimized power operation
+fun <T> RingOperations<T>.power(arg: T, power: Int): T {
+    var res = arg
+    repeat(power - 1) {
+        res *= arg
+    }
+    return res
+}

kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/BigInt.kt

@@ -0,0 +1,484 @@
+package scientifik.kmath.operations
+
+import scientifik.kmath.operations.BigInt.Companion.BASE
+import scientifik.kmath.operations.BigInt.Companion.BASE_SIZE
+import kotlin.math.log2
+import kotlin.math.max
+import kotlin.math.min
+import kotlin.math.sign
+
+
+typealias Magnitude = UIntArray
+typealias TBase = ULong
+
+/**
+ * Kotlin Multiplatform implementation of Big Integer numbers (KBigInteger).
+ *
+ * @author Robert Drynkin (https://github.com/robdrynkin) and Peter Klimai (https://github.com/pklimai)
+ */
+object BigIntField : Field<BigInt> {
+    override val zero: BigInt = BigInt.ZERO
+    override val one: BigInt = BigInt.ONE
+
+    override fun add(a: BigInt, b: BigInt): BigInt = a.plus(b)
+
+    override fun multiply(a: BigInt, k: Number): BigInt = a.times(k.toLong())
+
+    override fun multiply(a: BigInt, b: BigInt): BigInt = a.times(b)
+
+    operator fun String.unaryPlus(): BigInt = this.parseBigInteger() ?: error("Can't parse $this as big integer")
+
+    operator fun String.unaryMinus(): BigInt =
+        -(this.parseBigInteger() ?: error("Can't parse $this as big integer"))
+
+    override fun divide(a: BigInt, b: BigInt): BigInt = a.div(b)
+}
+
+class BigInt internal constructor(
+    private val sign: Byte,
+    private val magnitude: Magnitude
+) : Comparable<BigInt> {
+
+    override fun compareTo(other: BigInt): Int {
+        return when {
+            (this.sign == 0.toByte()) and (other.sign == 0.toByte()) -> 0
+            this.sign < other.sign -> -1
+            this.sign > other.sign -> 1
+            else -> this.sign * compareMagnitudes(this.magnitude, other.magnitude)
+        }
+    }
+
+    override fun equals(other: Any?): Boolean {
+        if (other is BigInt) {
+            return this.compareTo(other) == 0
+        } else error("Can't compare KBigInteger to a different type")
+    }
+
+    override fun hashCode(): Int {
+        return magnitude.hashCode() + this.sign
+    }
+
+    fun abs(): BigInt = if (sign == 0.toByte()) this else BigInt(1, magnitude)
+
+    operator fun unaryMinus(): BigInt {
+        return if (this.sign == 0.toByte()) this else BigInt((-this.sign).toByte(), this.magnitude)
+    }
+
+    operator fun plus(b: BigInt): BigInt {
+        return when {
+            b.sign == 0.toByte() -> this
+            this.sign == 0.toByte() -> b
+            this == -b -> ZERO
+            this.sign == b.sign -> BigInt(this.sign, addMagnitudes(this.magnitude, b.magnitude))
+            else -> {
+                val comp: Int = compareMagnitudes(this.magnitude, b.magnitude)
+
+                if (comp == 1) {
+                    BigInt(this.sign, subtractMagnitudes(this.magnitude, b.magnitude))
+                } else {
+                    BigInt((-this.sign).toByte(), subtractMagnitudes(b.magnitude, this.magnitude))
+                }
+            }
+        }
+    }
+
+    operator fun minus(b: BigInt): BigInt {
+        return this + (-b)
+    }
+
+    operator fun times(b: BigInt): BigInt {
+        return when {
+            this.sign == 0.toByte() -> ZERO
+            b.sign == 0.toByte() -> ZERO
+//          TODO: Karatsuba
+            else -> BigInt((this.sign * b.sign).toByte(), multiplyMagnitudes(this.magnitude, b.magnitude))
+        }
+    }
+
+    operator fun times(other: UInt): BigInt {
+        return when {
+            this.sign == 0.toByte() -> ZERO
+            other == 0U -> ZERO
+            else -> BigInt(this.sign, multiplyMagnitudeByUInt(this.magnitude, other))
+        }
+    }
+
+    operator fun times(other: Int): BigInt {
+        return if (other > 0)
+            this * kotlin.math.abs(other).toUInt()
+        else
+            -this * kotlin.math.abs(other).toUInt()
+    }
+
+    operator fun div(other: UInt): BigInt {
+        return BigInt(this.sign, divideMagnitudeByUInt(this.magnitude, other))
+    }
+
+    operator fun div(other: Int): BigInt {
+        return BigInt(
+            (this.sign * other.sign).toByte(),
+            divideMagnitudeByUInt(this.magnitude, kotlin.math.abs(other).toUInt())
+        )
+    }
+
+    private fun division(other: BigInt): Pair<BigInt, BigInt> {
+        // Long division algorithm:
+        //     https://en.wikipedia.org/wiki/Division_algorithm#Integer_division_(unsigned)_with_remainder
+        // TODO: Implement more effective algorithm
+        var q: BigInt = ZERO
+        var r: BigInt = ZERO
+
+        val bitSize =
+            (BASE_SIZE * (this.magnitude.size - 1) + log2(this.magnitude.lastOrNull()?.toFloat() ?: 0f + 1)).toInt()
+        for (i in bitSize downTo 0) {
+            r = r shl 1
+            r = r or ((abs(this) shr i) and ONE)
+            if (r >= abs(other)) {
+                r -= abs(other)
+                q += (ONE shl i)
+            }
+        }
+
+        return Pair(BigInt((this.sign * other.sign).toByte(), q.magnitude), r)
+    }
+
+    operator fun div(other: BigInt): BigInt {
+        return this.division(other).first
+    }
+
+    infix fun shl(i: Int): BigInt {
+        if (this == ZERO) return ZERO
+        if (i == 0) return this
+
+        val fullShifts = i / BASE_SIZE + 1
+        val relShift = i % BASE_SIZE
+        val shiftLeft = { x: UInt -> if (relShift >= 32) 0U else x shl relShift }
+        val shiftRight = { x: UInt -> if (BASE_SIZE - relShift >= 32) 0U else x shr (BASE_SIZE - relShift) }
+
+        val newMagnitude: Magnitude = Magnitude(this.magnitude.size + fullShifts)
+
+        for (j in this.magnitude.indices) {
+            newMagnitude[j + fullShifts - 1] = shiftLeft(this.magnitude[j])
+            if (j != 0) {
+                newMagnitude[j + fullShifts - 1] = newMagnitude[j + fullShifts - 1] or shiftRight(this.magnitude[j - 1])
+            }
+        }
+
+        newMagnitude[this.magnitude.size + fullShifts - 1] = shiftRight(this.magnitude.last())
+
+        return BigInt(this.sign, stripLeadingZeros(newMagnitude))
+    }
+
+    infix fun shr(i: Int): BigInt {
+        if (this == ZERO) return ZERO
+        if (i == 0) return this
+
+        val fullShifts = i / BASE_SIZE
+        val relShift = i % BASE_SIZE
+        val shiftRight = { x: UInt -> if (relShift >= 32) 0U else x shr relShift }
+        val shiftLeft = { x: UInt -> if (BASE_SIZE - relShift >= 32) 0U else x shl (BASE_SIZE - relShift) }
+        if (this.magnitude.size - fullShifts <= 0) {
+            return ZERO
+        }
+        val newMagnitude: Magnitude = Magnitude(this.magnitude.size - fullShifts)
+
+        for (j in fullShifts until this.magnitude.size) {
+            newMagnitude[j - fullShifts] = shiftRight(this.magnitude[j])
+            if (j != this.magnitude.size - 1) {
+                newMagnitude[j - fullShifts] = newMagnitude[j - fullShifts] or shiftLeft(this.magnitude[j + 1])
+            }
+        }
+
+        return BigInt(this.sign, stripLeadingZeros(newMagnitude))
+    }
+
+    infix fun or(other: BigInt): BigInt {
+        if (this == ZERO) return other;
+        if (other == ZERO) return this;
+        val resSize = max(this.magnitude.size, other.magnitude.size)
+        val newMagnitude: Magnitude = Magnitude(resSize)
+        for (i in 0 until resSize) {
+            if (i < this.magnitude.size) {
+                newMagnitude[i] = newMagnitude[i] or this.magnitude[i]
+            }
+            if (i < other.magnitude.size) {
+                newMagnitude[i] = newMagnitude[i] or other.magnitude[i]
+            }
+        }
+        return BigInt(1, stripLeadingZeros(newMagnitude))
+    }
+
+    infix fun and(other: BigInt): BigInt {
+        if ((this == ZERO) or (other == ZERO)) return ZERO;
+        val resSize = min(this.magnitude.size, other.magnitude.size)
+        val newMagnitude: Magnitude = Magnitude(resSize)
+        for (i in 0 until resSize) {
+            newMagnitude[i] = this.magnitude[i] and other.magnitude[i]
+        }
+        return BigInt(1, stripLeadingZeros(newMagnitude))
+    }
+
+    operator fun rem(other: Int): Int {
+        val res = this - (this / other) * other
+        return if (res == ZERO) 0 else res.sign * res.magnitude[0].toInt()
+    }
+
+    operator fun rem(other: BigInt): BigInt {
+        return this - (this / other) * other
+    }
+
+    fun modPow(exponent: BigInt, m: BigInt): BigInt {
+        return when {
+            exponent == ZERO -> ONE
+            exponent % 2 == 1 -> (this * modPow(exponent - ONE, m)) % m
+            else -> {
+                val sqRoot = modPow(exponent / 2, m)
+                (sqRoot * sqRoot) % m
+            }
+        }
+    }
+
+    override fun toString(): String {
+        if (this.sign == 0.toByte()) {
+            return "0x0"
+        }
+        var res: String = if (this.sign == (-1).toByte()) "-0x" else "0x"
+        var numberStarted = false
+
+        for (i in this.magnitude.size - 1 downTo 0) {
+            for (j in BASE_SIZE / 4 - 1 downTo 0) {
+                val curByte = (this.magnitude[i] shr 4 * j) and 0xfU
+                if (numberStarted or (curByte != 0U)) {
+                    numberStarted = true
+                    res += hexMapping[curByte]
+                }
+            }
+        }
+
+        return res
+    }
+
+    companion object {
+        const val BASE = 0xffffffffUL
+        const val BASE_SIZE: Int = 32
+        val ZERO: BigInt = BigInt(0, uintArrayOf())
+        val ONE: BigInt = BigInt(1, uintArrayOf(1u))
+
+        private val hexMapping: HashMap<UInt, String> = hashMapOf(
+            0U to "0", 1U to "1", 2U to "2", 3U to "3",
+            4U to "4", 5U to "5", 6U to "6", 7U to "7",
+            8U to "8", 9U to "9", 10U to "a", 11U to "b",
+            12U to "c", 13U to "d", 14U to "e", 15U to "f"
+        )
+
+        private fun compareMagnitudes(mag1: Magnitude, mag2: Magnitude): Int {
+            when {
+                mag1.size > mag2.size -> return 1
+                mag1.size < mag2.size -> return -1
+                else -> {
+                    for (i in mag1.size - 1 downTo 0) {
+                        if (mag1[i] > mag2[i]) {
+                            return 1
+                        } else if (mag1[i] < mag2[i]) {
+                            return -1
+                        }
+                    }
+                    return 0
+                }
+            }
+        }
+
+        private fun addMagnitudes(mag1: Magnitude, mag2: Magnitude): Magnitude {
+            val resultLength: Int = max(mag1.size, mag2.size) + 1
+            val result = Magnitude(resultLength)
+            var carry: TBase = 0UL
+
+            for (i in 0 until resultLength - 1) {
+                val res = when {
+                    i >= mag1.size -> mag2[i].toULong() + carry
+                    i >= mag2.size -> mag1[i].toULong() + carry
+                    else -> mag1[i].toULong() + mag2[i].toULong() + carry
+                }
+                result[i] = (res and BASE).toUInt()
+                carry = (res shr BASE_SIZE)
+            }
+            result[resultLength - 1] = carry.toUInt()
+            return stripLeadingZeros(result)
+        }
+
+        private fun subtractMagnitudes(mag1: Magnitude, mag2: Magnitude): Magnitude {
+            val resultLength: Int = mag1.size
+            val result = Magnitude(resultLength)
+            var carry = 0L
+
+            for (i in 0 until resultLength) {
+                var res: Long =
+                    if (i < mag2.size) mag1[i].toLong() - mag2[i].toLong() - carry
+                    else mag1[i].toLong() - carry
+
+                carry = if (res < 0) 1 else 0
+                res += carry * (BASE + 1UL).toLong()
+
+                result[i] = res.toUInt()
+            }
+
+            return stripLeadingZeros(result)
+        }
+
+        private fun multiplyMagnitudeByUInt(mag: Magnitude, x: UInt): Magnitude {
+            val resultLength: Int = mag.size + 1
+            val result = Magnitude(resultLength)
+            var carry: ULong = 0UL
+
+            for (i in mag.indices) {
+                val cur: ULong = carry + mag[i].toULong() * x.toULong()
+                result[i] = (cur and BASE.toULong()).toUInt()
+                carry = cur shr BASE_SIZE
+            }
+            result[resultLength - 1] = (carry and BASE).toUInt()
+
+            return stripLeadingZeros(result)
+        }
+
+        private fun multiplyMagnitudes(mag1: Magnitude, mag2: Magnitude): Magnitude {
+            val resultLength: Int = mag1.size + mag2.size
+            val result = Magnitude(resultLength)
+
+            for (i in mag1.indices) {
+                var carry: ULong = 0UL
+                for (j in mag2.indices) {
+                    val cur: ULong = result[i + j].toULong() + mag1[i].toULong() * mag2[j].toULong() + carry
+                    result[i + j] = (cur and BASE.toULong()).toUInt()
+                    carry = cur shr BASE_SIZE
+                }
+                result[i + mag2.size] = (carry and BASE).toUInt()
+            }
+
+            return stripLeadingZeros(result)
+        }
+
+        private fun divideMagnitudeByUInt(mag: Magnitude, x: UInt): Magnitude {
+            val resultLength: Int = mag.size
+            val result = Magnitude(resultLength)
+            var carry: ULong = 0UL
+
+            for (i in mag.size - 1 downTo 0) {
+                val cur: ULong = mag[i].toULong() + (carry shl BASE_SIZE)
+                result[i] = (cur / x).toUInt()
+                carry = cur % x
+            }
+            return stripLeadingZeros(result)
+        }
+
+    }
+
+}
+
+
+private fun stripLeadingZeros(mag: Magnitude): Magnitude {
+    if (mag.isEmpty() || mag.last() != 0U) {
+        return mag
+    }
+    var resSize: Int = mag.size - 1
+    while (mag[resSize] == 0U) {
+        if (resSize == 0)
+            break
+        resSize -= 1
+    }
+    return mag.sliceArray(IntRange(0, resSize))
+}
+
+fun abs(x: BigInt): BigInt = x.abs()
+
+/**
+ * Convert this [Int] to [BigInt]
+ */
+fun Int.toBigInt() = BigInt(sign.toByte(), uintArrayOf(kotlin.math.abs(this).toUInt()))
+
+/**
+ * Convert this [Long] to [BigInt]
+ */
+fun Long.toBigInt() = BigInt(
+    sign.toByte(), stripLeadingZeros(
+        uintArrayOf(
+            (kotlin.math.abs(this).toULong() and BASE).toUInt(),
+            ((kotlin.math.abs(this).toULong() shr BASE_SIZE) and BASE).toUInt()
+        )
+    )
+)
+
+/**
+ * Convert UInt to [BigInt]
+ */
+fun UInt.toBigInt() = BigInt(1, uintArrayOf(this))
+
+/**
+ * Convert ULong to [BigInt]
+ */
+fun ULong.toBigInt() = BigInt(
+    1,
+    stripLeadingZeros(
+        uintArrayOf(
+            (this and BigInt.BASE).toUInt(),
+            ((this shr BigInt.BASE_SIZE) and BigInt.BASE).toUInt()
+        )
+    )
+)
+
+/**
+ * Create a [BigInt] with this array of magnitudes with protective copy
+ */
+fun UIntArray.toBigInt(sign: Byte): BigInt {
+    if (sign == 0.toByte() && isNotEmpty()) error("")
+    return BigInt(sign, this.copyOf())
+}
+
+val hexChToInt = hashMapOf(
+    '0' to 0, '1' to 1, '2' to 2, '3' to 3,
+    '4' to 4, '5' to 5, '6' to 6, '7' to 7,
+    '8' to 8, '9' to 9, 'A' to 10, 'B' to 11,
+    'C' to 12, 'D' to 13, 'E' to 14, 'F' to 15
+)
+
+/**
+ * Returns null if a valid number can not be read from a string
+ */
+fun String.parseBigInteger(): BigInt? {
+    val sign: Int
+    val sPositive: String
+    when {
+        this[0] == '+' -> {
+            sign = +1
+            sPositive = this.substring(1)
+        }
+        this[0] == '-' -> {
+            sign = -1
+            sPositive = this.substring(1)
+        }
+        else -> {
+            sPositive = this
+            sign = +1
+        }
+    }
+    var res = BigInt.ZERO
+    var digitValue = BigInt.ONE
+    val sPositiveUpper = sPositive.toUpperCase()
+    if (sPositiveUpper.startsWith("0X")) {  // hex representation
+        val sHex = sPositiveUpper.substring(2)
+        for (ch in sHex.reversed()) {
+            if (ch == '_') continue
+            res += digitValue * (hexChToInt[ch] ?: return null)
+            digitValue *= 16.toBigInt()
+        }
+    } else { // decimal representation
+        for (ch in sPositiveUpper.reversed()) {
+            if (ch == '_') continue
+            if (ch !in '0'..'9') {
+                return null
+            }
+            res += digitValue * (ch.toInt() - '0'.toInt())
+            digitValue *= 10.toBigInt()
+        }
+    }
+    return res * sign
+}

kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt

@@ -29,7 +29,7 @@ fun <T : MathElement<out TrigonometricOperations<T>>> ctg(arg: T): T = arg.conte
 /**
  * A context extension to include power operations like square roots, etc
  */
-interface PowerOperations<T> {
+interface PowerOperations<T> : Algebra<T> {
     fun power(arg: T, pow: Number): T
     fun sqrt(arg: T) = power(arg, 0.5)
 
@@ -42,7 +42,7 @@ fun <T : MathElement<out PowerOperations<T>>> sqr(arg: T): T = arg pow 2.0
 
 /* Exponential */
 
-interface ExponentialOperations<T> {
+interface ExponentialOperations<T>: Algebra<T> {
     fun exp(arg: T): T
     fun ln(arg: T): T
 }

kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt

@@ -73,6 +73,8 @@ fun <T> Buffer<T>.asSequence(): Sequence<T> = Sequence(::iterator)
 
 fun <T> Buffer<T>.asIterable(): Iterable<T> = asSequence().asIterable()
 
+val Buffer<*>.indices: IntRange get() = IntRange(0, size - 1)
+
 interface MutableBuffer<T> : Buffer<T> {
     operator fun set(index: Int, value: T)
 

kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt

@@ -2,7 +2,6 @@ package scientifik.kmath.structures
 
 import scientifik.kmath.operations.*
 
-
 interface ExtendedNDField<T : Any, F, N : NDStructure<T>> :
     NDField<T, F, N>,
     TrigonometricOperations<N>,

kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt

@@ -1,26 +1,13 @@
 package scientifik.kmath.linear
 
 import scientifik.kmath.structures.Matrix
+import scientifik.kmath.structures.NDStructure
+import scientifik.kmath.structures.as2D
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
 class MatrixTest {
 
-    @Test
-    fun testSum() {
-        val vector1 = RealVector(5) { it.toDouble() }
-        val vector2 = RealVector(5) { 5 - it.toDouble() }
-        val sum = vector1 + vector2
-        assertEquals(5.0, sum[2])
-    }
-
-    @Test
-    fun testVectorToMatrix() {
-        val vector = RealVector(5) { it.toDouble() }
-        val matrix = vector.asMatrix()
-        assertEquals(4.0, matrix[4, 0])
-    }
-
     @Test
     fun testTranspose() {
         val matrix = MatrixContext.real.one(3, 3)
@@ -28,21 +15,6 @@ class MatrixTest {
         assertEquals(matrix, transposed)
     }
 
-
-    @Test
-    fun testDot() {
-        val vector1 = RealVector(5) { it.toDouble() }
-        val vector2 = RealVector(5) { 5 - it.toDouble() }
-
-        val matrix1 = vector1.asMatrix()
-        val matrix2 = vector2.asMatrix().transpose()
-        val product = MatrixContext.real.run { matrix1 dot matrix2 }
-
-
-        assertEquals(5.0, product[1, 0])
-        assertEquals(6.0, product[2, 2])
-    }
-
     @Test
     fun testBuilder() {
         val matrix = Matrix.build<Double>(2, 3)(
@@ -74,4 +46,20 @@ class MatrixTest {
 
         val toTenthPower = transitionMatrix pow 10
     }
+
+    @Test
+    fun test2DDot() {
+        val firstMatrix = NDStructure.auto(2,3){ (i, j) -> (i + j).toDouble() }.as2D()
+        val secondMatrix = NDStructure.auto(3,2){ (i, j) -> (i + j).toDouble() }.as2D()
+        MatrixContext.real.run {
+//            val firstMatrix = produce(2, 3) { i, j -> (i + j).toDouble() }
+//            val secondMatrix = produce(3, 2) { i, j -> (i + j).toDouble() }
+            val result = firstMatrix dot secondMatrix
+            assertEquals(2, result.rowNum)
+            assertEquals(2, result.colNum)
+            assertEquals(8.0, result[0,1])
+            assertEquals(8.0, result[1,0])
+            assertEquals(14.0, result[1,1])
+        }
+    }
 }

kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntAlgebraTest.kt

@@ -0,0 +1,50 @@
+package scientifik.kmath.operations
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class BigIntAlgebraTest {
+    @Test
+    fun testKBigIntegerRingSum() {
+        val res = BigIntField {
+            1_000L.toBigInt() * 1_000L.toBigInt()
+        }
+        assertEquals(res, 1_000_000.toBigInt())
+    }
+
+    @Test
+    fun testKBigIntegerRingSum_100_000_000__100_000_000() {
+        BigIntField {
+            val sum = +"100_000_000" + +"100_000_000"
+            assertEquals(sum, "200_000_000".parseBigInteger())
+        }
+    }
+
+    @Test
+    fun test_mul_3__4() {
+        BigIntField {
+            val prod = +"0x3000_0000_0000" * +"0x4000_0000_0000_0000_0000"
+            assertEquals(prod, "0xc00_0000_0000_0000_0000_0000_0000_0000".parseBigInteger())
+        }
+    }
+
+    @Test
+    fun test_div_big_1() {
+        BigIntField {
+            val res = +"1_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000" /
+                    +"555_000_444_000_333_000_222_000_111_000_999_001"
+            assertEquals(res, +"1801800360360432432518919022699")
+        }
+    }
+
+    @Test
+    fun test_rem_big_1() {
+        BigIntField {
+            val res = +"1_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000_000" %
+                    +"555_000_444_000_333_000_222_000_111_000_999_001"
+            assertEquals(res, +"324121220440768000291647788404676301")
+        }
+    }
+
+}
+

kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntConstructorTest.kt

@@ -0,0 +1,26 @@
+package scientifik.kmath.operations
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class BigIntConstructorTest {
+    @Test
+    fun testConstructorZero() {
+        assertEquals(0.toBigInt(), uintArrayOf().toBigInt(0))
+    }
+
+    @Test
+    fun testConstructor8() {
+        assertEquals(
+            8.toBigInt(),
+            uintArrayOf(8U).toBigInt(1)
+        )
+    }
+
+    @Test
+    fun testConstructor_0xffffffffaL() {
+        val x = -0xffffffffaL.toBigInt()
+        val y = uintArrayOf(0xfffffffaU, 0xfU).toBigInt(-1)
+        assertEquals(x, y)
+    }
+}

kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntConversionsTest.kt

@@ -0,0 +1,43 @@
+package scientifik.kmath.operations
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+@kotlin.ExperimentalUnsignedTypes
+class BigIntConversionsTest {
+    @Test
+    fun testToString0x10() {
+        val x = 0x10.toBigInt()
+        assertEquals("0x10", x.toString())
+    }
+
+    @Test
+    fun testToString0x17ffffffd() {
+        val x = 0x17ffffffdL.toBigInt()
+        assertEquals("0x17ffffffd", x.toString())
+    }
+
+    @Test
+    fun testToString_0x17ead2ffffd() {
+        val x = -0x17ead2ffffdL.toBigInt()
+        assertEquals("-0x17ead2ffffd", x.toString())
+    }
+
+    @Test
+    fun testToString_0x17ead2ffffd11223344() {
+        val x = uintArrayOf(0x11223344U, 0xad2ffffdU, 0x17eU).toBigInt(-1)
+        assertEquals("-0x17ead2ffffd11223344", x.toString())
+    }
+
+    @Test
+    fun testFromString_0x17ead2ffffd11223344() {
+        val x = "0x17ead2ffffd11223344".parseBigInteger()
+        assertEquals("0x17ead2ffffd11223344", x.toString())
+    }
+
+    @Test
+    fun testFromString_7059135710711894913860() {
+        val x = "-7059135710711894913860".parseBigInteger()
+        assertEquals("-0x17ead2ffffd11223344", x.toString())
+    }
+}

kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/BigIntOperationsTest.kt

@@ -0,0 +1,381 @@
+package scientifik.kmath.operations
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+@kotlin.ExperimentalUnsignedTypes
+class BigIntOperationsTest {
+    @Test
+    fun testPlus_1_1() {
+        val x = 1.toBigInt()
+        val y = 1.toBigInt()
+
+        val res = x + y
+        val sum = 2.toBigInt()
+
+        assertEquals(sum, res)
+    }
+
+    @Test
+    fun testPlusBigNumbers() {
+        val x = 0x7fffffff.toBigInt()
+        val y = 0x7fffffff.toBigInt()
+        val z = 0x7fffffff.toBigInt()
+
+        val res = x + y + z
+        val sum = uintArrayOf(0x7ffffffdU, 0x1U).toBigInt(1)
+
+        assertEquals(sum, res)
+    }
+
+    @Test
+    fun testUnaryMinus() {
+        val x = 1234.toBigInt()
+        val y = -1234.toBigInt()
+        assertEquals(-x, y)
+    }
+
+    @Test
+    fun testMinus_2_1() {
+        val x = 2.toBigInt()
+        val y = 1.toBigInt()
+
+        val res = x - y
+        val sum = 1.toBigInt()
+
+        assertEquals(sum, res)
+    }
+
+    @Test
+    fun testMinus__2_1() {
+        val x = -2.toBigInt()
+        val y = 1.toBigInt()
+
+        val res = x - y
+        val sum = -3.toBigInt()
+
+        assertEquals(sum, res)
+    }
+
+    @Test
+    fun testMinus___2_1() {
+        val x = -2.toBigInt()
+        val y = 1.toBigInt()
+
+        val res = -x - y
+        val sum = 1.toBigInt()
+
+        assertEquals(sum, res)
+    }
+
+    @Test
+    fun testMinusBigNumbers() {
+        val x = 12345.toBigInt()
+        val y = 0xffffffffaL.toBigInt()
+
+        val res = x - y
+        val sum = -0xfffffcfc1L.toBigInt()
+
+        assertEquals(sum, res)
+    }
+
+    @Test
+    fun testMultiply_2_3() {
+        val x = 2.toBigInt()
+        val y = 3.toBigInt()
+
+        val res = x * y
+        val prod = 6.toBigInt()
+
+        assertEquals(prod, res)
+    }
+
+    @Test
+    fun testMultiply__2_3() {
+        val x = -2.toBigInt()
+        val y = 3.toBigInt()
+
+        val res = x * y
+        val prod = -6.toBigInt()
+
+        assertEquals(prod, res)
+    }
+
+    @Test
+    fun testMultiply_0xfff123_0xfff456() {
+        val x = 0xfff123.toBigInt()
+        val y = 0xfff456.toBigInt()
+
+        val res = x * y
+        val prod = 0xffe579ad5dc2L.toBigInt()
+
+        assertEquals(prod, res)
+    }
+
+    @Test
+    fun testMultiplyUInt_0xfff123_0xfff456() {
+        val x = 0xfff123.toBigInt()
+        val y = 0xfff456U
+
+        val res = x * y
+        val prod = 0xffe579ad5dc2L.toBigInt()
+
+        assertEquals(prod, res)
+    }
+
+    @Test
+    fun testMultiplyInt_0xfff123__0xfff456() {
+        val x = 0xfff123.toBigInt()
+        val y = -0xfff456
+
+        val res = x * y
+        val prod = -0xffe579ad5dc2L.toBigInt()
+
+        assertEquals(prod, res)
+    }
+
+    @Test
+    fun testMultiply_0xffffffff_0xffffffff() {
+        val x = 0xffffffffL.toBigInt()
+        val y = 0xffffffffL.toBigInt()
+
+        val res = x * y
+        val prod = 0xfffffffe00000001UL.toBigInt()
+
+        assertEquals(prod, res)
+    }
+
+    @Test
+    fun test_shr_20() {
+        val x = 20.toBigInt()
+        assertEquals(10.toBigInt(), x shr 1)
+    }
+
+    @Test
+    fun test_shl_20() {
+        val x = 20.toBigInt()
+        assertEquals(40.toBigInt(), x shl 1)
+    }
+
+    @Test
+    fun test_shl_1_0() {
+        assertEquals(
+            BigInt.ONE,
+            BigInt.ONE shl 0
+        )
+    }
+
+    @Test
+    fun test_shl_1_32() {
+        assertEquals(
+            0x100000000UL.toBigInt(),
+            BigInt.ONE shl 32
+        )
+    }
+
+    @Test
+    fun test_shl_1_33() {
+        assertEquals(
+            0x200000000UL.toBigInt(),
+            BigInt.ONE shl 33
+        )
+    }
+
+    @Test
+    fun test_shr_1_33_33() {
+        assertEquals(
+            BigInt.ONE,
+            (BigInt.ONE shl 33) shr 33
+        )
+    }
+
+    @Test
+    fun test_shr_1_32() {
+        assertEquals(
+            BigInt.ZERO,
+            BigInt.ONE shr 32
+        )
+    }
+
+    @Test
+    fun test_and_123_456() {
+        val x = 123.toBigInt()
+        val y = 456.toBigInt()
+        assertEquals(72.toBigInt(), x and y)
+    }
+
+    @Test
+    fun test_or_123_456() {
+        val x = 123.toBigInt()
+        val y = 456.toBigInt()
+        assertEquals(507.toBigInt(), x or y)
+    }
+
+    @Test
+    fun test_asd() {
+        assertEquals(
+            BigInt.ONE,
+            BigInt.ZERO or ((20.toBigInt() shr 4) and BigInt.ONE)
+        )
+    }
+
+    @Test
+    fun test_square_0x11223344U_0xad2ffffdU_0x17eU() {
+        val num =
+            uintArrayOf(0x11223344U, 0xad2ffffdU, 0x17eU).toBigInt(-1)
+        println(num)
+        val res = num * num
+        assertEquals(
+            res,
+            uintArrayOf(0xb0542a10U, 0xbbd85bc8U, 0x2a1fa515U, 0x5069e03bU, 0x23c09U).toBigInt(1)
+        )
+    }
+
+    @Test
+    fun testDivision_6_3() {
+        val x = 6.toBigInt()
+        val y = 3U
+
+        val res = x / y
+        val div = 2.toBigInt()
+
+        assertEquals(div, res)
+    }
+
+    @Test
+    fun testBigDivision_6_3() {
+        val x = 6.toBigInt()
+        val y = 3.toBigInt()
+
+        val res = x / y
+        val div = 2.toBigInt()
+
+        assertEquals(div, res)
+    }
+
+    @Test
+    fun testDivision_20__3() {
+        val x = 20.toBigInt()
+        val y = -3
+
+        val res = x / y
+        val div = -6.toBigInt()
+
+        assertEquals(div, res)
+    }
+
+    @Test
+    fun testBigDivision_20__3() {
+        val x = 20.toBigInt()
+        val y = -3.toBigInt()
+
+        val res = x / y
+        val div = -6.toBigInt()
+
+        assertEquals(div, res)
+    }
+
+    @Test
+    fun testDivision_0xfffffffe00000001_0xffffffff() {
+        val x = 0xfffffffe00000001UL.toBigInt()
+        val y = 0xffffffffU
+
+        val res = x / y
+        val div = 0xffffffffL.toBigInt()
+
+        assertEquals(div, res)
+    }
+
+    @Test
+    fun testBigDivision_0xfffffffeabcdef01UL_0xfffffffeabc() {
+        val res = 0xfffffffeabcdef01UL.toBigInt() / 0xfffffffeabc.toBigInt()
+        assertEquals(res, 0x100000.toBigInt())
+    }
+
+    @Test
+    fun testBigDivision_0xfffffffe00000001_0xffffffff() {
+        val x = 0xfffffffe00000001UL.toBigInt()
+        val y = 0xffffffffU.toBigInt()
+
+        val res = x / y
+        val div = 0xffffffffL.toBigInt()
+
+        assertEquals(div, res)
+    }
+
+    @Test
+    fun testMod_20_3() {
+        val x = 20.toBigInt()
+        val y = 3
+
+        val res = x % y
+        val mod = 2
+
+        assertEquals(mod, res)
+    }
+
+    @Test
+    fun testBigMod_20_3() {
+        val x = 20.toBigInt()
+        val y = 3.toBigInt()
+
+        val res = x % y
+        val mod = 2.toBigInt()
+
+        assertEquals(mod, res)
+    }
+
+    @Test
+    fun testMod_0xfffffffe00000001_12345() {
+        val x = 0xfffffffe00000001UL.toBigInt()
+        val y = 12345
+
+        val res = x % y
+        val mod = 1980
+
+        assertEquals(mod, res)
+    }
+
+    @Test
+    fun testBigMod_0xfffffffe00000001_12345() {
+        val x = 0xfffffffe00000001UL.toBigInt()
+        val y = 12345.toBigInt()
+
+        val res = x % y
+        val mod = 1980.toBigInt()
+
+        assertEquals(mod, res)
+    }
+
+    @Test
+    fun testModPow_3_10_17() {
+        val x = 3.toBigInt()
+        val exp = 10.toBigInt()
+        val mod = 17.toBigInt()
+
+        val res = 8.toBigInt()
+
+        return assertEquals(res, x.modPow(exp, mod))
+    }
+
+    @Test
+    fun testModPowBigNumbers() {
+        val x = 0xfffffffeabcdef01UL.toBigInt()
+        val exp = 2.toBigInt()
+        val mod = 0xfffffffeabcUL.toBigInt()
+
+        val res = 0xc2253cde01.toBigInt()
+
+        return assertEquals(res, x.modPow(exp, mod))
+    }
+
+    @Test
+    fun testModBigNumbers() {
+        val x = 0xfffffffeabcdef01UL.toBigInt()
+        val mod = 0xfffffffeabcUL.toBigInt()
+
+        val res = 0xdef01.toBigInt()
+
+        return assertEquals(res, x % mod)
+    }
+}

kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt

@@ -16,7 +16,9 @@
 
 package scientifik.kmath.chains
 
+import kotlinx.coroutines.InternalCoroutinesApi
 import kotlinx.coroutines.flow.Flow
+import kotlinx.coroutines.flow.FlowCollector
 import kotlinx.coroutines.sync.Mutex
 import kotlinx.coroutines.sync.withLock
 
@@ -25,7 +27,7 @@ import kotlinx.coroutines.sync.withLock
  * A not-necessary-Markov chain of some type
  * @param R - the chain element type
  */
-interface Chain<out R> {
+interface Chain<out R>: Flow<R> {
     /**
      * Generate next value, changing state if needed
      */
@@ -36,14 +38,15 @@ interface Chain<out R> {
      */
     fun fork(): Chain<R>
 
+    @InternalCoroutinesApi
+    override suspend fun collect(collector: FlowCollector<R>) {
+        kotlinx.coroutines.flow.flow { while (true) emit(next()) }.collect(collector)
+    }
+
     companion object
 
 }
 
-/**
- * Chain as a coroutine flow. The flow emit affects chain state and vice versa
- */
-fun <R> Chain<R>.flow(): Flow<R> = kotlinx.coroutines.flow.flow { while (true) emit(next()) }
 
 fun <T> Iterator<T>.asChain(): Chain<T> = SimpleChain { next() }
 fun <T> Sequence<T>.asChain(): Chain<T> = iterator().asChain()
@@ -127,6 +130,21 @@ fun <T, R> Chain<T>.map(func: suspend (T) -> R): Chain<R> = object : Chain<R> {
     override fun fork(): Chain<R> = this@map.fork().map(func)
 }
 
+/**
+ * [block] must be a pure function or at least not use external random variables, otherwise fork could be broken
+ */
+fun <T> Chain<T>.filter(block: (T) -> Boolean): Chain<T> = object : Chain<T> {
+    override suspend fun next(): T {
+        var next: T
+        do {
+            next = this@filter.next()
+        } while (!block(next))
+        return next
+    }
+
+    override fun fork(): Chain<T> = this@filter.fork().filter(block)
+}
+
 /**
  * Map the whole chain
  */

kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/BufferFlowTest.kt

@@ -3,11 +3,12 @@ package scientifik.kmath.streaming
 import kotlinx.coroutines.*
 import kotlinx.coroutines.flow.asFlow
 import kotlinx.coroutines.flow.collect
-import org.junit.Test
+import org.junit.jupiter.api.Timeout
 import scientifik.kmath.coroutines.async
 import scientifik.kmath.coroutines.collect
 import scientifik.kmath.coroutines.mapParallel
 import java.util.concurrent.Executors
+import kotlin.test.Test
 
 
 @ExperimentalCoroutinesApi
@@ -17,7 +18,8 @@ class BufferFlowTest {
 
     val dispatcher = Executors.newFixedThreadPool(4).asCoroutineDispatcher()
 
-    @Test(timeout = 2000)
+    @Test
+    @Timeout(2000)
     fun map() {
         runBlocking {
             (1..20).asFlow().mapParallel( dispatcher) {
@@ -31,7 +33,8 @@ class BufferFlowTest {
         }
     }
 
-    @Test(timeout = 2000)
+    @Test
+    @Timeout(2000)
     fun async() {
         runBlocking {
             (1..20).asFlow().async(dispatcher) {

kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/RingBufferTest.kt

@@ -2,8 +2,8 @@ package scientifik.kmath.streaming
 
 import kotlinx.coroutines.flow.*
 import kotlinx.coroutines.runBlocking
-import org.junit.Test
 import scientifik.kmath.structures.asSequence
+import kotlin.test.Test
 import kotlin.test.assertEquals
 
 class RingBufferTest {

kmath-dimensions/build.gradle.kts

@@ -0,0 +1,19 @@
+plugins {
+    id("scientifik.mpp")
+}
+
+description = "A proof of concept module for adding typ-safe dimensions to structures"
+
+kotlin.sourceSets {
+    commonMain {
+        dependencies {
+            api(project(":kmath-core"))
+        }
+    }
+
+    jvmMain{
+        dependencies{
+            api(kotlin("reflect"))
+        }
+    }
+}

kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Dimensions.kt

@@ -0,0 +1,35 @@
+package scientifik.kmath.dimensions
+
+import kotlin.reflect.KClass
+
+/**
+ * An abstract class which is not used in runtime. Designates a size of some structure.
+ * Could be replaced later by fully inline constructs
+ */
+interface Dimension {
+
+    val dim: UInt
+    companion object {
+
+    }
+}
+
+fun <D : Dimension> KClass<D>.dim(): UInt = Dimension.resolve(this).dim
+
+expect fun <D : Dimension> Dimension.Companion.resolve(type: KClass<D>): D
+
+expect fun Dimension.Companion.of(dim: UInt): Dimension
+
+inline fun <reified D : Dimension> Dimension.Companion.dim(): UInt = D::class.dim()
+
+object D1 : Dimension {
+    override val dim: UInt get() = 1U
+}
+
+object D2 : Dimension {
+    override val dim: UInt get() = 2U
+}
+
+object D3 : Dimension {
+    override val dim: UInt get() = 31U
+}

kmath-dimensions/src/commonMain/kotlin/scientifik/kmath/dimensions/Wrappers.kt

@@ -0,0 +1,161 @@
+package scientifik.kmath.dimensions
+
+import scientifik.kmath.linear.GenericMatrixContext
+import scientifik.kmath.linear.MatrixContext
+import scientifik.kmath.linear.Point
+import scientifik.kmath.linear.transpose
+import scientifik.kmath.operations.Ring
+import scientifik.kmath.structures.Matrix
+import scientifik.kmath.structures.Structure2D
+
+/**
+ * A matrix with compile-time controlled dimension
+ */
+interface DMatrix<T, R : Dimension, C : Dimension> : Structure2D<T> {
+    companion object {
+        /**
+         * Coerces a regular matrix to a matrix with type-safe dimensions and throws a error if coercion failed
+         */
+        inline fun <T, reified R : Dimension, reified C : Dimension> coerce(structure: Structure2D<T>): DMatrix<T, R, C> {
+            if (structure.rowNum != Dimension.dim<R>().toInt()) {
+                error("Row number mismatch: expected ${Dimension.dim<R>()} but found ${structure.rowNum}")
+            }
+            if (structure.colNum != Dimension.dim<C>().toInt()) {
+                error("Column number mismatch: expected ${Dimension.dim<C>()} but found ${structure.colNum}")
+            }
+            return DMatrixWrapper(structure)
+        }
+
+        /**
+         * The same as [coerce] but without dimension checks. Use with caution
+         */
+        fun <T, R : Dimension, C : Dimension> coerceUnsafe(structure: Structure2D<T>): DMatrix<T, R, C> {
+            return DMatrixWrapper(structure)
+        }
+    }
+}
+
+/**
+ * An inline wrapper for a Matrix
+ */
+inline class DMatrixWrapper<T, R : Dimension, C : Dimension>(
+    val structure: Structure2D<T>
+) : DMatrix<T, R, C> {
+    override val shape: IntArray get() = structure.shape
+    override fun get(i: Int, j: Int): T = structure[i, j]
+}
+
+/**
+ * Dimension-safe point
+ */
+interface DPoint<T, D : Dimension> : Point<T> {
+    companion object {
+        inline fun <T, reified D : Dimension> coerce(point: Point<T>): DPoint<T, D> {
+            if (point.size != Dimension.dim<D>().toInt()) {
+                error("Vector dimension mismatch: expected ${Dimension.dim<D>()}, but found ${point.size}")
+            }
+            return DPointWrapper(point)
+        }
+
+        fun <T, D : Dimension> coerceUnsafe(point: Point<T>): DPoint<T, D> {
+            return DPointWrapper(point)
+        }
+    }
+}
+
+/**
+ * Dimension-safe point wrapper
+ */
+inline class DPointWrapper<T, D : Dimension>(val point: Point<T>) :
+    DPoint<T, D> {
+    override val size: Int get() = point.size
+
+    override fun get(index: Int): T = point[index]
+
+    override fun iterator(): Iterator<T> = point.iterator()
+}
+
+
+/**
+ * Basic operations on dimension-safe matrices. Operates on [Matrix]
+ */
+inline class DMatrixContext<T : Any, Ri : Ring<T>>(val context: GenericMatrixContext<T, Ri>) {
+
+    inline fun <reified R : Dimension, reified C : Dimension> Matrix<T>.coerce(): DMatrix<T, R, C> {
+        if (rowNum != Dimension.dim<R>().toInt()) {
+            error("Row number mismatch: expected ${Dimension.dim<R>()} but found $rowNum")
+        }
+        if (colNum != Dimension.dim<C>().toInt()) {
+            error("Column number mismatch: expected ${Dimension.dim<C>()} but found $colNum")
+        }
+        return DMatrix.coerceUnsafe(this)
+    }
+
+    /**
+     * Produce a matrix with this context and given dimensions
+     */
+    inline fun <reified R : Dimension, reified C : Dimension> produce(noinline initializer: (i: Int, j: Int) -> T): DMatrix<T, R, C> {
+        val rows = Dimension.dim<R>()
+        val cols = Dimension.dim<C>()
+        return context.produce(rows.toInt(), cols.toInt(), initializer).coerce<R,C>()
+    }
+
+    inline fun <reified D : Dimension> point(noinline initializer: (Int) -> T): DPoint<T, D> {
+        val size = Dimension.dim<D>()
+        return DPoint.coerceUnsafe(
+            context.point(
+                size.toInt(),
+                initializer
+            )
+        )
+    }
+
+    inline infix fun <reified R1 : Dimension, reified C1 : Dimension, reified C2 : Dimension> DMatrix<T, R1, C1>.dot(
+        other: DMatrix<T, C1, C2>
+    ): DMatrix<T, R1, C2> {
+        return context.run { this@dot dot other }.coerce()
+    }
+
+    inline infix fun <reified R : Dimension, reified C : Dimension> DMatrix<T, R, C>.dot(vector: DPoint<T, C>): DPoint<T, R> {
+        return DPoint.coerceUnsafe(context.run { this@dot dot vector })
+    }
+
+    inline operator fun <reified R : Dimension, reified C : Dimension> DMatrix<T, R, C>.times(value: T): DMatrix<T, R, C> {
+        return context.run { this@times.times(value) }.coerce()
+    }
+
+    inline operator fun <reified R : Dimension, reified C : Dimension> T.times(m: DMatrix<T, R, C>): DMatrix<T, R, C> =
+        m * this
+
+
+    inline operator fun <reified R : Dimension, reified C : Dimension> DMatrix<T, C, R>.plus(other: DMatrix<T, C, R>): DMatrix<T, C, R> {
+        return context.run { this@plus + other }.coerce()
+    }
+
+    inline operator fun <reified R : Dimension, reified C : Dimension> DMatrix<T, C, R>.minus(other: DMatrix<T, C, R>): DMatrix<T, C, R> {
+        return context.run { this@minus + other }.coerce()
+    }
+
+    inline operator fun <reified R : Dimension, reified C : Dimension> DMatrix<T, C, R>.unaryMinus(): DMatrix<T, C, R> {
+        return context.run { this@unaryMinus.unaryMinus() }.coerce()
+    }
+
+    inline fun <reified R : Dimension, reified C : Dimension> DMatrix<T, C, R>.transpose(): DMatrix<T, R, C> {
+        return context.run { (this@transpose as Matrix<T>).transpose() }.coerce()
+    }
+
+    /**
+     * A square unit matrix
+     */
+    inline fun <reified D : Dimension> one(): DMatrix<T, D, D> = produce { i, j ->
+        if (i == j) context.elementContext.one else context.elementContext.zero
+    }
+
+    inline fun <reified R : Dimension, reified C : Dimension> zero(): DMatrix<T, R, C> = produce { _, _ ->
+        context.elementContext.zero
+    }
+
+    companion object {
+        val real = DMatrixContext(MatrixContext.real)
+    }
+}

kmath-dimensions/src/commonTest/kotlin/scientifik/dimensions/DMatrixContextTest.kt

@@ -0,0 +1,30 @@
+package scientifik.dimensions
+
+import scientifik.kmath.dimensions.D2
+import scientifik.kmath.dimensions.D3
+import scientifik.kmath.dimensions.DMatrixContext
+import kotlin.test.Test
+
+
+class DMatrixContextTest {
+    @Test
+    fun testDimensionSafeMatrix() {
+        val res = DMatrixContext.real.run {
+            val m = produce<D2, D2> { i, j -> (i + j).toDouble() }
+
+            //The dimension of `one()` is inferred from type
+            (m + one())
+        }
+    }
+
+    @Test
+    fun testTypeCheck() {
+        val res = DMatrixContext.real.run {
+            val m1 = produce<D2, D3> { i, j -> (i + j).toDouble() }
+            val m2 = produce<D3, D2> { i, j -> (i + j).toDouble() }
+
+            //Dimension-safe addition
+            m1.transpose() + m2
+        }
+    }
+}

kmath-dimensions/src/jsMain/kotlin/scientifik/kmath/dimensions/dim.kt

@@ -0,0 +1,22 @@
+package scientifik.kmath.dimensions
+
+import kotlin.reflect.KClass
+
+private val dimensionMap = hashMapOf<UInt, Dimension>(
+    1u to D1,
+    2u to D2,
+    3u to D3
+)
+
+@Suppress("UNCHECKED_CAST")
+actual fun <D : Dimension> Dimension.Companion.resolve(type: KClass<D>): D {
+    return dimensionMap.entries.find { it.value::class == type }?.value as? D ?: error("Can't resolve dimension $type")
+}
+
+actual fun Dimension.Companion.of(dim: UInt): Dimension {
+    return dimensionMap.getOrPut(dim) {
+        object : Dimension {
+            override val dim: UInt get() = dim
+        }
+    }
+}

kmath-dimensions/src/jvmMain/kotlin/scientifik/kmath/dimensions/dim.kt

@@ -0,0 +1,18 @@
+package scientifik.kmath.dimensions
+
+import kotlin.reflect.KClass
+
+actual fun <D:Dimension> Dimension.Companion.resolve(type: KClass<D>): D{
+    return type.objectInstance ?: error("No object instance for dimension class")
+}
+
+actual fun Dimension.Companion.of(dim: UInt): Dimension{
+    return when(dim){
+        1u -> D1
+        2u -> D2
+        3u -> D3
+        else -> object : Dimension {
+            override val dim: UInt get() = dim
+        }
+    }
+}

kmath-for-real/build.gradle.kts

@@ -0,0 +1,11 @@
+plugins {
+    id("scientifik.mpp")
+}
+
+kotlin.sourceSets {
+    commonMain {
+        dependencies {
+            api(project(":kmath-core"))
+        }
+    }
+}

kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt

@@ -0,0 +1,146 @@
+package scientifik.kmath.real
+
+import scientifik.kmath.linear.MatrixContext
+import scientifik.kmath.linear.RealMatrixContext.elementContext
+import scientifik.kmath.linear.VirtualMatrix
+import scientifik.kmath.operations.sum
+import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.Matrix
+import scientifik.kmath.structures.asSequence
+import kotlin.math.pow
+
+/*
+ *  Functions for convenient "numpy-like" operations with Double matrices.
+ *
+ *  Initial implementation of these functions is taken from:
+ *    https://github.com/thomasnield/numky/blob/master/src/main/kotlin/org/nield/numky/linear/DoubleOperators.kt
+ *
+ */
+
+/*
+ *  Functions that help create a real (Double) matrix
+ */
+
+fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double) =
+        MatrixContext.real.produce(rowNum, colNum, initializer)
+
+fun Sequence<DoubleArray>.toMatrix() = toList().let {
+    MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] }
+}
+
+fun Matrix<Double>.repeatStackVertical(n: Int) = VirtualMatrix(rowNum*n, colNum) {
+    row, col -> get(if (row == 0) 0 else row % rowNum, col)
+}
+
+/*
+ *  Operations for matrix and real number
+ */
+
+operator fun Matrix<Double>.times(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
+    row, col -> this[row, col] * double
+}
+
+operator fun Matrix<Double>.plus(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
+    row, col -> this[row, col] + double
+}
+
+operator fun Matrix<Double>.minus(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
+    row, col -> this[row, col] - double
+}
+
+operator fun Matrix<Double>.div(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
+    row, col -> this[row, col] / double
+}
+
+operator fun Double.times(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
+    row, col -> this * matrix[row, col]
+}
+
+operator fun Double.plus(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
+    row, col -> this * matrix[row, col]
+}
+
+operator fun Double.minus(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
+    row, col -> this - matrix[row, col]
+}
+
+// TODO: does this operation make sense? Should it be 'this/matrix[row, col]'?
+//operator fun Double.div(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
+//    row, col -> matrix[row, col] / this
+//}
+
+/*
+ *  Per-element (!) square and power operations
+ */
+
+fun Matrix<Double>.square() = MatrixContext.real.produce(rowNum, colNum) {
+    row, col -> this[row, col].pow(2)
+}
+
+fun Matrix<Double>.pow(n: Int) = MatrixContext.real.produce(rowNum, colNum) {
+    i, j -> this[i, j].pow(n)
+}
+
+/*
+ * Operations on two matrices (per-element!)
+ */
+
+operator fun Matrix<Double>.times(other: Matrix<Double>) = MatrixContext.real.produce(rowNum, colNum) {
+    row, col -> this[row, col] * other[row, col]
+}
+
+operator fun Matrix<Double>.plus(other: Matrix<Double>) = MatrixContext.real.add(this, other)
+
+operator fun Matrix<Double>.minus(other: Matrix<Double>) = MatrixContext.real.produce(rowNum, colNum) {
+    row, col -> this[row,col] - other[row,col]
+}
+
+/*
+ *  Operations on columns
+ */
+
+inline fun Matrix<Double>.appendColumn(crossinline mapper: (Buffer<Double>) -> Double) =
+    MatrixContext.real.produce(rowNum,colNum+1) {
+        row, col ->
+            if (col < colNum)
+                this[row, col]
+            else
+                mapper(rows[row])
+    }
+
+fun Matrix<Double>.extractColumns(columnRange: IntRange) = MatrixContext.real.produce(rowNum, columnRange.count()) {
+    row, col -> this[row, columnRange.first + col]
+}
+
+fun Matrix<Double>.extractColumn(columnIndex: Int) = extractColumns(columnIndex..columnIndex)
+
+fun Matrix<Double>.sumByColumn() = MatrixContext.real.produce(1, colNum) { _, j ->
+    val column = columns[j]
+    with(elementContext) {
+        sum(column.asSequence())
+    }
+}
+
+fun Matrix<Double>.minByColumn() = MatrixContext.real.produce(1, colNum) {
+    _, j -> columns[j].asSequence().min() ?: throw Exception("Cannot produce min on empty column")
+}
+
+fun Matrix<Double>.maxByColumn() = MatrixContext.real.produce(1, colNum) {
+    _, j -> columns[j].asSequence().max() ?: throw Exception("Cannot produce min on empty column")
+}
+
+fun Matrix<Double>.averageByColumn() = MatrixContext.real.produce(1, colNum) {
+    _, j -> columns[j].asSequence().average()
+}
+
+/*
+ * Operations processing all elements
+ */
+
+fun Matrix<Double>.sum() = elements().map { (_, value) -> value }.sum()
+
+fun Matrix<Double>.min() = elements().map { (_, value) -> value }.min()
+
+fun Matrix<Double>.max() = elements().map { (_, value) -> value }.max()
+
+fun Matrix<Double>.average() = elements().map { (_, value) -> value }.average()

kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt

@@ -0,0 +1,16 @@
+package scientific.kmath.real
+
+import scientifik.kmath.real.average
+import scientifik.kmath.real.realMatrix
+import scientifik.kmath.real.sum
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class RealMatrixTest {
+    @Test
+    fun testSum() {
+        val m = realMatrix(10, 10) { i, j -> (i + j).toDouble() }
+        assertEquals(m.sum(), 900.0)
+        assertEquals(m.average(), 9.0)
+    }
+}

kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt

@@ -0,0 +1,36 @@
+package scientifik.kmath.linear
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class VectorTest {
+    @Test
+    fun testSum() {
+        val vector1 = RealVector(5) { it.toDouble() }
+        val vector2 = RealVector(5) { 5 - it.toDouble() }
+        val sum = vector1 + vector2
+        assertEquals(5.0, sum[2])
+    }
+
+    @Test
+    fun testVectorToMatrix() {
+        val vector = RealVector(5) { it.toDouble() }
+        val matrix = vector.asMatrix()
+        assertEquals(4.0, matrix[4, 0])
+    }
+
+    @Test
+    fun testDot() {
+        val vector1 = RealVector(5) { it.toDouble() }
+        val vector2 = RealVector(5) { 5 - it.toDouble() }
+
+        val matrix1 = vector1.asMatrix()
+        val matrix2 = vector2.asMatrix().transpose()
+        val product = MatrixContext.real.run { matrix1 dot matrix2 }
+
+
+        assertEquals(5.0, product[1, 0])
+        assertEquals(6.0, product[2, 2])
+    }
+
+}

kmath-functions/build.gradle.kts

@@ -0,0 +1,11 @@
+plugins {
+    id("scientifik.mpp")
+}
+
+kotlin.sourceSets {
+    commonMain {
+        dependencies {
+            api(project(":kmath-core"))
+        }
+    }
+}

kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Piecewise.kt

@@ -0,0 +1,57 @@
+package scientifik.kmath.functions
+
+import scientifik.kmath.operations.Ring
+
+interface Piecewise<T, R> {
+    fun findPiece(arg: T): R?
+}
+
+interface PiecewisePolynomial<T : Any> :
+    Piecewise<T, Polynomial<T>>
+
+/**
+ * Ordered list of pieces in piecewise function
+ */
+class OrderedPiecewisePolynomial<T : Comparable<T>>(delimeter: T) :
+    PiecewisePolynomial<T> {
+
+    private val delimiters: ArrayList<T> = arrayListOf(delimeter)
+    private val pieces: ArrayList<Polynomial<T>> = ArrayList()
+
+    /**
+     * Dynamically add a piece to the "right" side (beyond maximum argument value of previous piece)
+     * @param right new rightmost position. If is less then current rightmost position, a error is thrown.
+     */
+    fun putRight(right: T, piece: Polynomial<T>) {
+        require(right > delimiters.last()) { "New delimiter should be to the right of old one" }
+        delimiters.add(right)
+        pieces.add(piece)
+    }
+
+    fun putLeft(left: T, piece: Polynomial<T>) {
+        require(left < delimiters.first()) { "New delimiter should be to the left of old one" }
+        delimiters.add(0, left)
+        pieces.add(0, piece)
+    }
+
+    override fun findPiece(arg: T): Polynomial<T>? {
+        if (arg < delimiters.first() || arg >= delimiters.last()) {
+            return null
+        } else {
+            for (index in 1 until delimiters.size) {
+                if (arg < delimiters[index]) {
+                    return pieces[index - 1]
+                }
+            }
+            error("Piece not found")
+        }
+    }
+}
+
+/**
+ * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition.
+ */
+fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.value(ring: C, arg: T): T? =
+    findPiece(arg)?.value(ring, arg)
+
+fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { value(ring, it) }

kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/Polynomial.kt

@@ -0,0 +1,73 @@
+package scientifik.kmath.functions
+
+import scientifik.kmath.operations.Ring
+import scientifik.kmath.operations.Space
+import kotlin.math.max
+import kotlin.math.pow
+
+/**
+ * Polynomial coefficients without fixation on specific context they are applied to
+ * @param coefficients constant is the leftmost coefficient
+ */
+inline class Polynomial<T : Any>(val coefficients: List<T>) {
+    constructor(vararg coefficients: T) : this(coefficients.toList())
+}
+
+fun Polynomial<Double>.value() =
+    coefficients.reduceIndexed { index: Int, acc: Double, d: Double -> acc + d.pow(index) }
+
+
+fun <T : Any, C : Ring<T>> Polynomial<T>.value(ring: C, arg: T): T = ring.run {
+    if (coefficients.isEmpty()) return@run zero
+    var res = coefficients.first()
+    var powerArg = arg
+    for (index in 1 until coefficients.size) {
+        res += coefficients[index] * powerArg
+        //recalculating power on each step to avoid power costs on long polynomials
+        powerArg *= arg
+    }
+    return@run res
+}
+
+/**
+ * Represent a polynomial as a context-dependent function
+ */
+fun <T : Any, C : Ring<T>> Polynomial<T>.asMathFunction(): MathFunction<T, out C, T> = object :
+    MathFunction<T, C, T> {
+    override fun C.invoke(arg: T): T = value(this, arg)
+}
+
+/**
+ * Represent the polynomial as a regular context-less function
+ */
+fun <T : Any, C : Ring<T>> Polynomial<T>.asFunction(ring: C): (T) -> T = { value(ring, it) }
+
+/**
+ * An algebra for polynomials
+ */
+class PolynomialSpace<T : Any, C : Ring<T>>(val ring: C) : Space<Polynomial<T>> {
+
+    override fun add(a: Polynomial<T>, b: Polynomial<T>): Polynomial<T> {
+        val dim = max(a.coefficients.size, b.coefficients.size)
+        ring.run {
+            return Polynomial(List(dim) { index ->
+                a.coefficients.getOrElse(index) { zero } + b.coefficients.getOrElse(index) { zero }
+            })
+        }
+    }
+
+    override fun multiply(a: Polynomial<T>, k: Number): Polynomial<T> {
+        ring.run {
+            return Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * k })
+        }
+    }
+
+    override val zero: Polynomial<T> =
+        Polynomial(emptyList())
+
+    operator fun Polynomial<T>.invoke(arg: T): T = value(ring, arg)
+}
+
+fun <T : Any, C : Ring<T>, R> C.polynomial(block: PolynomialSpace<T, C>.() -> R): R {
+    return PolynomialSpace(this).run(block)
+}

kmath-functions/src/commonMain/kotlin/scientifik/kmath/functions/functions.kt

@@ -0,0 +1,33 @@
+package scientifik.kmath.functions
+
+import scientifik.kmath.operations.Algebra
+import scientifik.kmath.operations.RealField
+
+/**
+ * A regular function that could be called only inside specific algebra context
+ * @param T source type
+ * @param C source algebra constraint
+ * @param R result type
+ */
+interface MathFunction<T, C : Algebra<T>, R> {
+    operator fun C.invoke(arg: T): R
+}
+
+fun <R> MathFunction<Double, RealField, R>.invoke(arg: Double): R = RealField.invoke(arg)
+
+/**
+ * A suspendable function defined in algebraic context
+ */
+interface SuspendableMathFunction<T, C : Algebra<T>, R> {
+    suspend operator fun C.invoke(arg: T): R
+}
+
+suspend fun <R> SuspendableMathFunction<Double, RealField, R>.invoke(arg: Double) = RealField.invoke(arg)
+
+
+/**
+ * A parametric function with parameter
+ */
+interface ParametricFunction<T, P, C : Algebra<T>> {
+    operator fun C.invoke(arg: T, parameter: P): T
+}

kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/Interpolator.kt

@@ -0,0 +1,45 @@
+package scientifik.kmath.interpolation
+
+import scientifik.kmath.functions.PiecewisePolynomial
+import scientifik.kmath.functions.value
+import scientifik.kmath.operations.Ring
+import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.asBuffer
+
+interface Interpolator<X, Y> {
+    fun interpolate(points: XYPointSet<X, Y>): (X) -> Y
+}
+
+interface PolynomialInterpolator<T : Comparable<T>> : Interpolator<T, T> {
+    val algebra: Ring<T>
+
+    fun getDefaultValue(): T = error("Out of bounds")
+
+    fun interpolatePolynomials(points: XYPointSet<T, T>): PiecewisePolynomial<T>
+
+    override fun interpolate(points: XYPointSet<T, T>): (T) -> T = { x ->
+        interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue()
+    }
+}
+
+fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
+    x: Buffer<T>,
+    y: Buffer<T>
+): PiecewisePolynomial<T> {
+    val pointSet = BufferXYPointSet(x, y)
+    return interpolatePolynomials(pointSet)
+}
+
+fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
+    data: Map<T, T>
+): PiecewisePolynomial<T> {
+    val pointSet = BufferXYPointSet(data.keys.toList().asBuffer(), data.values.toList().asBuffer())
+    return interpolatePolynomials(pointSet)
+}
+
+fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
+    data: List<Pair<T, T>>
+): PiecewisePolynomial<T> {
+    val pointSet = BufferXYPointSet(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer())
+    return interpolatePolynomials(pointSet)
+}

kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt

@@ -0,0 +1,26 @@
+package scientifik.kmath.interpolation
+
+import scientifik.kmath.functions.OrderedPiecewisePolynomial
+import scientifik.kmath.functions.PiecewisePolynomial
+import scientifik.kmath.functions.Polynomial
+import scientifik.kmath.operations.Field
+
+/**
+ * Reference JVM implementation: https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/analysis/interpolation/LinearInterpolator.java
+ */
+class LinearInterpolator<T : Comparable<T>>(override val algebra: Field<T>) : PolynomialInterpolator<T> {
+
+    override fun interpolatePolynomials(points: XYPointSet<T, T>): PiecewisePolynomial<T> = algebra.run {
+        require(points.size > 0) { "Point array should not be empty" }
+        insureSorted(points)
+
+        OrderedPiecewisePolynomial(points.x[0]).apply {
+            for (i in 0 until points.size - 1) {
+                val slope = (points.y[i + 1] - points.y[i]) / (points.x[i + 1] - points.x[i])
+                val const = points.y[i] - slope * points.x[i]
+                val polynomial = Polynomial(const, slope)
+                putRight(points.x[i + 1], polynomial)
+            }
+        }
+    }
+}

kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LoessInterpolator.kt

@@ -0,0 +1,296 @@
+package scientifik.kmath.interpolation
+//
+//import scientifik.kmath.functions.PiecewisePolynomial
+//import scientifik.kmath.operations.Ring
+//import scientifik.kmath.structures.Buffer
+//import kotlin.math.abs
+//import kotlin.math.sqrt
+//
+//
+///**
+// * Original code: https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/LoessInterpolator.java
+// */
+//class LoessInterpolator<T : Comparable<T>>(override val algebra: Ring<T>) : PolynomialInterpolator<T> {
+//    /**
+//     * The bandwidth parameter: when computing the loess fit at
+//     * a particular point, this fraction of source points closest
+//     * to the current point is taken into account for computing
+//     * a least-squares regression.
+//     *
+//     *
+//     * A sensible value is usually 0.25 to 0.5.
+//     */
+//    private var bandwidth = 0.0
+//
+//    /**
+//     * The number of robustness iterations parameter: this many
+//     * robustness iterations are done.
+//     *
+//     *
+//     * A sensible value is usually 0 (just the initial fit without any
+//     * robustness iterations) to 4.
+//     */
+//    private var robustnessIters = 0
+//
+//    /**
+//     * If the median residual at a certain robustness iteration
+//     * is less than this amount, no more iterations are done.
+//     */
+//    private var accuracy = 0.0
+//
+//    /**
+//     * Constructs a new [LoessInterpolator]
+//     * with a bandwidth of [.DEFAULT_BANDWIDTH],
+//     * [.DEFAULT_ROBUSTNESS_ITERS] robustness iterations
+//     * and an accuracy of {#link #DEFAULT_ACCURACY}.
+//     * See [.LoessInterpolator] for an explanation of
+//     * the parameters.
+//     */
+//    fun LoessInterpolator() {
+//        bandwidth = DEFAULT_BANDWIDTH
+//        robustnessIters = DEFAULT_ROBUSTNESS_ITERS
+//        accuracy = DEFAULT_ACCURACY
+//    }
+//
+//    fun LoessInterpolator(bandwidth: Double, robustnessIters: Int) {
+//        this(bandwidth, robustnessIters, DEFAULT_ACCURACY)
+//    }
+//
+//    fun LoessInterpolator(bandwidth: Double, robustnessIters: Int, accuracy: Double) {
+//        if (bandwidth < 0 ||
+//            bandwidth > 1
+//        ) {
+//            throw OutOfRangeException(LocalizedFormats.BANDWIDTH, bandwidth, 0, 1)
+//        }
+//        this.bandwidth = bandwidth
+//        if (robustnessIters < 0) {
+//            throw NotPositiveException(LocalizedFormats.ROBUSTNESS_ITERATIONS, robustnessIters)
+//        }
+//        this.robustnessIters = robustnessIters
+//        this.accuracy = accuracy
+//    }
+//
+//    fun interpolate(
+//        xval: DoubleArray,
+//        yval: DoubleArray
+//    ): PolynomialSplineFunction {
+//        return SplineInterpolator().interpolate(xval, smooth(xval, yval))
+//    }
+//
+//    fun XYZPointSet<Double, Double, Double>.smooth(): XYPointSet<Double, Double> {
+//        checkAllFiniteReal(x)
+//        checkAllFiniteReal(y)
+//        checkAllFiniteReal(z)
+//        MathArrays.checkOrder(xval)
+//        if (size == 1) {
+//            return doubleArrayOf(y[0])
+//        }
+//        if (size == 2) {
+//            return doubleArrayOf(y[0], y[1])
+//        }
+//        val bandwidthInPoints = (bandwidth * size).toInt()
+//        if (bandwidthInPoints < 2) {
+//            throw NumberIsTooSmallException(
+//                LocalizedFormats.BANDWIDTH,
+//                bandwidthInPoints, 2, true
+//            )
+//        }
+//        val res = DoubleArray(size)
+//        val residuals = DoubleArray(size)
+//        val sortedResiduals = DoubleArray(size)
+//        val robustnessWeights = DoubleArray(size)
+//        // Do an initial fit and 'robustnessIters' robustness iterations.
+//        // This is equivalent to doing 'robustnessIters+1' robustness iterations
+//        // starting with all robustness weights set to 1.
+//        Arrays.fill(robustnessWeights, 1.0)
+//        for (iter in 0..robustnessIters) {
+//            val bandwidthInterval = intArrayOf(0, bandwidthInPoints - 1)
+//            // At each x, compute a local weighted linear regression
+//            for (i in 0 until size) {
+////                val x = x[i]
+//                // Find out the interval of source points on which
+//                // a regression is to be made.
+//                if (i > 0) {
+//                    updateBandwidthInterval(x, z, i, bandwidthInterval)
+//                }
+//                val ileft = bandwidthInterval[0]
+//                val iright = bandwidthInterval[1]
+//                // Compute the point of the bandwidth interval that is
+//                // farthest from x
+//                val edge: Int
+//                edge = if (x[i] - x[ileft] > x[iright] - x[i]) {
+//                    ileft
+//                } else {
+//                    iright
+//                }
+//                // Compute a least-squares linear fit weighted by
+//                // the product of robustness weights and the tricube
+//                // weight function.
+//                // See http://en.wikipedia.org/wiki/Linear_regression
+//                // (section "Univariate linear case")
+//                // and http://en.wikipedia.org/wiki/Weighted_least_squares
+//                // (section "Weighted least squares")
+//                var sumWeights = 0.0
+//                var sumX = 0.0
+//                var sumXSquared = 0.0
+//                var sumY = 0.0
+//                var sumXY = 0.0
+//                val denom: Double = abs(1.0 / (x[edge] - x[i]))
+//                for (k in ileft..iright) {
+//                    val xk = x[k]
+//                    val yk = y[k]
+//                    val dist = if (k < i) x - xk else xk - x[i]
+//                    val w = tricube(dist * denom) * robustnessWeights[k] * z[k]
+//                    val xkw = xk * w
+//                    sumWeights += w
+//                    sumX += xkw
+//                    sumXSquared += xk * xkw
+//                    sumY += yk * w
+//                    sumXY += yk * xkw
+//                }
+//                val meanX = sumX / sumWeights
+//                val meanY = sumY / sumWeights
+//                val meanXY = sumXY / sumWeights
+//                val meanXSquared = sumXSquared / sumWeights
+//                val beta: Double
+//                beta = if (sqrt(abs(meanXSquared - meanX * meanX)) < accuracy) {
+//                    0.0
+//                } else {
+//                    (meanXY - meanX * meanY) / (meanXSquared - meanX * meanX)
+//                }
+//                val alpha = meanY - beta * meanX
+//                res[i] = beta * x[i] + alpha
+//                residuals[i] = abs(y[i] - res[i])
+//            }
+//            // No need to recompute the robustness weights at the last
+//            // iteration, they won't be needed anymore
+//            if (iter == robustnessIters) {
+//                break
+//            }
+//            // Recompute the robustness weights.
+//            // Find the median residual.
+//            // An arraycopy and a sort are completely tractable here,
+//            // because the preceding loop is a lot more expensive
+//            java.lang.System.arraycopy(residuals, 0, sortedResiduals, 0, size)
+//            Arrays.sort(sortedResiduals)
+//            val medianResidual = sortedResiduals[size / 2]
+//            if (abs(medianResidual) < accuracy) {
+//                break
+//            }
+//            for (i in 0 until size) {
+//                val arg = residuals[i] / (6 * medianResidual)
+//                if (arg >= 1) {
+//                    robustnessWeights[i] = 0.0
+//                } else {
+//                    val w = 1 - arg * arg
+//                    robustnessWeights[i] = w * w
+//                }
+//            }
+//        }
+//        return res
+//    }
+//
+//    fun smooth(xval: DoubleArray, yval: DoubleArray): DoubleArray {
+//        if (xval.size != yval.size) {
+//            throw DimensionMismatchException(xval.size, yval.size)
+//        }
+//        val unitWeights = DoubleArray(xval.size)
+//        Arrays.fill(unitWeights, 1.0)
+//        return smooth(xval, yval, unitWeights)
+//    }
+//
+//    /**
+//     * Given an index interval into xval that embraces a certain number of
+//     * points closest to `xval[i-1]`, update the interval so that it
+//     * embraces the same number of points closest to `xval[i]`,
+//     * ignoring zero weights.
+//     *
+//     * @param xval Arguments array.
+//     * @param weights Weights array.
+//     * @param i Index around which the new interval should be computed.
+//     * @param bandwidthInterval a two-element array {left, right} such that:
+//     * `(left==0 or xval[i] - xval[left-1] > xval[right] - xval[i])`
+//     * and
+//     * `(right==xval.length-1 or xval[right+1] - xval[i] > xval[i] - xval[left])`.
+//     * The array will be updated.
+//     */
+//    private fun updateBandwidthInterval(
+//        xval: Buffer<Double>, weights: Buffer<Double>,
+//        i: Int,
+//        bandwidthInterval: IntArray
+//    ) {
+//        val left = bandwidthInterval[0]
+//        val right = bandwidthInterval[1]
+//        // The right edge should be adjusted if the next point to the right
+//        // is closer to xval[i] than the leftmost point of the current interval
+//        val nextRight = nextNonzero(weights, right)
+//        if (nextRight < xval.size && xval[nextRight] - xval[i] < xval[i] - xval[left]) {
+//            val nextLeft = nextNonzero(weights, bandwidthInterval[0])
+//            bandwidthInterval[0] = nextLeft
+//            bandwidthInterval[1] = nextRight
+//        }
+//    }
+//
+//    /**
+//     * Return the smallest index `j` such that
+//     * `j > i && (j == weights.length || weights[j] != 0)`.
+//     *
+//     * @param weights Weights array.
+//     * @param i Index from which to start search.
+//     * @return the smallest compliant index.
+//     */
+//    private fun nextNonzero(weights: Buffer<Double>, i: Int): Int {
+//        var j = i + 1
+//        while (j < weights.size && weights[j] == 0.0) {
+//            ++j
+//        }
+//        return j
+//    }
+//
+//    /**
+//     * Compute the
+//     * [tricube](http://en.wikipedia.org/wiki/Local_regression#Weight_function)
+//     * weight function
+//     *
+//     * @param x Argument.
+//     * @return `(1 - |x|<sup>3</sup>)<sup>3</sup>` for |x| &lt; 1, 0 otherwise.
+//     */
+//    private fun tricube(x: Double): Double {
+//        val absX: Double = FastMath.abs(x)
+//        if (absX >= 1.0) {
+//            return 0.0
+//        }
+//        val tmp = 1 - absX * absX * absX
+//        return tmp * tmp * tmp
+//    }
+//
+//    /**
+//     * Check that all elements of an array are finite real numbers.
+//     *
+//     * @param values Values array.
+//     * @throws org.apache.commons.math4.exception.NotFiniteNumberException
+//     * if one of the values is not a finite real number.
+//     */
+//    private fun checkAllFiniteReal(values: DoubleArray) {
+//        for (i in values.indices) {
+//            MathUtils.checkFinite(values[i])
+//        }
+//    }
+//
+//    override fun interpolatePolynomials(points: Collection<Pair<T, T>>): PiecewisePolynomial<T> {
+//        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
+//    }
+//
+//    companion object {
+//        /** Default value of the bandwidth parameter.  */
+//        const val DEFAULT_BANDWIDTH = 0.3
+//
+//        /** Default value of the number of robustness iterations.  */
+//        const val DEFAULT_ROBUSTNESS_ITERS = 2
+//
+//        /**
+//         * Default value for accuracy.
+//         */
+//        const val DEFAULT_ACCURACY = 1e-12
+//    }
+//}

kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/SplineInterpolator.kt

@@ -0,0 +1,58 @@
+package scientifik.kmath.interpolation
+
+import scientifik.kmath.functions.OrderedPiecewisePolynomial
+import scientifik.kmath.functions.PiecewisePolynomial
+import scientifik.kmath.functions.Polynomial
+import scientifik.kmath.operations.Field
+import scientifik.kmath.structures.MutableBufferFactory
+
+/**
+ * Generic spline interpolator. Not recommended for performance critical places, use platform-specific and type specific ones.
+ * Based on https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/SplineInterpolator.java
+ */
+class SplineInterpolator<T : Comparable<T>>(
+    override val algebra: Field<T>,
+    val bufferFactory: MutableBufferFactory<T>
+) : PolynomialInterpolator<T> {
+
+    //TODO possibly optimize zeroed buffers
+
+    override fun interpolatePolynomials(points: XYPointSet<T, T>): PiecewisePolynomial<T> = algebra.run {
+        if (points.size < 3) {
+            error("Can't use spline interpolator with less than 3 points")
+        }
+        insureSorted(points)
+
+        // Number of intervals.  The number of data points is n + 1.
+        val n = points.size - 1
+        // Differences between knot points
+        val h = bufferFactory(points.size) { i -> points.x[i + 1] - points.x[i] }
+        val mu = bufferFactory(points.size - 1) { zero }
+        val z = bufferFactory(points.size) { zero }
+
+        for (i in 1 until n) {
+            val g = 2.0 * (points.x[i + 1] - points.x[i - 1]) - h[i - 1] * mu[i - 1]
+            mu[i] = h[i] / g
+            z[i] =
+                (3.0 * (points.y[i + 1] * h[i - 1] - points.x[i] * (points.x[i + 1] - points.x[i - 1]) + points.y[i - 1] * h[i]) / (h[i - 1] * h[i])
+                        - h[i - 1] * z[i - 1]) / g
+        }
+
+        // cubic spline coefficients --  b is linear, c quadratic, d is cubic (original y's are constants)
+
+        OrderedPiecewisePolynomial<T>(points.x[points.size - 1]).apply {
+            var cOld = zero
+            for (j in n - 1 downTo 0) {
+                val c = z[j] - mu[j] * cOld
+                val a = points.y[j]
+                val b = (points.y[j + 1] - points.y[j]) / h[j] - h[j] * (cOld + 2.0 * c) / 3.0
+                val d = (cOld - c) / (3.0 * h[j])
+                val polynomial = Polynomial(a, b, c, d)
+                cOld = c
+                putLeft(points.x[j], polynomial)
+            }
+        }
+
+    }
+
+}

kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/XYPointSet.kt

@@ -0,0 +1,54 @@
+package scientifik.kmath.interpolation
+
+import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.Structure2D
+
+interface XYPointSet<X, Y> {
+    val size: Int
+    val x: Buffer<X>
+    val y: Buffer<Y>
+}
+
+interface XYZPointSet<X, Y, Z> : XYPointSet<X, Y> {
+    val z: Buffer<Z>
+}
+
+internal fun <T : Comparable<T>> insureSorted(points: XYPointSet<T, *>) {
+    for (i in 0 until points.size - 1) {
+        if (points.x[i + 1] <= points.x[i]) error("Input data is not sorted at index $i")
+    }
+}
+
+class NDStructureColumn<T>(val structure: Structure2D<T>, val column: Int) : Buffer<T> {
+    init {
+        require(column < structure.colNum) { "Column index is outside of structure column range" }
+    }
+
+    override val size: Int get() = structure.rowNum
+
+    override fun get(index: Int): T = structure[index, column]
+
+    override fun iterator(): Iterator<T> = sequence {
+        repeat(size) {
+            yield(get(it))
+        }
+    }.iterator()
+}
+
+class BufferXYPointSet<X, Y>(override val x: Buffer<X>, override val y: Buffer<Y>) : XYPointSet<X, Y> {
+    init {
+        require(x.size == y.size) { "Sizes of x and y buffers should be the same" }
+    }
+
+    override val size: Int
+        get() = x.size
+}
+
+fun <T> Structure2D<T>.asXYPointSet(): XYPointSet<T, T> {
+    require(shape[1] == 2) { "Structure second dimension should be of size 2" }
+    return object : XYPointSet<T, T> {
+        override val size: Int get() = this@asXYPointSet.shape[0]
+        override val x: Buffer<T> get() = NDStructureColumn(this@asXYPointSet, 0)
+        override val y: Buffer<T> get() = NDStructureColumn(this@asXYPointSet, 1)
+    }
+}

kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt

@@ -0,0 +1,27 @@
+package scientifik.kmath.interpolation
+
+import scientifik.kmath.functions.PiecewisePolynomial
+import scientifik.kmath.functions.asFunction
+import scientifik.kmath.operations.RealField
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+
+class LinearInterpolatorTest {
+    @Test
+    fun testInterpolation() {
+        val data = listOf(
+            0.0 to 0.0,
+            1.0 to 1.0,
+            2.0 to 3.0,
+            3.0 to 4.0
+        )
+        val polynomial: PiecewisePolynomial<Double> = LinearInterpolator(RealField).interpolatePolynomials(data)
+        val function = polynomial.asFunction(RealField)
+
+        assertEquals(null, function(-1.0))
+        assertEquals(0.5, function(0.5))
+        assertEquals(2.0, function(1.5))
+        assertEquals(3.0, function(2.0))
+    }
+}

kmath-geometry/build.gradle.kts

@@ -0,0 +1,9 @@
+plugins {
+    id("scientifik.mpp")
+}
+
+kotlin.sourceSets.commonMain {
+    dependencies {
+        api(project(":kmath-core"))
+    }
+}

kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Euclidean2DSpace.kt

@@ -0,0 +1,58 @@
+package scientifik.kmath.geometry
+
+import scientifik.kmath.linear.Point
+import scientifik.kmath.operations.SpaceElement
+import scientifik.kmath.operations.invoke
+import kotlin.math.sqrt
+
+
+interface Vector2D : Point<Double>, Vector, SpaceElement<Vector2D, Vector2D, Euclidean2DSpace> {
+    val x: Double
+    val y: Double
+
+    override val size: Int get() = 2
+
+    override fun get(index: Int): Double = when (index) {
+        1 -> x
+        2 -> y
+        else -> error("Accessing outside of point bounds")
+    }
+
+    override fun iterator(): Iterator<Double> = listOf(x, y).iterator()
+
+    override val context: Euclidean2DSpace get() = Euclidean2DSpace
+
+    override fun unwrap(): Vector2D = this
+
+    override fun Vector2D.wrap(): Vector2D = this
+}
+
+val Vector2D.r: Double get() = Euclidean2DSpace.run { sqrt(norm()) }
+
+@Suppress("FunctionName")
+fun Vector2D(x: Double, y: Double): Vector2D = Vector2DImpl(x, y)
+
+private data class Vector2DImpl(
+    override val x: Double,
+    override val y: Double
+) : Vector2D
+
+/**
+ * 2D Euclidean space
+ */
+object Euclidean2DSpace : GeometrySpace<Vector2D> {
+    fun Vector2D.norm(): Double = sqrt(x * x + y * y)
+
+    override fun Vector2D.distanceTo(other: Vector2D): Double = (this - other).norm()
+
+    override fun add(a: Vector2D, b: Vector2D): Vector2D =
+        Vector2D(a.x + b.x, a.y + b.y)
+
+    override fun multiply(a: Vector2D, k: Number): Vector2D =
+        Vector2D(a.x * k.toDouble(), a.y * k.toDouble())
+
+    override val zero: Vector2D = Vector2D(0.0, 0.0)
+
+    override fun Vector2D.dot(other: Vector2D): Double =
+        x * other.x + y * other.y
+}

kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Euclidean3DSpace.kt

@@ -0,0 +1,57 @@
+package scientifik.kmath.geometry
+
+import scientifik.kmath.linear.Point
+import scientifik.kmath.operations.SpaceElement
+import kotlin.math.sqrt
+
+
+interface Vector3D : Point<Double>, Vector, SpaceElement<Vector3D, Vector3D, Euclidean3DSpace> {
+    val x: Double
+    val y: Double
+    val z: Double
+
+    override val size: Int get() = 3
+
+    override fun get(index: Int): Double = when (index) {
+        1 -> x
+        2 -> y
+        3 -> z
+        else -> error("Accessing outside of point bounds")
+    }
+
+    override fun iterator(): Iterator<Double> = listOf(x, y, z).iterator()
+
+    override val context: Euclidean3DSpace get() = Euclidean3DSpace
+
+    override fun unwrap(): Vector3D = this
+
+    override fun Vector3D.wrap(): Vector3D = this
+}
+
+@Suppress("FunctionName")
+fun Vector3D(x: Double, y: Double, z: Double): Vector3D = Vector3DImpl(x, y, z)
+
+val Vector3D.r: Double get() = Euclidean3DSpace.run { sqrt(norm()) }
+
+private data class Vector3DImpl(
+    override val x: Double,
+    override val y: Double,
+    override val z: Double
+) : Vector3D
+
+object Euclidean3DSpace : GeometrySpace<Vector3D> {
+    override val zero: Vector3D = Vector3D(0.0, 0.0, 0.0)
+
+    fun Vector3D.norm(): Double = sqrt(x * x + y * y + z * z)
+
+    override fun Vector3D.distanceTo(other: Vector3D): Double = (this - other).norm()
+
+    override fun add(a: Vector3D, b: Vector3D): Vector3D =
+        Vector3D(a.x + b.x, a.y + b.y, a.z + b.z)
+
+    override fun multiply(a: Vector3D, k: Number): Vector3D =
+        Vector3D(a.x * k.toDouble(), a.y * k.toDouble(), a.z * k.toDouble())
+
+    override fun Vector3D.dot(other: Vector3D): Double =
+        x * other.x + y * other.y + z * other.z
+}

kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/GeometrySpace.kt

@@ -0,0 +1,17 @@
+package scientifik.kmath.geometry
+
+import scientifik.kmath.operations.Space
+
+interface Vector
+
+interface GeometrySpace<V: Vector>: Space<V> {
+    /**
+     * L2 distance
+     */
+    fun V.distanceTo(other: V): Double
+
+    /**
+     * Scalar product
+     */
+    infix fun V.dot(other: V): Double
+}

kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/Line.kt

@@ -0,0 +1,6 @@
+package scientifik.kmath.geometry
+
+data class Line<V: Vector>(val base: V, val direction: V)
+
+typealias Line2D = Line<Vector2D>
+typealias Line3D = Line<Vector3D>

kmath-geometry/src/commonMain/kotlin/scientifik/kmath/geometry/ReferenceFrame.kt

@@ -0,0 +1,4 @@
+package scientifik.kmath.geometry
+
+interface ReferenceFrame {
+}

kmath-histograms/build.gradle.kts

@@ -5,5 +5,6 @@ plugins {
 kotlin.sourceSets.commonMain {
     dependencies {
         api(project(":kmath-core"))
+        api(project(":kmath-for-real"))
     }
 }

kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt

@@ -45,8 +45,7 @@ class RealHistogram(
 
     private val values: NDStructure<LongCounter> = NDStructure.auto(strides) { LongCounter() }
 
-    //private val weight: NDStructure<DoubleCounter?> = ndStructure(strides){null}
+    private val weights: NDStructure<DoubleCounter> = NDStructure.auto(strides) { DoubleCounter() }
-
 
     override val dimension: Int get() = lower.size
 
@@ -102,21 +101,33 @@ class RealHistogram(
         return MultivariateBin(getDef(index), getValue(index))
     }
 
+//    fun put(point: Point<out Double>){
+//        val index = getIndex(point)
+//        values[index].increment()
+//    }
+
     override fun putWithWeight(point: Buffer<out Double>, weight: Double) {
-        if (weight != 1.0) TODO("Implement weighting")
         val index = getIndex(point)
         values[index].increment()
+        weights[index].add(weight)
     }
 
-    override fun iterator(): Iterator<MultivariateBin<Double>> = values.elements().map { (index, value) ->
+    override fun iterator(): Iterator<MultivariateBin<Double>> = weights.elements().map { (index, value) ->
         MultivariateBin(getDef(index), value.sum())
     }.iterator()
 
     /**
      * Convert this histogram into NDStructure containing bin values but not bin descriptions
      */
-    fun asNDStructure(): NDStructure<Number> {
+    fun values(): NDStructure<Number> {
-        return NDStructure.auto(this.values.shape) { values[it].sum() }
+        return NDStructure.auto(values.shape) { values[it].sum() }
+    }
+
+    /**
+     * Sum of weights
+     */
+    fun weights():NDStructure<Double>{
+        return NDStructure.auto(weights.shape) { weights[it].sum() }
     }
 
     companion object {

kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt

@@ -1,6 +1,10 @@
 package scientifik.memory
 
 import java.nio.ByteBuffer
+import java.nio.channels.FileChannel
+import java.nio.file.Files
+import java.nio.file.Path
+import java.nio.file.StandardOpenOption
 
 
 /**
@@ -11,7 +15,7 @@ actual fun Memory.Companion.allocate(length: Int): Memory {
     return ByteBufferMemory(buffer)
 }
 
-class ByteBufferMemory(
+private class ByteBufferMemory(
     val buffer: ByteBuffer,
     val startOffset: Int = 0,
     override val size: Int = buffer.limit()
@@ -91,3 +95,12 @@ class ByteBufferMemory(
 
     override fun writer(): MemoryWriter = writer
 }
+
+/**
+ * Use direct memory-mapped buffer from file to read something and close it afterwards.
+ */
+fun <R> Path.readAsMemory(position: Long = 0, size: Long = Files.size(this), block: Memory.() -> R): R {
+    return FileChannel.open(this, StandardOpenOption.READ).use {
+        ByteBufferMemory(it.map(FileChannel.MapMode.READ_ONLY, position, size)).block()
+    }
+}

kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Distribution.kt

@@ -58,9 +58,13 @@ fun <T : Any> Sampler<T>.sampleBuffer(
     //creating temporary storage once
     val tmp = ArrayList<T>(size)
     return sample(generator).collect { chain ->
-        for (i in tmp.indices) {
+        //clear list from previous run
-            tmp[i] = chain.next()
+        tmp.clear()
+        //Fill list
+        repeat(size){
+            tmp.add(chain.next())
         }
+        //return new buffer with elements from tmp
         bufferFactory(size) { tmp[it] }
     }
 }

kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomChain.kt

@@ -12,4 +12,3 @@ class RandomChain<out R>(val generator: RandomGenerator, private val gen: suspen
 }
 
 fun <R> RandomGenerator.chain(gen: suspend RandomGenerator.() -> R): RandomChain<R> = RandomChain(this, gen)
-fun <R> RandomGenerator.flow(gen: suspend RandomGenerator.() -> R) = chain(gen).fork()

kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Statistic.kt

@@ -11,6 +11,7 @@ import scientifik.kmath.coroutines.mapParallel
 import scientifik.kmath.operations.*
 import scientifik.kmath.structures.Buffer
 import scientifik.kmath.structures.asIterable
+import scientifik.kmath.structures.asSequence
 
 /**
  * A function, that transforms a buffer of random quantities to some resulting value
@@ -83,9 +84,9 @@ class Mean<T>(val space: Space<T>) : ComposableStatistic<T, Pair<T, Int>, T> {
 /**
  * Non-composable median
  */
-class Median<T>(comparator: Comparator<T>) : Statistic<T, T> {
+class Median<T>(private val comparator: Comparator<T>) : Statistic<T, T> {
     override suspend fun invoke(data: Buffer<T>): T {
-        return data.asIterable().toList()[data.size / 2] //TODO check if this is correct
+        return data.asSequence().sortedWith(comparator).toList()[data.size / 2] //TODO check if this is correct
     }
 
     companion object {

kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/SamplerTest.kt

@@ -0,0 +1,17 @@
+package scientifik.kmath.prob
+
+import kotlinx.coroutines.runBlocking
+import kotlin.test.Test
+
+class SamplerTest {
+
+    @Test
+    fun bufferSamplerTest(){
+        val sampler: Sampler<Double> =
+            BasicSampler { it.chain { nextDouble() } }
+        val data = sampler.sampleBuffer(RandomGenerator.default, 100)
+        runBlocking {
+            println(data.next())
+        }
+    }
+}

kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt

@@ -3,7 +3,7 @@ package scientifik.kmath.prob
 import kotlinx.coroutines.flow.drop
 import kotlinx.coroutines.flow.first
 import kotlinx.coroutines.runBlocking
-import scientifik.kmath.chains.flow
+
 import scientifik.kmath.streaming.chunked
 import kotlin.test.Test
 
@@ -13,7 +13,7 @@ class StatisticTest {
     //Create a stateless chain from generator.
     val data = generator.chain { nextDouble() }
     //Convert a chaint to Flow and break it into chunks.
-    val chunked = data.flow().chunked(1000)
+    val chunked = data.chunked(1000)
 
     @Test
     fun testParallelMean() {

kmath-viktor/build.gradle.kts

@@ -0,0 +1,10 @@
+plugins {
+    id("scientifik.jvm")
+}
+
+description = "Binding for https://github.com/JetBrains-Research/viktor"
+
+dependencies {
+    api(project(":kmath-core"))
+    api("org.jetbrains.bio:viktor:1.0.1")
+}

kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorBuffer.kt

@@ -0,0 +1,20 @@
+package scientifik.kmath.viktor
+
+import org.jetbrains.bio.viktor.F64FlatArray
+import scientifik.kmath.structures.MutableBuffer
+
+@Suppress("NOTHING_TO_INLINE", "OVERRIDE_BY_INLINE")
+inline class ViktorBuffer(val flatArray: F64FlatArray) : MutableBuffer<Double> {
+    override val size: Int get() = flatArray.size
+
+    override inline fun get(index: Int): Double = flatArray[index]
+    override inline fun set(index: Int, value: Double) {
+        flatArray[index] = value
+    }
+
+    override fun copy(): MutableBuffer<Double> {
+        return ViktorBuffer(flatArray.copy().flatten())
+    }
+
+    override fun iterator(): Iterator<Double> = flatArray.data.iterator()
+}

kmath-viktor/src/main/kotlin/scientifik/kmath/viktor/ViktorNDStructure.kt

@@ -0,0 +1,86 @@
+package scientifik.kmath.viktor
+
+import org.jetbrains.bio.viktor.F64Array
+import scientifik.kmath.operations.RealField
+import scientifik.kmath.structures.DefaultStrides
+import scientifik.kmath.structures.MutableNDStructure
+import scientifik.kmath.structures.NDField
+
+@Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
+inline class ViktorNDStructure(val f64Buffer: F64Array) : MutableNDStructure<Double> {
+
+    override val shape: IntArray get() = f64Buffer.shape
+
+    override inline fun get(index: IntArray): Double = f64Buffer.get(*index)
+
+    override inline fun set(index: IntArray, value: Double) {
+        f64Buffer.set(*index, value = value)
+    }
+
+    override fun elements(): Sequence<Pair<IntArray, Double>> {
+        return DefaultStrides(shape).indices().map { it to get(it) }
+    }
+}
+
+fun F64Array.asStructure(): ViktorNDStructure = ViktorNDStructure(this)
+
+@Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
+class ViktorNDField(override val shape: IntArray) : NDField<Double, RealField, ViktorNDStructure> {
+    override val zero: ViktorNDStructure
+        get() = F64Array.full(init = 0.0, shape = *shape).asStructure()
+    override val one: ViktorNDStructure
+        get() = F64Array.full(init = 1.0, shape = *shape).asStructure()
+
+    val strides = DefaultStrides(shape)
+
+    override val elementContext: RealField get() = RealField
+
+    override fun produce(initializer: RealField.(IntArray) -> Double): ViktorNDStructure = F64Array(*shape).apply {
+        this@ViktorNDField.strides.indices().forEach { index ->
+            set(value = RealField.initializer(index), indices = *index)
+        }
+    }.asStructure()
+
+    override fun map(arg: ViktorNDStructure, transform: RealField.(Double) -> Double): ViktorNDStructure =
+        F64Array(*shape).apply {
+            this@ViktorNDField.strides.indices().forEach { index ->
+                set(value = RealField.transform(arg[index]), indices = *index)
+            }
+        }.asStructure()
+
+    override fun mapIndexed(
+        arg: ViktorNDStructure,
+        transform: RealField.(index: IntArray, Double) -> Double
+    ): ViktorNDStructure = F64Array(*shape).apply {
+        this@ViktorNDField.strides.indices().forEach { index ->
+            set(value = RealField.transform(index, arg[index]), indices = *index)
+        }
+    }.asStructure()
+
+    override fun combine(
+        a: ViktorNDStructure,
+        b: ViktorNDStructure,
+        transform: RealField.(Double, Double) -> Double
+    ): ViktorNDStructure = F64Array(*shape).apply {
+        this@ViktorNDField.strides.indices().forEach { index ->
+            set(value = RealField.transform(a[index], b[index]), indices = *index)
+        }
+    }.asStructure()
+
+    override fun add(a: ViktorNDStructure, b: ViktorNDStructure): ViktorNDStructure {
+        return (a.f64Buffer + b.f64Buffer).asStructure()
+    }
+
+    override fun multiply(a: ViktorNDStructure, k: Number): ViktorNDStructure =
+        (a.f64Buffer * k.toDouble()).asStructure()
+
+    override inline fun ViktorNDStructure.plus(b: ViktorNDStructure): ViktorNDStructure =
+        (f64Buffer + b.f64Buffer).asStructure()
+
+    override inline fun ViktorNDStructure.minus(b: ViktorNDStructure): ViktorNDStructure =
+        (f64Buffer - b.f64Buffer).asStructure()
+
+    override inline fun ViktorNDStructure.times(k: Number): ViktorNDStructure = (f64Buffer * k.toDouble()).asStructure()
+
+    override inline fun ViktorNDStructure.plus(arg: Double): ViktorNDStructure = (f64Buffer.plus(arg)).asStructure()
+}

settings.gradle.kts

@@ -1,10 +1,10 @@
 pluginManagement {
 
     plugins {
-        id("scientifik.mpp") version "0.2.5"
+        id("scientifik.mpp") version "0.4.1"
-        id("scientifik.jvm") version "0.2.5"
+        id("scientifik.jvm") version "0.4.1"
-        id("scientifik.atomic") version "0.2.5"
+        id("scientifik.atomic") version "0.4.1"
-        id("scientifik.publish") version "0.2.5"
+        id("scientifik.publish") version "0.4.1"
     }
 
     repositories {
@@ -13,13 +13,14 @@ pluginManagement {
         gradlePluginPortal()
         maven("https://dl.bintray.com/kotlin/kotlin-eap")
         maven("https://dl.bintray.com/mipt-npm/scientifik")
+        maven("https://dl.bintray.com/mipt-npm/dev")
         maven("https://dl.bintray.com/kotlin/kotlinx")
     }
 
     resolutionStrategy {
         eachPlugin {
             when (requested.id.id) {
-                "scientifik.mpp", "scientifik.publish" -> useModule("scientifik:gradle-tools:${requested.version}")
+                "scientifik.mpp", "scientifik.jvm", "scientifik.publish" -> useModule("scientifik:gradle-tools:${requested.version}")
             }
         }
     }
@@ -29,12 +30,17 @@ rootProject.name = "kmath"
 include(
     ":kmath-memory",
     ":kmath-core",
+    ":kmath-functions",
 //    ":kmath-io",
     ":kmath-coroutines",
     ":kmath-histograms",
     ":kmath-commons",
+    ":kmath-viktor",
     ":kmath-koma",
     ":kmath-prob",
     ":kmath-io",
+    ":kmath-dimensions",
+    ":kmath-for-real",
+    ":kmath-geometry",
     ":examples"
 )