Merge branch 'dev' into beta/kotlin-2.0.20

# Conflicts: # benchmarks/build.gradle.kts # gradle.properties # kmath-ast/src/jsMain/kotlin/space/kscience/kmath/internal/base64/base64.kt
2024-09-22 16:46:45 +03:00 · 2024-09-22 16:46:45 +03:00 · 8a453bf0b9
commit 8a453bf0b9
parent 0a90d2e8c9 ab16bd16ac
176 changed files with 1851 additions and 1284 deletions
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@ -3,6 +3,7 @@
 ## Unreleased

 ### Added
+- Metropolis-Hastings sampler

 ### Changed

@ -200,7 +201,7 @@
 - Full autodiff refactoring based on `Symbol`
 - `kmath-prob` renamed to `kmath-stat`
 - Grid generators moved to `kmath-for-real`
- Use `Point<Double>` instead of specialized type in `kmath-for-real`
+- Use `Point<Float64>` instead of specialized type in `kmath-for-real`
 - Optimized dot product for buffer matrices moved to `kmath-for-real`
 - EjmlMatrix context is an object
 - Matrix LUP `inverse` renamed to `inverseWithLup`
--- a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attributes.kt
+++ b/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attributes.kt
@ -140,4 +140,18 @@ public fun <A : Attribute<Unit>> Attributes(
    attribute: A,
 ): Attributes = MapAttributes(mapOf(attribute to Unit))

-public operator fun Attributes.plus(other: Attributes): Attributes = MapAttributes(content + other.content)
+/**
+ * Create a new [Attributes] that overlays [other] on top of this set of attributes. New attributes are added.
+ * Existing attribute keys are replaced.
+ */
+public operator fun Attributes.plus(other: Attributes): Attributes = when {
+    isEmpty() -> other
+    other.isEmpty() -> this
+    else -> MapAttributes(content + other.content)
+}
+
+/**
+ * Create a new [Attributes] with removed [key] (if it is present).
+ */
+public operator fun Attributes.minus(key: Attribute<*>): Attributes =
+    if (content.contains(key)) MapAttributes(content.minus(key)) else this
--- a/benchmarks/README.md
+++ b/benchmarks/README.md
@ -1,4 +1,119 @@
-# Module benchmarks
+# BenchmarksResult
+
+## Report for benchmark configuration <code>main</code>
+
+* Run on Java HotSpot(TM) 64-Bit Server VM (build 21.0.4+8-LTS-jvmci-23.1-b41) with Java process:
+
+```
+C:\Users\altavir\scoop\apps\graalvm-oracle-21jdk\current\bin\java.exe -XX:ThreadPriorityPolicy=1 -XX:+UnlockExperimentalVMOptions -XX:+EnableJVMCIProduct -XX:-UnlockExperimentalVMOptions -Dfile.encoding=UTF-8 -Duser.country=US -Duser.language=en -Duser.variant
+```
+* JMH 1.21 was used in `thrpt` mode with 5 warmup iterations by 10 s and 5 measurement iterations by 10 s.
+### [ArrayBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ArrayBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`benchmarkArrayRead`|1.9E+07 &plusmn; 2.3E+05 ops/s|
+|`benchmarkBufferRead`|1.4E+07 &plusmn; 8.7E+05 ops/s|
+|`nativeBufferRead`|1.4E+07 &plusmn; 1.3E+06 ops/s|
+### [BigIntBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/BigIntBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`jvmAdd`|5.1E+07 &plusmn; 1.3E+06 ops/s|
+|`jvmAddLarge`|5.1E+04 &plusmn; 8.2E+02 ops/s|
+|`jvmMultiply`|8.5E+07 &plusmn; 9.7E+06 ops/s|
+|`jvmMultiplyLarge`|2.5E+02 &plusmn; 15 ops/s|
+|`jvmParsing10`|8.7E+06 &plusmn; 5.1E+05 ops/s|
+|`jvmParsing16`|6.4E+06 &plusmn; 1.8E+05 ops/s|
+|`jvmPower`|28 &plusmn; 0.79 ops/s|
+|`jvmSmallAdd`|7.0E+07 &plusmn; 4.3E+06 ops/s|
+|`kmAdd`|4.8E+07 &plusmn; 2.2E+06 ops/s|
+|`kmAddLarge`|3.5E+04 &plusmn; 3.7E+03 ops/s|
+|`kmMultiply`|6.7E+07 &plusmn; 1.5E+07 ops/s|
+|`kmMultiplyLarge`|54 &plusmn; 4.2 ops/s|
+|`kmParsing10`|4.5E+06 &plusmn; 8.3E+04 ops/s|
+|`kmParsing16`|4.9E+06 &plusmn; 1.1E+05 ops/s|
+|`kmPower`|10 &plusmn; 0.96 ops/s|
+|`kmSmallAdd`|4.1E+07 &plusmn; 5.9E+05 ops/s|
+### [BufferBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/BufferBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`bufferViewReadWrite`|5.8E+06 &plusmn; 1.6E+05 ops/s|
+|`bufferViewReadWriteSpecialized`|5.6E+06 &plusmn; 2.6E+05 ops/s|
+|`complexBufferReadWrite`|6.6E+06 &plusmn; 2.7E+05 ops/s|
+|`doubleArrayReadWrite`|7.5E+06 &plusmn; 1.0E+06 ops/s|
+|`doubleBufferReadWrite`|8.0E+06 &plusmn; 6.7E+05 ops/s|
+### [DotBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`bufferedDot`|1.3 &plusmn; 0.020 ops/s|
+|`cmDot`|0.47 &plusmn; 0.42 ops/s|
+|`cmDotWithConversion`|0.76 &plusmn; 0.13 ops/s|
+|`ejmlDot`|6.7 &plusmn; 0.091 ops/s|
+|`ejmlDotWithConversion`|6.4 &plusmn; 0.82 ops/s|
+|`multikDot`|40 &plusmn; 6.7 ops/s|
+|`parallelDot`|12 &plusmn; 1.8 ops/s|
+|`tensorDot`|1.2 &plusmn; 0.041 ops/s|
+|`tfDot`|5.9 &plusmn; 0.49 ops/s|
+### [ExpressionsInterpretersBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`asmGenericExpression`|29 &plusmn; 1.2 ops/s|
+|`asmPrimitiveExpression`|43 &plusmn; 1.3 ops/s|
+|`asmPrimitiveExpressionArray`|71 &plusmn; 0.38 ops/s|
+|`functionalExpression`|5.6 &plusmn; 0.11 ops/s|
+|`justCalculate`|69 &plusmn; 9.0 ops/s|
+|`mstExpression`|7.1 &plusmn; 0.020 ops/s|
+|`rawExpression`|41 &plusmn; 1.5 ops/s|
+### [IntegrationBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/IntegrationBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`complexIntegration`|3.6E+03 &plusmn; 1.9E+02 ops/s|
+|`doubleIntegration`|3.7E+03 &plusmn; 12 ops/s|
+### [JafamaBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/JafamaBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`core`|38 &plusmn; 0.64 ops/s|
+|`jafama`|52 &plusmn; 0.36 ops/s|
+|`strictJafama`|52 &plusmn; 4.0 ops/s|
+### [MatrixInverseBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`cmLUPInversion`|2.2E+03 &plusmn; 76 ops/s|
+|`ejmlInverse`|1.3E+03 &plusmn; 5.7 ops/s|
+|`kmathLupInversion`|9.5E+02 &plusmn; 1.8E+02 ops/s|
+|`kmathParallelLupInversion`|9.1E+02 &plusmn; 1.4E+02 ops/s|
+### [NDFieldBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/NDFieldBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`boxingFieldAdd`|7.7 &plusmn; 0.79 ops/s|
+|`multikAdd`|6.5 &plusmn; 0.33 ops/s|
+|`multikInPlaceAdd`|64 &plusmn; 0.79 ops/s|
+|`specializedFieldAdd`|8.0 &plusmn; 0.090 ops/s|
+|`tensorAdd`|9.2 &plusmn; 0.053 ops/s|
+|`tensorInPlaceAdd`|17 &plusmn; 10 ops/s|
+|`viktorAdd`|7.6 &plusmn; 1.2 ops/s|
+### [ViktorBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`doubleFieldAddition`|7.7 &plusmn; 0.34 ops/s|
+|`rawViktor`|5.9 &plusmn; 1.1 ops/s|
+|`viktorFieldAddition`|7.3 &plusmn; 1.1 ops/s|
+### [ViktorLogBenchmark](src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorLogBenchmark.kt)
+
+| Benchmark | Score |
+|:---------:|:-----:|
+|`rawViktorLog`|1.4 &plusmn; 0.076 ops/s|
+|`realFieldLog`|1.3 &plusmn; 0.069 ops/s|
+|`viktorFieldLog`|1.3 &plusmn; 0.032 ops/s|



--- a/benchmarks/build.gradle.kts
+++ b/benchmarks/build.gradle.kts
@ -1,13 +1,14 @@
-@file:Suppress("UNUSED_VARIABLE")
-
-import org.jetbrains.kotlin.gradle.dsl.JvmTarget
-import org.jetbrains.kotlin.gradle.tasks.KotlinJvmCompile
-import space.kscience.kmath.benchmarks.addBenchmarkProperties
+import com.fasterxml.jackson.module.kotlin.jacksonObjectMapper
+import com.fasterxml.jackson.module.kotlin.readValue
+import kotlinx.benchmark.gradle.BenchmarksExtension
+import java.time.LocalDateTime
+import java.time.ZoneId
+import java.util.*

 plugins {
    kotlin("multiplatform")
    alias(spclibs.plugins.kotlin.plugin.allopen)
-    id("org.jetbrains.kotlinx.benchmark")
+    alias(spclibs.plugins.kotlinx.benchmark)
 }

 allOpen.annotation("org.openjdk.jmh.annotations.State")
@ -17,8 +18,6 @@ repositories {
    mavenCentral()
 }

-val multikVersion: String by rootProject.extra
-
 kotlin {
    jvm()

@ -47,7 +46,7 @@ kotlin {
                implementation(project(":kmath-for-real"))
                implementation(project(":kmath-tensors"))
                implementation(project(":kmath-multik"))
-                implementation("org.jetbrains.kotlinx:multik-default:$multikVersion")
+                implementation(libs.multik.default)
                implementation(spclibs.kotlinx.benchmark.runtime)
            }
        }
@ -153,23 +152,130 @@ benchmark {
    }
 }

-kotlin.sourceSets.all {
-    with(languageSettings) {
-        optIn("kotlin.contracts.ExperimentalContracts")
-        optIn("kotlin.ExperimentalUnsignedTypes")
-        optIn("space.kscience.kmath.UnstableKMathAPI")
+kotlin {
+    jvmToolchain(11)
+    compilerOptions {
+        optIn.addAll(
+            "space.kscience.kmath.UnstableKMathAPI"
+        )
    }
 }

-tasks.withType<KotlinJvmCompile> {
-    compilerOptions {
-        jvmTarget.set(JvmTarget.JVM_11)
-        freeCompilerArgs.addAll("-Xjvm-default=all", "-Xlambdas=indy")
+
+private data class JmhReport(
+    val jmhVersion: String,
+    val benchmark: String,
+    val mode: String,
+    val threads: Int,
+    val forks: Int,
+    val jvm: String,
+    val jvmArgs: List<String>,
+    val jdkVersion: String,
+    val vmName: String,
+    val vmVersion: String,
+    val warmupIterations: Int,
+    val warmupTime: String,
+    val warmupBatchSize: Int,
+    val measurementIterations: Int,
+    val measurementTime: String,
+    val measurementBatchSize: Int,
+    val params: Map<String, String> = emptyMap(),
+    val primaryMetric: PrimaryMetric,
+    val secondaryMetrics: Map<String, SecondaryMetric>,
+) {
+    interface Metric {
+        val score: Double
+        val scoreError: Double
+        val scoreConfidence: List<Double>
+        val scorePercentiles: Map<Double, Double>
+        val scoreUnit: String
    }
+
+    data class PrimaryMetric(
+        override val score: Double,
+        override val scoreError: Double,
+        override val scoreConfidence: List<Double>,
+        override val scorePercentiles: Map<Double, Double>,
+        override val scoreUnit: String,
+        val rawDataHistogram: List<List<List<List<Double>>>>? = null,
+        val rawData: List<List<Double>>? = null,
+    ) : Metric
+
+    data class SecondaryMetric(
+        override val score: Double,
+        override val scoreError: Double,
+        override val scoreConfidence: List<Double>,
+        override val scorePercentiles: Map<Double, Double>,
+        override val scoreUnit: String,
+        val rawData: List<List<Double>>,
+    ) : Metric
 }

 readme {
    maturity = space.kscience.gradle.Maturity.EXPERIMENTAL
-}

-addBenchmarkProperties()
+    val jsonMapper = jacksonObjectMapper()
+
+    fun noun(number: Number, singular: String, plural: String) = if (number.toLong() == 1L) singular else plural
+
+    extensions.findByType(BenchmarksExtension::class.java)?.configurations?.forEach { cfg ->
+        val propertyName =
+            "benchmark${cfg.name.replaceFirstChar { if (it.isLowerCase()) it.titlecase(Locale.getDefault()) else it.toString() }}"
+
+        logger.info("Processing benchmark data from benchmark ${cfg.name} into readme property $propertyName")
+
+        val launches = layout.buildDirectory.dir("reports/benchmarks/${cfg.name}").get().asFile
+        if (!launches.exists()) return@forEach
+
+        property(propertyName) {
+            val resDirectory = launches.listFiles()?.maxByOrNull {
+                LocalDateTime.parse(it.name).atZone(ZoneId.systemDefault()).toInstant()
+            }
+
+            if (resDirectory == null || !(resDirectory.resolve("jvm.json")).exists()) {
+                "> **Can't find appropriate benchmark data. Try generating readme files after running benchmarks**."
+            } else {
+                val reports: List<JmhReport> =
+                    jsonMapper.readValue<List<JmhReport>>(resDirectory.resolve("jvm.json"))
+
+                buildString {
+                    appendLine("## Report for benchmark configuration <code>${cfg.name}</code>")
+                    appendLine()
+                    val first = reports.first()
+
+                    appendLine("* Run on ${first.vmName} (build ${first.vmVersion}) with Java process:")
+                    appendLine()
+                    appendLine("```")
+                    appendLine(
+                        "${first.jvm} ${
+                            first.jvmArgs.joinToString(" ")
+                        }"
+                    )
+                    appendLine("```")
+
+                    appendLine(
+                        "* JMH ${first.jmhVersion} was used in `${first.mode}` mode with ${first.warmupIterations} warmup ${
+                            noun(first.warmupIterations, "iteration", "iterations")
+                        } by ${first.warmupTime} and ${first.measurementIterations} measurement ${
+                            noun(first.measurementIterations, "iteration", "iterations")
+                        } by ${first.measurementTime}."
+                    )
+
+                    reports.groupBy { it.benchmark.substringBeforeLast(".") }.forEach { (cl, compare) ->
+                        appendLine("### [${cl.substringAfterLast(".")}](src/jvmMain/kotlin/${cl.replace(".","/")}.kt)")
+                        appendLine()
+                        appendLine("| Benchmark | Score |")
+                        appendLine("|:---------:|:-----:|")
+                        compare.forEach { report ->
+                            val benchmarkName = report.benchmark.substringAfterLast(".")
+                            val score = String.format("%.2G", report.primaryMetric.score)
+                            val error = String.format("%.2G", report.primaryMetric.scoreError)
+
+                            appendLine("|`$benchmarkName`|$score &plusmn; $error ${report.primaryMetric.scoreUnit}|")
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
--- a/benchmarks/docs/README-TEMPLATE.md
+++ b/benchmarks/docs/README-TEMPLATE.md
@ -0,0 +1,5 @@
+# BenchmarksResult
+
+${benchmarkMain}
+
+
--- a/benchmarks/src/jsMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt
+++ b/benchmarks/src/jsMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt
@ -14,6 +14,7 @@ import space.kscience.kmath.expressions.*
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.bindSymbol
 import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.Float64
 import kotlin.math.sin
 import kotlin.random.Random
 import space.kscience.kmath.estree.compileToExpression as estreeCompileToExpression
@ -67,7 +68,7 @@ class ExpressionsInterpretersBenchmark {
        blackhole.consume(sum)
    }

-    private fun invokeAndSum(expr: Expression<Double>, blackhole: Blackhole) {
+    private fun invokeAndSum(expr: Expression<Float64>, blackhole: Blackhole) {
        val random = Random(0)
        var sum = 0.0
        val m = HashMap<Symbol, Double>()
@ -99,7 +100,7 @@ class ExpressionsInterpretersBenchmark {
        private val wasm = node.wasmCompileToExpression(Float64Field)
        private val estree = node.estreeCompileToExpression(Float64Field)

-        private val raw = Expression<Double> { args ->
+        private val raw = Expression<Float64> { args ->
            val x = args.getValue(x)
            x * 2.0 + 2.0 / x - 16.0 / sin(x)
        }
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt
+++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt
@ -15,6 +15,7 @@ import space.kscience.kmath.operations.Algebra
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.bindSymbol
 import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.Float64
 import kotlin.math.sin
 import kotlin.random.Random

@ -83,7 +84,7 @@ internal class ExpressionsInterpretersBenchmark {
        blackhole.consume(sum)
    }

-    private fun invokeAndSum(expr: Expression<Double>, blackhole: Blackhole) {
+    private fun invokeAndSum(expr: Expression<Float64>, blackhole: Blackhole) {
        val random = Random(0)
        var sum = 0.0
        val m = HashMap<Symbol, Double>()
@ -114,9 +115,9 @@ internal class ExpressionsInterpretersBenchmark {
        private val asmPrimitive = node.compileToExpression(Float64Field)
        private val xIdx = asmPrimitive.indexer.indexOf(x)

-        private val asmGeneric = node.compileToExpression(Float64Field as Algebra<Double>)
+        private val asmGeneric = node.compileToExpression(Float64Field as Algebra<Float64>)

-        private val raw = Expression<Double> { args ->
+        private val raw = Expression<Float64> { args ->
            val x = args[x]!!
            x * 2.0 + 2.0 / x - 16.0 / sin(x)
        }
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/NDFieldBenchmark.kt
+++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/NDFieldBenchmark.kt
@ -17,6 +17,7 @@ import space.kscience.kmath.UnsafeKMathAPI
 import space.kscience.kmath.nd.*
 import space.kscience.kmath.nd4j.nd4j
 import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.tensors.core.DoubleTensor
 import space.kscience.kmath.tensors.core.one
 import space.kscience.kmath.tensors.core.tensorAlgebra
@ -37,28 +38,28 @@ internal class NDFieldBenchmark {

    @Benchmark
    fun specializedFieldAdd(blackhole: Blackhole) = with(specializedField) {
-        var res: StructureND<Double> = one(shape)
+        var res: StructureND<Float64> = one(shape)
        repeat(n) { res += 1.0 }
        blackhole.consume(res)
    }

    @Benchmark
    fun boxingFieldAdd(blackhole: Blackhole) = with(genericField) {
-        var res: StructureND<Double> = one(shape)
+        var res: StructureND<Float64> = one(shape)
        repeat(n) { res += 1.0 }
        blackhole.consume(res)
    }

    @Benchmark
    fun multikAdd(blackhole: Blackhole) = with(multikAlgebra) {
-        var res: StructureND<Double> = one(shape)
+        var res: StructureND<Float64> = one(shape)
        repeat(n) { res += 1.0 }
        blackhole.consume(res)
    }

    @Benchmark
    fun viktorAdd(blackhole: Blackhole) = with(viktorField) {
-        var res: StructureND<Double> = one(shape)
+        var res: StructureND<Float64> = one(shape)
        repeat(n) { res += 1.0 }
        blackhole.consume(res)
    }
@ -87,7 +88,7 @@ internal class NDFieldBenchmark {

 //    @Benchmark
 //    fun nd4jAdd(blackhole: Blackhole) = with(nd4jField) {
-//        var res: StructureND<Double> = one(dim, dim)
+//        var res: StructureND<Float64> = one(dim, dim)
 //        repeat(n) { res += 1.0 }
 //        blackhole.consume(res)
 //    }
--- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorBenchmark.kt
+++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorBenchmark.kt
@ -15,6 +15,7 @@ import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.nd.ndAlgebra
 import space.kscience.kmath.nd.one
 import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.viktor.ViktorFieldND

@State(Scope.Benchmark)
@ -23,7 +24,7 @@ internal class ViktorBenchmark {
    @Benchmark
    fun doubleFieldAddition(blackhole: Blackhole) {
        with(doubleField) {
-            var res: StructureND<Double> = one(shape)
+            var res: StructureND<Float64> = one(shape)
            repeat(n) { res += 1.0 }
            blackhole.consume(res)
        }
--- a/build.gradle.kts
+++ b/build.gradle.kts
@ -2,8 +2,8 @@ import space.kscience.gradle.useApache2Licence
 import space.kscience.gradle.useSPCTeam

 plugins {
-    id("space.kscience.gradle.project")
-    id("org.jetbrains.kotlinx.kover") version "0.7.6"
+    alias(spclibs.plugins.kscience.project)
+    alias(spclibs.plugins.kotlinx.kover)
 }

 val attributesVersion by extra("0.2.0")
@ -70,5 +70,3 @@ ksciencePublish {
 }

 apiValidation.nonPublicMarkers.add("space.kscience.kmath.UnstableKMathAPI")
-
-val multikVersion by extra("0.2.3")
--- a/buildSrc/build.gradle.kts
+++ b/buildSrc/build.gradle.kts
@ -1,34 +0,0 @@
-plugins {
-    `kotlin-dsl`
-    `version-catalog`
-}
-
-repositories {
-    mavenLocal()
-    maven("https://repo.kotlin.link")
-    mavenCentral()
-    gradlePluginPortal()
-}
-
-val toolsVersion = spclibs.versions.tools.get()
-val kotlinVersion = spclibs.versions.kotlin.asProvider().get()
-val benchmarksVersion = spclibs.versions.kotlinx.benchmark.get()
-
-dependencies {
-    api("space.kscience:gradle-tools:$toolsVersion")
-    //plugins form benchmarks
-    api("org.jetbrains.kotlinx:kotlinx-benchmark-plugin:$benchmarksVersion")
-    //api("org.jetbrains.kotlin:kotlin-allopen:$kotlinVersion")
-    //to be used inside build-script only
-    //implementation(spclibs.kotlinx.serialization.json)
-    implementation("com.fasterxml.jackson.module:jackson-module-kotlin:2.14.+")
-}
-
-kotlin {
-    jvmToolchain {
-        languageVersion.set(JavaLanguageVersion.of(11))
-    }
-    sourceSets.all {
-        languageSettings.optIn("kotlin.OptIn")
-    }
-}
--- a/buildSrc/settings.gradle.kts
+++ b/buildSrc/settings.gradle.kts
@ -1,38 +0,0 @@
-/*
- * Copyright 2018-2021 KMath contributors.
- * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
- */
-
-enableFeaturePreview("TYPESAFE_PROJECT_ACCESSORS")
-
-plugins {
-    id("org.gradle.toolchains.foojay-resolver-convention") version "0.8.0"
-}
-
-dependencyResolutionManagement {
-    val projectProperties = java.util.Properties()
-    file("../gradle.properties").inputStream().use {
-        projectProperties.load(it)
-    }
-
-    projectProperties.forEach { key, value ->
-        extra.set(key.toString(), value)
-    }
-
-
-    val toolsVersion: String = projectProperties["toolsVersion"].toString()
-
-    @Suppress("UnstableApiUsage")
-    repositories {
-        mavenLocal()
-        maven("https://repo.kotlin.link")
-        mavenCentral()
-        gradlePluginPortal()
-    }
-
-    versionCatalogs {
-        create("spclibs") {
-            from("space.kscience:version-catalog:$toolsVersion")
-        }
-    }
-}
--- a/buildSrc/src/main/kotlin/space/kscience/kmath/benchmarks/JmhReport.kt
+++ b/buildSrc/src/main/kotlin/space/kscience/kmath/benchmarks/JmhReport.kt
@ -1,55 +0,0 @@
-/*
- * Copyright 2018-2024 KMath contributors.
- * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
- */
-
-package space.kscience.kmath.benchmarks
-
-data class JmhReport(
-    val jmhVersion: String,
-    val benchmark: String,
-    val mode: String,
-    val threads: Int,
-    val forks: Int,
-    val jvm: String,
-    val jvmArgs: List<String>,
-    val jdkVersion: String,
-    val vmName: String,
-    val vmVersion: String,
-    val warmupIterations: Int,
-    val warmupTime: String,
-    val warmupBatchSize: Int,
-    val measurementIterations: Int,
-    val measurementTime: String,
-    val measurementBatchSize: Int,
-    val params: Map<String, String> = emptyMap(),
-    val primaryMetric: PrimaryMetric,
-    val secondaryMetrics: Map<String, SecondaryMetric>,
-) {
-    interface Metric {
-        val score: Double
-        val scoreError: Double
-        val scoreConfidence: List<Double>
-        val scorePercentiles: Map<Double, Double>
-        val scoreUnit: String
-    }
-
-    data class PrimaryMetric(
-        override val score: Double,
-        override val scoreError: Double,
-        override val scoreConfidence: List<Double>,
-        override val scorePercentiles: Map<Double, Double>,
-        override val scoreUnit: String,
-        val rawDataHistogram: List<List<List<List<Double>>>>? = null,
-        val rawData: List<List<Double>>? = null,
-    ) : Metric
-
-    data class SecondaryMetric(
-        override val score: Double,
-        override val scoreError: Double,
-        override val scoreConfidence: List<Double>,
-        override val scorePercentiles: Map<Double, Double>,
-        override val scoreUnit: String,
-        val rawData: List<List<Double>>,
-    ) : Metric
-}
--- a/buildSrc/src/main/kotlin/space/kscience/kmath/benchmarks/addBenchmarkProperties.kt
+++ b/buildSrc/src/main/kotlin/space/kscience/kmath/benchmarks/addBenchmarkProperties.kt
@ -1,110 +0,0 @@
-/*
- * Copyright 2018-2024 KMath contributors.
- * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
- */
-
-package space.kscience.kmath.benchmarks
-
-import com.fasterxml.jackson.module.kotlin.jacksonObjectMapper
-import com.fasterxml.jackson.module.kotlin.readValue
-import kotlinx.benchmark.gradle.BenchmarksExtension
-import org.gradle.api.Project
-import space.kscience.gradle.KScienceReadmeExtension
-import java.time.LocalDateTime
-import java.time.ZoneId
-import java.time.format.DateTimeFormatter
-import java.time.format.DateTimeFormatterBuilder
-import java.time.format.SignStyle
-import java.time.temporal.ChronoField.*
-import java.util.*
-
-private val ISO_DATE_TIME: DateTimeFormatter = DateTimeFormatterBuilder().run {
-    parseCaseInsensitive()
-    appendValue(YEAR, 4, 10, SignStyle.EXCEEDS_PAD)
-    appendLiteral('-')
-    appendValue(MONTH_OF_YEAR, 2)
-    appendLiteral('-')
-    appendValue(DAY_OF_MONTH, 2)
-    appendLiteral('T')
-    appendValue(HOUR_OF_DAY, 2)
-    appendLiteral('.')
-    appendValue(MINUTE_OF_HOUR, 2)
-    optionalStart()
-    appendLiteral('.')
-    appendValue(SECOND_OF_MINUTE, 2)
-    optionalStart()
-    appendFraction(NANO_OF_SECOND, 0, 9, true)
-    optionalStart()
-    appendOffsetId()
-    optionalStart()
-    appendLiteral('[')
-    parseCaseSensitive()
-    appendZoneRegionId()
-    appendLiteral(']')
-    toFormatter()
-}
-
-private fun noun(number: Number, singular: String, plural: String) = if (number.toLong() == 1L) singular else plural
-
-private val jsonMapper = jacksonObjectMapper()
-
-fun Project.addBenchmarkProperties() {
-    val benchmarksProject = this
-    rootProject.subprojects.forEach { p ->
-        p.extensions.findByType(KScienceReadmeExtension::class.java)?.run {
-            benchmarksProject.extensions.findByType(BenchmarksExtension::class.java)?.configurations?.forEach { cfg ->
-                property("benchmark${cfg.name.replaceFirstChar { if (it.isLowerCase()) it.titlecase(Locale.getDefault()) else it.toString() }}") {
-                    val launches = benchmarksProject.layout.buildDirectory.dir("reports/benchmarks/${cfg.name}").get()
-
-                    val resDirectory = launches.files().maxByOrNull {
-                        LocalDateTime.parse(it.name, ISO_DATE_TIME).atZone(ZoneId.systemDefault()).toInstant()
-                    }
-
-                    if (resDirectory == null || !(resDirectory.resolve("jvm.json")).exists()) {
-                        "> **Can't find appropriate benchmark data. Try generating readme files after running benchmarks**."
-                    } else {
-                        val reports: List<JmhReport> =
-                            jsonMapper.readValue<List<JmhReport>>(resDirectory.resolve("jvm.json"))
-
-                        buildString {
-                            appendLine("<details>")
-                            appendLine("<summary>")
-                            appendLine("Report for benchmark configuration <code>${cfg.name}</code>")
-                            appendLine("</summary>")
-                            appendLine()
-                            val first = reports.first()
-
-                            appendLine("* Run on ${first.vmName} (build ${first.vmVersion}) with Java process:")
-                            appendLine()
-                            appendLine("```")
-                            appendLine(
-                                "${first.jvm} ${
-                                    first.jvmArgs.joinToString(" ")
-                                }"
-                            )
-                            appendLine("```")
-
-                            appendLine(
-                                "* JMH ${first.jmhVersion} was used in `${first.mode}` mode with ${first.warmupIterations} warmup ${
-                                    noun(first.warmupIterations, "iteration", "iterations")
-                                } by ${first.warmupTime} and ${first.measurementIterations} measurement ${
-                                    noun(first.measurementIterations, "iteration", "iterations")
-                                } by ${first.measurementTime}."
-                            )
-
-                            appendLine()
-                            appendLine("| Benchmark | Score |")
-                            appendLine("|:---------:|:-----:|")
-
-                            reports.forEach { report ->
-                                appendLine("|`${report.benchmark}`|${report.primaryMetric.score} &plusmn; ${report.primaryMetric.scoreError} ${report.primaryMetric.scoreUnit}|")
-                            }
-
-                            appendLine("</details>")
-                        }
-                    }
-                }
-            }
-        }
-    }
-}
--- a/docs/nd-structure.md
+++ b/docs/nd-structure.md
@ -31,7 +31,7 @@ The code to run this looks like:

 ```kotlin
    specializedField.run {
-        var res: NDBuffer<Double> = one
+        var res: NDBuffer<Float64> = one
        repeat(n) {
            res += 1.0
        }
@ -103,7 +103,7 @@ The boxing field produced by

 ```kotlin
    genericField.run {
-        var res: NDBuffer<Double> = one
+        var res: NDBuffer<Float64> = one
        repeat(n) {
            res += 1.0
        }
--- a/examples/build.gradle.kts
+++ b/examples/build.gradle.kts
@ -10,8 +10,6 @@ repositories {
    maven("https://maven.pkg.jetbrains.space/kotlin/p/kotlin/kotlin-js-wrappers")
 }

-val multikVersion: String by rootProject.extra
-
 dependencies {
    implementation(project(":kmath-ast"))
    implementation(project(":kmath-kotlingrad"))
@ -33,7 +31,7 @@ dependencies {
    implementation(project(":kmath-jafama"))
    //multik
    implementation(project(":kmath-multik"))
-    implementation("org.jetbrains.kotlinx:multik-default:$multikVersion")
+    implementation(libs.multik.default)

    //datetime
    implementation("org.jetbrains.kotlinx:kotlinx-datetime:0.4.0")
--- a/examples/src/main/kotlin/space/kscience/kmath/functions/integrate.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/functions/integrate.kt
@ -14,11 +14,12 @@ import space.kscience.kmath.integration.gaussIntegrator
 import space.kscience.kmath.integration.integrate
 import space.kscience.kmath.integration.value
 import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
 import kotlin.math.pow

 fun main() {
    //Define a function
-    val function: Function1D<Double> = { x -> 3 * x.pow(2) + 2 * x + 1 }
+    val function: Function1D<Float64> = { x -> 3 * x.pow(2) + 2 * x + 1 }

    //get the result of the integration
    val result = Float64Field.gaussIntegrator.integrate(0.0..10.0, function = function)
--- a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt
@ -8,6 +8,7 @@ package space.kscience.kmath.functions
 import space.kscience.kmath.interpolation.SplineInterpolator
 import space.kscience.kmath.interpolation.interpolatePolynomials
 import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
 import space.kscience.plotly.Plotly
 import space.kscience.plotly.UnstablePlotlyAPI
 import space.kscience.plotly.makeFile
@ -23,7 +24,7 @@ fun main() {
        x to sin(x)
    }

-    val polynomial: PiecewisePolynomial<Double> = SplineInterpolator(Float64Field).interpolatePolynomials(data)
+    val polynomial: PiecewisePolynomial<Float64> = SplineInterpolator(Float64Field).interpolatePolynomials(data)

    val function = polynomial.asFunction(Float64Field, 0.0)

--- a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt
@ -10,6 +10,7 @@ import space.kscience.kmath.interpolation.splineInterpolator
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
+import space.kscience.kmath.structures.Float64
 import space.kscience.plotly.Plotly
 import space.kscience.plotly.UnstablePlotlyAPI
 import space.kscience.plotly.makeFile
@ -18,7 +19,7 @@ import space.kscience.plotly.scatter

@OptIn(UnstablePlotlyAPI::class)
 fun main() {
-    val function: Function1D<Double> = { x ->
+    val function: Function1D<Float64> = { x ->
        if (x in 30.0..50.0) {
            1.0
        } else {
@ -28,7 +29,7 @@ fun main() {
    val xs = 0.0..100.0 step 0.5
    val ys = xs.map(function)

-    val polynomial: PiecewisePolynomial<Double> = Float64Field.splineInterpolator.interpolatePolynomials(xs, ys)
+    val polynomial: PiecewisePolynomial<Float64> = Float64Field.splineInterpolator.interpolatePolynomials(xs, ys)

    val polyFunction = polynomial.asFunction(Float64Field, 0.0)

--- a/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt
@ -12,6 +12,7 @@ import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.nd.structureND
 import space.kscience.kmath.nd.withNdAlgebra
 import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.structures.Float64
 import kotlin.math.pow

 fun main(): Unit = Double.algebra.withNdAlgebra(2, 2) {
@ -22,7 +23,7 @@ fun main(): Unit = Double.algebra.withNdAlgebra(2, 2) {
    }

    //Define a function in a nd space
-    val function: (Double) -> StructureND<Double> = { x: Double -> 3 * x.pow(2) + 2 * diagonal(x) + 1 }
+    val function: (Double) -> StructureND<Float64> = { x: Double -> 3 * x.pow(2) + 2 * diagonal(x) + 1 }

    //get the result of the integration
    val result = gaussIntegrator.integrate(0.0..10.0, function = function)
--- a/examples/src/main/kotlin/space/kscience/kmath/linear/eigenValueDecomposition.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/linear/eigenValueDecomposition.kt
@ -0,0 +1,40 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linear
+
+import space.kscience.kmath.commons.linear.CMLinearSpace
+import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.structures.Float64
+import kotlin.random.Random
+
+fun main() {
+    val dim = 46
+
+    val random = Random(123)
+
+    val u = Float64.algebra.linearSpace.buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+
+    listOf(CMLinearSpace, EjmlLinearSpaceDDRM).forEach { algebra ->
+        with(algebra) {
+            //create a simmetric matrix
+            val matrix = buildMatrix(dim, dim) { row, col ->
+                if (row >= col) u[row, col] else u[col, row]
+            }
+            val eigen = matrix.getOrComputeAttribute(EIG) ?: error("Failed to compute eigenvalue decomposition")
+            check(
+                StructureND.contentEquals(
+                    matrix,
+                    eigen.v dot eigen.d dot eigen.v.transposed(),
+                    1e-4
+                )
+            ) { "$algebra decomposition failed" }
+            println("$algebra eigenvalue decomposition complete and checked" )
+        }
+    }
+
+}
--- a/examples/src/main/kotlin/space/kscience/kmath/linear/gradient.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/linear/gradient.kt
@ -6,18 +6,19 @@
 package space.kscience.kmath.linear

 import space.kscience.kmath.real.*
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer

 fun main() {
    val x0 = DoubleVector(0.0, 0.0, 0.0)
    val sigma = DoubleVector(1.0, 1.0, 1.0)

-    val gaussian: (Point<Double>) -> Double = { x ->
+    val gaussian: (Point<Float64>) -> Double = { x ->
        require(x.size == x0.size)
        kotlin.math.exp(-((x - x0) / sigma).square().sum())
    }

-    fun ((Point<Double>) -> Double).grad(x: Point<Double>): Point<Double> {
+    fun ((Point<Float64>) -> Double).grad(x: Point<Float64>): Point<Float64> {
        require(x.size == x0.size)
        return Float64Buffer(x.size) { i ->
            val h = sigma[i] / 5
--- a/examples/src/main/kotlin/space/kscience/kmath/operations/mixedNDOperations.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/operations/mixedNDOperations.kt
@ -11,6 +11,7 @@ import space.kscience.kmath.nd.Float64BufferND
 import space.kscience.kmath.nd.Structure2D
 import space.kscience.kmath.nd.mutableStructureND
 import space.kscience.kmath.nd.ndAlgebra
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.viktor.viktorAlgebra
 import kotlin.collections.component1
 import kotlin.collections.component2
@ -20,7 +21,7 @@ fun main() {
        if (i == j) 2.0 else 0.0
    }

-    val cmMatrix: Structure2D<Double> = CMLinearSpace.matrix(2, 2)(0.0, 1.0, 0.0, 3.0)
+    val cmMatrix: Structure2D<Float64> = CMLinearSpace.matrix(2, 2)(0.0, 1.0, 0.0, 3.0)

    val res: Float64BufferND = Float64Field.ndAlgebra {
        exp(viktorStructure) + 2.0 * cmMatrix
--- a/examples/src/main/kotlin/space/kscience/kmath/series/analyzeDif.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/series/analyzeDif.kt
@ -14,6 +14,7 @@ import space.kscience.kmath.operations.toList
 import space.kscience.kmath.stat.KMComparisonResult
 import space.kscience.kmath.stat.ksComparisonStatistic
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.slice
 import space.kscience.plotly.*
 import kotlin.math.PI
@ -24,7 +25,7 @@ fun Double.Companion.seriesAlgebra() = Double.algebra.bufferAlgebra.seriesAlgebr
 fun main() = with(Double.seriesAlgebra()) {


-    fun Plot.plotSeries(name: String, buffer: Buffer<Double>) {
+    fun Plot.plotSeries(name: String, buffer: Buffer<Float64>) {
        scatter {
            this.name = name
            x.numbers = buffer.labels
@ -37,10 +38,10 @@ fun main() = with(Double.seriesAlgebra()) {

    val s2 = s1.slice(20..50).moveTo(40)

-    val s3: Buffer<Double> = s1.zip(s2) { l, r -> l + r } //s1 + s2
+    val s3: Buffer<Float64> = s1.zip(s2) { l, r -> l + r } //s1 + s2
    val s4 = s3.map { ln(it) }

-    val kmTest: KMComparisonResult<Double> = ksComparisonStatistic(s1, s2)
+    val kmTest: KMComparisonResult<Float64> = ksComparisonStatistic(s1, s2)

    Plotly.page {
        h1 { +"This is my plot" }
--- a/examples/src/main/kotlin/space/kscience/kmath/series/seriesBuilder.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/series/seriesBuilder.kt
@ -6,10 +6,7 @@
 package space.kscience.kmath.series


-import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.Float64Buffer
-import space.kscience.kmath.structures.asBuffer
-import space.kscience.kmath.structures.toDoubleArray
+import space.kscience.kmath.structures.*
 import space.kscience.plotly.*
 import space.kscience.plotly.models.Scatter
 import space.kscience.plotly.models.ScatterMode
@ -22,7 +19,7 @@ fun main(): Unit = with(Double.seriesAlgebra()) {
    val arrayOfRandoms = DoubleArray(20) { random.nextDouble() }

    val series1: Float64Buffer = arrayOfRandoms.asBuffer()
-    val series2: Series<Double> = series1.moveBy(3)
+    val series2: Series<Float64> = series1.moveBy(3)

    val res = series2 - series1

@ -30,7 +27,7 @@ fun main(): Unit = with(Double.seriesAlgebra()) {

    println(res)

-    fun Plot.series(name: String, buffer: Buffer<Double>, block: Scatter.() -> Unit = {}) {
+    fun Plot.series(name: String, buffer: Buffer<Float64>, block: Scatter.() -> Unit = {}) {
        scatter {
            this.name = name
            x.numbers = buffer.offsetIndices
--- a/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionBenchmark.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionBenchmark.kt
@ -36,8 +36,28 @@ private suspend fun runKMathChained(): Duration {
    return Duration.between(startTime, Instant.now())
 }

+private fun runKMathBlocking(): Duration {
+    val generator = RandomGenerator.fromSource(RandomSource.MT, 123L)
+    val normal = GaussianSampler(7.0, 2.0)
+    val chain = normal.sample(generator)
+    val startTime = Instant.now()
+    var sum = 0.0
+
+    repeat(10000001) { counter ->
+        sum += chain.nextBlocking()
+
+        if (counter % 100000 == 0) {
+            val duration = Duration.between(startTime, Instant.now())
+            val meanValue = sum / counter
+            println("Chain sampler completed $counter elements in $duration: $meanValue")
+        }
+    }
+
+    return Duration.between(startTime, Instant.now())
+}
+
 private fun runCMDirect(): Duration {
-    val rng = RandomSource.create(RandomSource.MT, 123L)
+    val rng = RandomSource.MT.create(123L)

    val sampler = CMGaussianSampler.of(
        BoxMullerNormalizedGaussianSampler.of(rng),
@ -67,6 +87,8 @@ private fun runCMDirect(): Duration {
 fun main(): Unit = runBlocking(Dispatchers.Default) {
    val directJob = async { runCMDirect() }
    val chainJob = async { runKMathChained() }
+    val blockingJob = async { runKMathBlocking() }
    println("KMath Chained: ${chainJob.await()}")
+    println("KMath Blocking: ${blockingJob.await()}")
    println("Apache Direct: ${directJob.await()}")
 }
--- a/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionDemo.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/stat/DistributionDemo.kt
@ -10,13 +10,14 @@ import space.kscience.kmath.chains.Chain
 import space.kscience.kmath.chains.combineWithState
 import space.kscience.kmath.distributions.NormalDistribution
 import space.kscience.kmath.random.RandomGenerator
+import space.kscience.kmath.structures.Float64

 private data class AveragingChainState(var num: Int = 0, var value: Double = 0.0)

 /**
 * Averaging.
 */
-private fun Chain<Double>.mean(): Chain<Double> = combineWithState(AveragingChainState(), { it.copy() }) { chain ->
+private fun Chain<Float64>.mean(): Chain<Float64> = combineWithState(AveragingChainState(), { it.copy() }) { chain ->
    val next = chain.next()
    num++
    value += next
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/ComplexND.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/ComplexND.kt
@ -26,7 +26,7 @@ fun main() {

    val realTime = measureTimeMillis {
        realField {
-            var res: StructureND<Double> = one
+            var res: StructureND<Float64> = one
            repeat(n) {
                res += 1.0
            }
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/NDField.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/NDField.kt
@ -45,35 +45,35 @@ fun main() {

    measureAndPrint("Boxing addition") {
        genericField {
-            var res: StructureND<Double> = one(shape)
+            var res: StructureND<Float64> = one(shape)
            repeat(n) { res += 1.0 }
        }
    }

    measureAndPrint("Specialized addition") {
        doubleField {
-            var res: StructureND<Double> = one(shape)
+            var res: StructureND<Float64> = one(shape)
            repeat(n) { res += 1.0 }
        }
    }

    measureAndPrint("Nd4j specialized addition") {
        nd4jField {
-            var res: StructureND<Double> = one(shape)
+            var res: StructureND<Float64> = one(shape)
            repeat(n) { res += 1.0 }
        }
    }

    measureAndPrint("Viktor addition") {
        viktorField {
-            var res: StructureND<Double> = one
+            var res: StructureND<Float64> = one
            repeat(n) { res += 1.0 }
        }
    }

    measureAndPrint("Parallel stream addition") {
        parallelField {
-            var res: StructureND<Double> = one
+            var res: StructureND<Float64> = one
            repeat(n) { res += 1.0 }
        }
    }
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/StreamDoubleFieldND.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/StreamDoubleFieldND.kt
@ -19,21 +19,21 @@ import java.util.stream.IntStream
 * execution.
 */
 class StreamDoubleFieldND(override val shape: ShapeND) : FieldND<Double, Float64Field>,
-    NumbersAddOps<StructureND<Double>>,
-    ExtendedField<StructureND<Double>> {
+    NumbersAddOps<StructureND<Float64>>,
+    ExtendedField<StructureND<Float64>> {

    private val strides = ColumnStrides(shape)
    override val elementAlgebra: Float64Field get() = Float64Field
-    override val zero: BufferND<Double> by lazy { structureND(shape) { zero } }
-    override val one: BufferND<Double> by lazy { structureND(shape) { one } }
+    override val zero: BufferND<Float64> by lazy { structureND(shape) { zero } }
+    override val one: BufferND<Float64> by lazy { structureND(shape) { one } }

-    override fun number(value: Number): BufferND<Double> {
+    override fun number(value: Number): BufferND<Float64> {
        val d = value.toDouble() // minimize conversions
        return structureND(shape) { d }
    }

    @OptIn(PerformancePitfall::class)
-    private val StructureND<Double>.buffer: Float64Buffer
+    private val StructureND<Float64>.buffer: Float64Buffer
        get() = when {
            shape != this@StreamDoubleFieldND.shape -> throw ShapeMismatchException(
                this@StreamDoubleFieldND.shape,
@ -44,7 +44,7 @@ class StreamDoubleFieldND(override val shape: ShapeND) : FieldND<Double, Float64
            else -> Float64Buffer(strides.linearSize) { offset -> get(strides.index(offset)) }
        }

-    override fun structureND(shape: ShapeND, initializer: Float64Field.(IntArray) -> Double): BufferND<Double> {
+    override fun structureND(shape: ShapeND, initializer: Float64Field.(IntArray) -> Double): BufferND<Float64> {
        val array = IntStream.range(0, strides.linearSize).parallel().mapToDouble { offset ->
            val index = strides.index(offset)
            Float64Field.initializer(index)
@ -56,7 +56,7 @@ class StreamDoubleFieldND(override val shape: ShapeND) : FieldND<Double, Float64
    override fun mutableStructureND(
        shape: ShapeND,
        initializer: DoubleField.(IntArray) -> Double,
-    ): MutableBufferND<Double> {
+    ): MutableBufferND<Float64> {
        val array = IntStream.range(0, strides.linearSize).parallel().mapToDouble { offset ->
            val index = strides.index(offset)
            DoubleField.initializer(index)
@ -66,17 +66,17 @@ class StreamDoubleFieldND(override val shape: ShapeND) : FieldND<Double, Float64
    }

    @OptIn(PerformancePitfall::class)
-    override fun StructureND<Double>.map(
+    override fun StructureND<Float64>.map(
        transform: Float64Field.(Double) -> Double,
-    ): BufferND<Double> {
+    ): BufferND<Float64> {
        val array = Arrays.stream(buffer.array).parallel().map { Float64Field.transform(it) }.toArray()
        return BufferND(strides, array.asBuffer())
    }

    @OptIn(PerformancePitfall::class)
-    override fun StructureND<Double>.mapIndexed(
+    override fun StructureND<Float64>.mapIndexed(
        transform: Float64Field.(index: IntArray, Double) -> Double,
-    ): BufferND<Double> {
+    ): BufferND<Float64> {
        val array = IntStream.range(0, strides.linearSize).parallel().mapToDouble { offset ->
            Float64Field.transform(
                strides.index(offset),
@ -89,39 +89,39 @@ class StreamDoubleFieldND(override val shape: ShapeND) : FieldND<Double, Float64

    @OptIn(PerformancePitfall::class)
    override fun zip(
-        left: StructureND<Double>,
-        right: StructureND<Double>,
+        left: StructureND<Float64>,
+        right: StructureND<Float64>,
        transform: Float64Field.(Double, Double) -> Double,
-    ): BufferND<Double> {
+    ): BufferND<Float64> {
        val array = IntStream.range(0, strides.linearSize).parallel().mapToDouble { offset ->
            Float64Field.transform(left.buffer.array[offset], right.buffer.array[offset])
        }.toArray()
        return BufferND(strides, array.asBuffer())
    }

-    override fun StructureND<Double>.unaryMinus(): StructureND<Double> = map { -it }
+    override fun StructureND<Float64>.unaryMinus(): StructureND<Float64> = map { -it }

-    override fun scale(a: StructureND<Double>, value: Double): StructureND<Double> = a.map { it * value }
+    override fun scale(a: StructureND<Float64>, value: Double): StructureND<Float64> = a.map { it * value }

-    override fun power(arg: StructureND<Double>, pow: Number): BufferND<Double> = arg.map { power(it, pow) }
+    override fun power(arg: StructureND<Float64>, pow: Number): BufferND<Float64> = arg.map { power(it, pow) }

-    override fun exp(arg: StructureND<Double>): BufferND<Double> = arg.map { exp(it) }
+    override fun exp(arg: StructureND<Float64>): BufferND<Float64> = arg.map { exp(it) }

-    override fun ln(arg: StructureND<Double>): BufferND<Double> = arg.map { ln(it) }
+    override fun ln(arg: StructureND<Float64>): BufferND<Float64> = arg.map { ln(it) }

-    override fun sin(arg: StructureND<Double>): BufferND<Double> = arg.map { sin(it) }
-    override fun cos(arg: StructureND<Double>): BufferND<Double> = arg.map { cos(it) }
-    override fun tan(arg: StructureND<Double>): BufferND<Double> = arg.map { tan(it) }
-    override fun asin(arg: StructureND<Double>): BufferND<Double> = arg.map { asin(it) }
-    override fun acos(arg: StructureND<Double>): BufferND<Double> = arg.map { acos(it) }
-    override fun atan(arg: StructureND<Double>): BufferND<Double> = arg.map { atan(it) }
+    override fun sin(arg: StructureND<Float64>): BufferND<Float64> = arg.map { sin(it) }
+    override fun cos(arg: StructureND<Float64>): BufferND<Float64> = arg.map { cos(it) }
+    override fun tan(arg: StructureND<Float64>): BufferND<Float64> = arg.map { tan(it) }
+    override fun asin(arg: StructureND<Float64>): BufferND<Float64> = arg.map { asin(it) }
+    override fun acos(arg: StructureND<Float64>): BufferND<Float64> = arg.map { acos(it) }
+    override fun atan(arg: StructureND<Float64>): BufferND<Float64> = arg.map { atan(it) }

-    override fun sinh(arg: StructureND<Double>): BufferND<Double> = arg.map { sinh(it) }
-    override fun cosh(arg: StructureND<Double>): BufferND<Double> = arg.map { cosh(it) }
-    override fun tanh(arg: StructureND<Double>): BufferND<Double> = arg.map { tanh(it) }
-    override fun asinh(arg: StructureND<Double>): BufferND<Double> = arg.map { asinh(it) }
-    override fun acosh(arg: StructureND<Double>): BufferND<Double> = arg.map { acosh(it) }
-    override fun atanh(arg: StructureND<Double>): BufferND<Double> = arg.map { atanh(it) }
+    override fun sinh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { sinh(it) }
+    override fun cosh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { cosh(it) }
+    override fun tanh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { tanh(it) }
+    override fun asinh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { asinh(it) }
+    override fun acosh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { acosh(it) }
+    override fun atanh(arg: StructureND<Float64>): BufferND<Float64> = arg.map { atanh(it) }
 }

 fun Float64Field.ndStreaming(vararg shape: Int): StreamDoubleFieldND = StreamDoubleFieldND(ShapeND(shape))
--- a/examples/src/main/kotlin/space/kscience/kmath/structures/mutableNd.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/structures/mutableNd.kt
@ -11,7 +11,7 @@ import space.kscience.kmath.operations.algebra

@OptIn(PerformancePitfall::class)
 fun main(): Unit = with(Double.algebra.ndAlgebra) {
-    val structure: MutableStructure2D<Double> = mutableStructureND(ShapeND(2, 2)) { (i, j) ->
+    val structure: MutableStructure2D<Float64> = mutableStructureND(ShapeND(2, 2)) { (i, j) ->
        i.toDouble() + j.toDouble()
    }.as2D()

--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/StreamingLm/streamLm.kt
@ -12,6 +12,7 @@ import space.kscience.kmath.nd.MutableStructure2D
 import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.nd.as2D
 import space.kscience.kmath.nd.component1
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.tensors.LevenbergMarquardt.StartDataLm
 import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.zeros
 import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
@ -20,9 +21,9 @@ import space.kscience.kmath.tensors.core.levenbergMarquardt
 import kotlin.random.Random

 fun streamLm(
-    lm_func: (MutableStructure2D<Double>, MutableStructure2D<Double>, Int) -> (MutableStructure2D<Double>),
+    lm_func: (MutableStructure2D<Float64>, MutableStructure2D<Float64>, Int) -> (MutableStructure2D<Float64>),
    startData: StartDataLm, launchFrequencyInMs: Long, numberOfLaunches: Int,
-): Flow<MutableStructure2D<Double>> = flow {
+): Flow<MutableStructure2D<Float64>> = flow {

    var example_number = startData.example_number
    var p_init = startData.p_init
@ -64,7 +65,7 @@ fun streamLm(
    }
 }

-fun generateNewYDat(y_dat: MutableStructure2D<Double>, delta: Double): MutableStructure2D<Double> {
+fun generateNewYDat(y_dat: MutableStructure2D<Float64>, delta: Double): MutableStructure2D<Float64> {
    val n = y_dat.shape.component1()
    val y_dat_new = zeros(ShapeND(intArrayOf(n, 1))).as2D()
    for (i in 0 until n) {
--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/LevenbergMarquardt/functionsToOptimize.kt
@ -9,6 +9,7 @@ import space.kscience.kmath.nd.MutableStructure2D
 import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.nd.as2D
 import space.kscience.kmath.nd.component1
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra
 import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra.div
 import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
@ -19,24 +20,24 @@ import space.kscience.kmath.tensors.core.DoubleTensorAlgebra.Companion.times
 import space.kscience.kmath.tensors.core.asDoubleTensor

 public data class StartDataLm(
-    var lm_matx_y_dat: MutableStructure2D<Double>,
+    var lm_matx_y_dat: MutableStructure2D<Float64>,
    var example_number: Int,
-    var p_init: MutableStructure2D<Double>,
-    var t: MutableStructure2D<Double>,
-    var y_dat: MutableStructure2D<Double>,
+    var p_init: MutableStructure2D<Float64>,
+    var t: MutableStructure2D<Float64>,
+    var y_dat: MutableStructure2D<Float64>,
    var weight: Double,
-    var dp: MutableStructure2D<Double>,
-    var p_min: MutableStructure2D<Double>,
-    var p_max: MutableStructure2D<Double>,
-    var consts: MutableStructure2D<Double>,
+    var dp: MutableStructure2D<Float64>,
+    var p_min: MutableStructure2D<Float64>,
+    var p_max: MutableStructure2D<Float64>,
+    var consts: MutableStructure2D<Float64>,
    var opts: DoubleArray,
 )

 fun funcEasyForLm(
-    t: MutableStructure2D<Double>,
-    p: MutableStructure2D<Double>,
+    t: MutableStructure2D<Float64>,
+    p: MutableStructure2D<Float64>,
    exampleNumber: Int,
-): MutableStructure2D<Double> {
+): MutableStructure2D<Float64> {
    val m = t.shape.component1()
    var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1)))

@ -59,10 +60,10 @@ fun funcEasyForLm(
 }

 fun funcMiddleForLm(
-    t: MutableStructure2D<Double>,
-    p: MutableStructure2D<Double>,
+    t: MutableStructure2D<Float64>,
+    p: MutableStructure2D<Float64>,
    exampleNumber: Int,
-): MutableStructure2D<Double> {
+): MutableStructure2D<Float64> {
    val m = t.shape.component1()
    var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1)))

@ -79,10 +80,10 @@ fun funcMiddleForLm(
 }

 fun funcDifficultForLm(
-    t: MutableStructure2D<Double>,
-    p: MutableStructure2D<Double>,
+    t: MutableStructure2D<Float64>,
+    p: MutableStructure2D<Float64>,
    exampleNumber: Int,
-): MutableStructure2D<Double> {
+): MutableStructure2D<Float64> {
    val m = t.shape.component1()
    var y_hat = DoubleTensorAlgebra.zeros(ShapeND(intArrayOf(m, 1)))

--- a/examples/src/main/kotlin/space/kscience/kmath/tensors/multik.kt
+++ b/examples/src/main/kotlin/space/kscience/kmath/tensors/multik.kt
@ -10,12 +10,13 @@ import org.jetbrains.kotlinx.multik.api.ndarray
 import org.jetbrains.kotlinx.multik.default.DefaultEngine
 import space.kscience.kmath.multik.MultikDoubleAlgebra
 import space.kscience.kmath.nd.one
+import space.kscience.kmath.structures.Float64


 val multikAlgebra = MultikDoubleAlgebra(DefaultEngine())

 fun main(): Unit = with(multikAlgebra) {
-    val a = Multik.ndarray(intArrayOf(1, 2, 3)).asType<Double>().wrap()
+    val a = Multik.ndarray(intArrayOf(1, 2, 3)).asType<Float64>().wrap()
    val b = Multik.ndarray(doubleArrayOf(1.0, 2.0, 3.0)).wrap()
    one(a.shape) - a + b * 3.0
 }
--- a/gradle.properties
+++ b/gradle.properties
@ -9,6 +9,5 @@ org.gradle.configureondemand=true
 org.gradle.jvmargs=-Xmx4096m
 org.gradle.parallel=true
 org.gradle.workers.max=4
-toolsVersion=0.15.5-kotlin-2.0.20-Beta2
-#kotlin.experimental.tryK2=true
-#kscience.wasm.disabled=true
+
+toolsVersion=0.15.4-kotlin-2.0.0
--- a/gradle/libs.versions.toml
+++ b/gradle/libs.versions.toml
@ -0,0 +1,14 @@
+[versions]
+
+commons-rng = "1.6"
+multik = "0.2.3"
+
+[libraries]
+
+commons-rng-simple = { module = "org.apache.commons:commons-rng-simple", version.ref = "commons-rng" }
+commons-rng-sampling = { module = "org.apache.commons:commons-rng-sampling", version.ref = "commons-rng" }
+
+multik-core = { module = "org.jetbrains.kotlinx:multik-core", version.ref = "multik" }
+multik-default = { module = "org.jetbrains.kotlinx:multik-default", version.ref = "multik" }
+
+[plugins]
--- a/kmath-ast/src/commonTest/kotlin/space/kscience/kmath/ast/TestFolding.kt
+++ b/kmath-ast/src/commonTest/kotlin/space/kscience/kmath/ast/TestFolding.kt
@ -9,6 +9,7 @@ import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.Int32Ring
 import space.kscience.kmath.operations.Int8Ring
 import space.kscience.kmath.operations.pi
+import space.kscience.kmath.structures.Float64
 import kotlin.test.Test
 import kotlin.test.assertEquals
 import kotlin.test.fail
@ -41,7 +42,7 @@ internal class TestFolding {
    @Test
    fun foldSymbol() = assertEquals(
        Float64Field.pi,
-        ("pi".parseMath().evaluateConstants(Float64Field) as? TypedMst.Constant<Double> ?: fail()).value,
+        ("pi".parseMath().evaluateConstants(Float64Field) as? TypedMst.Constant<Float64> ?: fail()).value,
    )

    @Test
--- a/kmath-ast/src/commonTest/kotlin/space/kscience/kmath/ast/utils.kt
+++ b/kmath-ast/src/commonTest/kotlin/space/kscience/kmath/ast/utils.kt
@ -10,12 +10,13 @@ import space.kscience.kmath.expressions.MST
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.Int32Ring
+import space.kscience.kmath.structures.Float64

 internal interface CompilerTestContext {
    fun MST.compileToExpression(algebra: Int32Ring): Expression<Int>
    fun MST.compile(algebra: Int32Ring, arguments: Map<Symbol, Int>): Int
    fun MST.compile(algebra: Int32Ring, vararg arguments: Pair<Symbol, Int>): Int = compile(algebra, mapOf(*arguments))
-    fun MST.compileToExpression(algebra: Float64Field): Expression<Double>
+    fun MST.compileToExpression(algebra: Float64Field): Expression<Float64>
    fun MST.compile(algebra: Float64Field, arguments: Map<Symbol, Double>): Double

    fun MST.compile(algebra: Float64Field, vararg arguments: Pair<Symbol, Double>): Double =
--- a/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/internal/base64/base64.kt
+++ b/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/internal/base64/base64.kt
@ -4,10 +4,6 @@
 */

@file:Suppress(
-    "INTERFACE_WITH_SUPERCLASS",
-    "OVERRIDING_FINAL_MEMBER",
-    "RETURN_TYPE_MISMATCH_ON_OVERRIDE",
-    "CONFLICTING_OVERLOADS",
    "ObjectPropertyName",
    "ClassName",
 )
--- a/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/internal/binaryen/index.binaryen.kt
+++ b/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/internal/binaryen/index.binaryen.kt
@ -4,10 +4,6 @@
 */

@file:Suppress(
-    "INTERFACE_WITH_SUPERCLASS",
-    "OVERRIDING_FINAL_MEMBER",
-    "RETURN_TYPE_MISMATCH_ON_OVERRIDE",
-    "CONFLICTING_OVERLOADS",
    "PropertyName",
    "ClassName", "ENUM_CLASS_IN_EXTERNAL_DECLARATION_WARNING",
 )
--- a/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/internal/webassembly/lib.dom.WebAssembly.module_dukat.kt
+++ b/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/internal/webassembly/lib.dom.WebAssembly.module_dukat.kt
@ -6,10 +6,6 @@
@file:JsQualifier("WebAssembly")

@file:Suppress(
-    "INTERFACE_WITH_SUPERCLASS",
-    "OVERRIDING_FINAL_MEMBER",
-    "RETURN_TYPE_MISMATCH_ON_OVERRIDE",
-    "NO_EXPLICIT_VISIBILITY_IN_API_MODE_WARNING",
    "ClassName",
 )

@ -22,78 +18,34 @@ import org.w3c.fetch.Response
 import space.kscience.kmath.internal.tsstdlib.PromiseLike
 import kotlin.js.Promise

-@Suppress("NESTED_CLASS_IN_EXTERNAL_INTERFACE")
-internal external interface CompileError {
-    companion object {
-        var prototype: CompileError
-    }
-}
+internal external interface CompileError 

-@Suppress("NESTED_CLASS_IN_EXTERNAL_INTERFACE")
 internal external interface Global {
    var value: Any
    fun valueOf(): Any
-
-    companion object {
-        var prototype: Global
-    }
 }

-@Suppress("NESTED_CLASS_IN_EXTERNAL_INTERFACE")
@JsName("Instance")
 internal external interface Instance1 {
    var exports: Exports
-
-    companion object {
-        var prototype: Instance
-    }
 }

-@Suppress("NESTED_CLASS_IN_EXTERNAL_INTERFACE")
-internal external interface LinkError {
-    companion object {
-        var prototype: LinkError
-    }
-}
+internal external interface LinkError

-@Suppress("NESTED_CLASS_IN_EXTERNAL_INTERFACE")
 internal external interface Memory {
    var buffer: ArrayBuffer
    fun grow(delta: Number): Number
-
-    companion object {
-        var prototype: Memory
-    }
 }

-@Suppress("NESTED_CLASS_IN_EXTERNAL_INTERFACE")
-@JsName("Module")
-internal external interface Module1 {
-    companion object {
-        var prototype: Module
-        fun customSections(moduleObject: Module, sectionName: String): Array<ArrayBuffer>
-        fun exports(moduleObject: Module): Array<ModuleExportDescriptor>
-        fun imports(moduleObject: Module): Array<ModuleImportDescriptor>
-    }
-}

-@Suppress("NESTED_CLASS_IN_EXTERNAL_INTERFACE")
-internal external interface RuntimeError {
-    companion object {
-        var prototype: RuntimeError
-    }
-}
+internal external interface RuntimeError
+

-@Suppress("NESTED_CLASS_IN_EXTERNAL_INTERFACE")
 internal external interface Table {
    var length: Number
    fun get(index: Number): Function<*>?
    fun grow(delta: Number): Number
    fun set(index: Number, value: Function<*>?)
-
-    companion object {
-        var prototype: Table
-    }
 }

 internal external interface GlobalDescriptor {
--- a/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/internal/webassembly/nonDeclarations.WebAssembly.kt
+++ b/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/internal/webassembly/nonDeclarations.WebAssembly.kt
@ -3,13 +3,6 @@
 * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
 */

-@file:Suppress(
-    "INTERFACE_WITH_SUPERCLASS",
-    "OVERRIDING_FINAL_MEMBER",
-    "RETURN_TYPE_MISMATCH_ON_OVERRIDE",
-    "CONFLICTING_OVERLOADS",
-    "NO_EXPLICIT_VISIBILITY_IN_API_MODE_WARNING",
-)

 package space.kscience.kmath.internal.webassembly

--- a/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/wasm/internal/WasmBuilder.kt
+++ b/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/wasm/internal/WasmBuilder.kt
@ -11,6 +11,7 @@ import space.kscience.kmath.expressions.*
 import space.kscience.kmath.internal.binaryen.*
 import space.kscience.kmath.internal.webassembly.Instance
 import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.internal.binaryen.Module as BinaryenModule
 import space.kscience.kmath.internal.webassembly.Module as WasmModule

@ -85,13 +86,13 @@ internal sealed class WasmBuilder<T : Number, out E : Expression<T>>(
 }

@UnstableKMathAPI
-internal class DoubleWasmBuilder(target: TypedMst<Double>) :
+internal class DoubleWasmBuilder(target: TypedMst<Float64>) :
    WasmBuilder<Double, DoubleExpression>(f64, Float64Field, target) {
    override val instance by lazy {
        object : DoubleExpression {
            override val indexer = SimpleSymbolIndexer(keys)

-            override fun invoke(arguments: DoubleArray) = spreader(executable, arguments).unsafeCast<Double>()
+            override fun invoke(arguments: DoubleArray) = spreader(executable, arguments).unsafeCast<Float64>()
        }
    }

@ -99,7 +100,7 @@ internal class DoubleWasmBuilder(target: TypedMst<Double>) :

    override fun visitNumber(number: Number) = ctx.f64.const(number.toDouble())

-    override fun visitUnary(node: TypedMst.Unary<Double>): ExpressionRef = when (node.operation) {
+    override fun visitUnary(node: TypedMst.Unary<Float64>): ExpressionRef = when (node.operation) {
        GroupOps.MINUS_OPERATION -> ctx.f64.neg(visit(node.value))
        GroupOps.PLUS_OPERATION -> visit(node.value)
        PowerOperations.SQRT_OPERATION -> ctx.f64.sqrt(visit(node.value))
@ -120,7 +121,7 @@ internal class DoubleWasmBuilder(target: TypedMst<Double>) :
        else -> super.visitUnary(node)
    }

-    override fun visitBinary(mst: TypedMst.Binary<Double>): ExpressionRef = when (mst.operation) {
+    override fun visitBinary(mst: TypedMst.Binary<Float64>): ExpressionRef = when (mst.operation) {
        GroupOps.PLUS_OPERATION -> ctx.f64.add(visit(mst.left), visit(mst.right))
        GroupOps.MINUS_OPERATION -> ctx.f64.sub(visit(mst.left), visit(mst.right))
        RingOps.TIMES_OPERATION -> ctx.f64.mul(visit(mst.left), visit(mst.right))
--- a/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/wasm/wasm.kt
+++ b/kmath-ast/src/jsMain/kotlin/space/kscience/kmath/wasm/wasm.kt
@ -13,6 +13,7 @@ import space.kscience.kmath.ast.evaluateConstants
 import space.kscience.kmath.expressions.*
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.Int32Ring
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.wasm.internal.DoubleWasmBuilder
 import space.kscience.kmath.wasm.internal.IntWasmBuilder

@ -58,7 +59,7 @@ public fun MST.compile(algebra: Int32Ring, vararg arguments: Pair<Symbol, Int>):
 * @author Iaroslav Postovalov
 */
@UnstableKMathAPI
-public fun MST.compileToExpression(algebra: Float64Field): Expression<Double> {
+public fun MST.compileToExpression(algebra: Float64Field): Expression<Float64> {
    val typed = evaluateConstants(algebra)

    return if (typed is TypedMst.Constant) object : DoubleExpression {
--- a/kmath-ast/src/jsTest/kotlin/space/kscience/kmath/ast/utils.kt
+++ b/kmath-ast/src/jsTest/kotlin/space/kscience/kmath/ast/utils.kt
@ -13,6 +13,7 @@ import space.kscience.kmath.expressions.MST
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.Int32Ring
+import space.kscience.kmath.structures.Float64
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract
 import space.kscience.kmath.estree.compile as estreeCompile
@ -24,7 +25,7 @@ import space.kscience.kmath.wasm.compileToExpression as wasmCompileToExpression
 private object WasmCompilerTestContext : CompilerTestContext {
    override fun MST.compileToExpression(algebra: Int32Ring): Expression<Int> = wasmCompileToExpression(algebra)
    override fun MST.compile(algebra: Int32Ring, arguments: Map<Symbol, Int>): Int = wasmCompile(algebra, arguments)
-    override fun MST.compileToExpression(algebra: Float64Field): Expression<Double> = wasmCompileToExpression(algebra)
+    override fun MST.compileToExpression(algebra: Float64Field): Expression<Float64> = wasmCompileToExpression(algebra)

    override fun MST.compile(algebra: Float64Field, arguments: Map<Symbol, Double>): Double =
        wasmCompile(algebra, arguments)
@ -33,7 +34,7 @@ private object WasmCompilerTestContext : CompilerTestContext {
 private object ESTreeCompilerTestContext : CompilerTestContext {
    override fun MST.compileToExpression(algebra: Int32Ring): Expression<Int> = estreeCompileToExpression(algebra)
    override fun MST.compile(algebra: Int32Ring, arguments: Map<Symbol, Int>): Int = estreeCompile(algebra, arguments)
-    override fun MST.compileToExpression(algebra: Float64Field): Expression<Double> = estreeCompileToExpression(algebra)
+    override fun MST.compileToExpression(algebra: Float64Field): Expression<Float64> = estreeCompileToExpression(algebra)

    override fun MST.compile(algebra: Float64Field, arguments: Map<Symbol, Double>): Double =
        estreeCompile(algebra, arguments)
--- a/kmath-ast/src/jvmMain/kotlin/space/kscience/kmath/asm/internal/PrimitiveAsmBuilder.kt
+++ b/kmath-ast/src/jvmMain/kotlin/space/kscience/kmath/asm/internal/PrimitiveAsmBuilder.kt
@ -15,6 +15,7 @@ import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.ast.TypedMst
 import space.kscience.kmath.expressions.*
 import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.Float64
 import java.lang.invoke.MethodHandles
 import java.lang.invoke.MethodType
 import java.nio.file.Paths
@ -381,7 +382,7 @@ internal sealed class PrimitiveAsmBuilder<T : Number, out E : Expression<T>>(
 }

@UnstableKMathAPI
-internal class DoubleAsmBuilder(target: TypedMst<Double>) : PrimitiveAsmBuilder<Double, DoubleExpression>(
+internal class DoubleAsmBuilder(target: TypedMst<Float64>) : PrimitiveAsmBuilder<Double, DoubleExpression>(
    Float64Field,
    java.lang.Double::class.java,
    java.lang.Double.TYPE,
@ -410,7 +411,7 @@ internal class DoubleAsmBuilder(target: TypedMst<Double>) : PrimitiveAsmBuilder<
        false,
    )

-    override fun visitUnary(node: TypedMst.Unary<Double>) {
+    override fun visitUnary(node: TypedMst.Unary<Float64>) {
        super.visitUnary(node)

        when (node.operation) {
@ -435,7 +436,7 @@ internal class DoubleAsmBuilder(target: TypedMst<Double>) : PrimitiveAsmBuilder<
        }
    }

-    override fun visitBinary(node: TypedMst.Binary<Double>) {
+    override fun visitBinary(node: TypedMst.Binary<Float64>) {
        super.visitBinary(node)

        when (node.operation) {
--- a/kmath-ast/src/jvmTest/kotlin/space/kscience/kmath/ast/utils.kt
+++ b/kmath-ast/src/jvmTest/kotlin/space/kscience/kmath/ast/utils.kt
@ -12,8 +12,7 @@ import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.operations.Algebra
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.Int32Ring
-import kotlin.contracts.InvocationKind
-import kotlin.contracts.contract
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.asm.compile as asmCompile
 import space.kscience.kmath.asm.compileToExpression as asmCompileToExpression

@ -24,18 +23,18 @@ private object GenericAsmCompilerTestContext : CompilerTestContext {
    override fun MST.compile(algebra: Int32Ring, arguments: Map<Symbol, Int>): Int =
        asmCompile(algebra as Algebra<Int>, arguments)

-    override fun MST.compileToExpression(algebra: Float64Field): Expression<Double> =
-        asmCompileToExpression(algebra as Algebra<Double>)
+    override fun MST.compileToExpression(algebra: Float64Field): Expression<Float64> =
+        asmCompileToExpression(algebra as Algebra<Float64>)

    override fun MST.compile(algebra: Float64Field, arguments: Map<Symbol, Double>): Double =
-        asmCompile(algebra as Algebra<Double>, arguments)
+        asmCompile(algebra as Algebra<Float64>, arguments)
 }

@OptIn(UnstableKMathAPI::class)
 private object PrimitiveAsmCompilerTestContext : CompilerTestContext {
    override fun MST.compileToExpression(algebra: Int32Ring): Expression<Int> = asmCompileToExpression(algebra)
    override fun MST.compile(algebra: Int32Ring, arguments: Map<Symbol, Int>): Int = asmCompile(algebra, arguments)
-    override fun MST.compileToExpression(algebra: Float64Field): Expression<Double> = asmCompileToExpression(algebra)
+    override fun MST.compileToExpression(algebra: Float64Field): Expression<Float64> = asmCompileToExpression(algebra)

    override fun MST.compile(algebra: Float64Field, arguments: Map<Symbol, Double>): Double =
        asmCompile(algebra, arguments)
@ -43,7 +42,6 @@ private object PrimitiveAsmCompilerTestContext : CompilerTestContext {


 internal actual inline fun runCompilerTest(action: CompilerTestContext.() -> Unit) {
-    contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
    action(GenericAsmCompilerTestContext)
    action(PrimitiveAsmCompilerTestContext)
 }
--- a/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/expressions/CmDsExpression.kt
+++ b/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/expressions/CmDsExpression.kt
@ -14,6 +14,7 @@ import space.kscience.kmath.expressions.*
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.ExtendedField
 import space.kscience.kmath.operations.NumbersAddOps
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.MutableBufferFactory

 /**
@ -132,9 +133,9 @@ public object CmDsProcessor : AutoDiffProcessor<Double, DerivativeStructure, CmD
@Deprecated("Use generic DSAlgebra from the core")
 public class CmDsExpression(
    public val function: CmDsField.() -> DerivativeStructure,
-) : DifferentiableExpression<Double> {
+) : DifferentiableExpression<Float64> {

-    override val type: SafeType<Double> get() = DoubleField.type
+    override val type: SafeType<Float64> get() = DoubleField.type

    override operator fun invoke(arguments: Map<Symbol, Double>): Double =
        CmDsField(0, arguments).function().value
@ -142,7 +143,7 @@ public class CmDsExpression(
    /**
     * Get the derivative expression with given orders
     */
-    override fun derivativeOrNull(symbols: List<Symbol>): Expression<Double> = Expression(type) { arguments ->
+    override fun derivativeOrNull(symbols: List<Symbol>): Expression<Float64> = Expression(type) { arguments ->
        with(CmDsField(symbols.size, arguments)) { function().derivative(symbols) }
    }
 }
--- a/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/integration/CMGaussRuleIntegrator.kt
+++ b/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/integration/CMGaussRuleIntegrator.kt
@ -7,6 +7,7 @@ package space.kscience.kmath.commons.integration
 import org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator
 import org.apache.commons.math3.analysis.integration.gauss.GaussIntegratorFactory
 import space.kscience.kmath.integration.*
+import space.kscience.kmath.structures.Float64

 /**
 * A simple one-pass integrator based on Gauss rule
@ -14,9 +15,9 @@ import space.kscience.kmath.integration.*
 public class CMGaussRuleIntegrator(
    private val numpoints: Int,
    private var type: GaussRule = GaussRule.LEGENDRE,
-) : UnivariateIntegrator<Double> {
+) : UnivariateIntegrator<Float64> {

-    override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
+    override fun integrate(integrand: UnivariateIntegrand<Float64>): UnivariateIntegrand<Float64> {
        val range = integrand[IntegrationRange]
            ?: error("Integration range is not provided")
        val integrator: GaussIntegrator = getIntegrator(range)
@ -28,7 +29,7 @@ public class CMGaussRuleIntegrator(
        }
    }

-    private fun getIntegrator(range: ClosedRange<Double>): GaussIntegrator {
+    private fun getIntegrator(range: ClosedRange<Float64>): GaussIntegrator {
        return when (type) {
            GaussRule.LEGENDRE -> factory.legendre(
                numpoints,
@ -77,7 +78,7 @@ public class CMGaussRuleIntegrator(
        private val factory: GaussIntegratorFactory = GaussIntegratorFactory()

        public fun integrate(
-            range: ClosedRange<Double>,
+            range: ClosedRange<Float64>,
            numPoints: Int = 100,
            type: GaussRule = GaussRule.LEGENDRE,
            function: (Double) -> Double,
--- a/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt
+++ b/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt
@ -10,6 +10,7 @@ import org.apache.commons.math3.analysis.integration.SimpsonIntegrator
 import space.kscience.attributes.Attributes
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.integration.*
+import space.kscience.kmath.structures.Float64
 import org.apache.commons.math3.analysis.integration.UnivariateIntegrator as CMUnivariateIntegrator

 /**
@ -17,10 +18,10 @@ import org.apache.commons.math3.analysis.integration.UnivariateIntegrator as CMU
 */
 public class CMIntegrator(
    private val defaultMaxCalls: Int = 200,
-    public val integratorBuilder: (Integrand<Double>) -> CMUnivariateIntegrator,
-) : UnivariateIntegrator<Double> {
+    public val integratorBuilder: (Integrand<Float64>) -> CMUnivariateIntegrator,
+) : UnivariateIntegrator<Float64> {

-    override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
+    override fun integrate(integrand: UnivariateIntegrand<Float64>): UnivariateIntegrand<Float64> {
        val integrator = integratorBuilder(integrand)
        val maxCalls = integrand[IntegrandMaxCalls] ?: defaultMaxCalls
        val remainingCalls = maxCalls - integrand.calls
--- a/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt
+++ b/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt
@ -12,6 +12,7 @@ import space.kscience.attributes.SafeType
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.linear.*
 import space.kscience.kmath.linear.CholeskyDecomposition
+import space.kscience.kmath.linear.EigenDecomposition
 import space.kscience.kmath.linear.QRDecomposition
 import space.kscience.kmath.nd.Structure2D
 import space.kscience.kmath.nd.StructureAttribute
@ -22,7 +23,7 @@ import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.IntBuffer
 import space.kscience.kmath.structures.asBuffer

-public class CMMatrix(public val origin: RealMatrix) : Matrix<Double> {
+public class CMMatrix(public val origin: RealMatrix) : Matrix<Float64> {

    override val rowNum: Int get() = origin.rowDimension
    override val colNum: Int get() = origin.columnDimension
@ -31,12 +32,12 @@ public class CMMatrix(public val origin: RealMatrix) : Matrix<Double> {
 }

@JvmInline
-public value class CMVector(public val origin: RealVector) : Point<Double> {
+public value class CMVector(public val origin: RealVector) : Point<Float64> {
    override val size: Int get() = origin.dimension

    override operator fun get(index: Int): Double = origin.getEntry(index)

-    override operator fun iterator(): Iterator<Double> = origin.toArray().iterator()
+    override operator fun iterator(): Iterator<Float64> = origin.toArray().iterator()

    override fun toString(): String = Buffer.toString(this)
 }
@ -46,7 +47,7 @@ public fun RealVector.toPoint(): CMVector = CMVector(this)
 public object CMLinearSpace : LinearSpace<Double, Float64Field> {
    override val elementAlgebra: Float64Field get() = Float64Field

-    override val type: SafeType<Double> get() = DoubleField.type
+    override val type: SafeType<Float64> get() = DoubleField.type

    override fun buildMatrix(
        rows: Int,
@ -58,7 +59,7 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {
    }

    @OptIn(UnstableKMathAPI::class)
-    public fun Matrix<Double>.toCM(): CMMatrix = when (val matrix = origin) {
+    public fun Matrix<Float64>.toCM(): CMMatrix = when (val matrix = origin) {
        is CMMatrix -> matrix
        else -> {
            //TODO add feature analysis
@ -67,7 +68,7 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {
        }
    }

-    public fun Point<Double>.toCM(): CMVector = if (this is CMVector) this else {
+    public fun Point<Float64>.toCM(): CMVector = if (this is CMVector) this else {
        val array = DoubleArray(size) { this[it] }
        ArrayRealVector(array).wrap()
    }
@ -75,40 +76,41 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {
    internal fun RealMatrix.wrap(): CMMatrix = CMMatrix(this)
    internal fun RealVector.wrap(): CMVector = CMVector(this)

-    override fun buildVector(size: Int, initializer: Float64Field.(Int) -> Double): Point<Double> =
+    override fun buildVector(size: Int, initializer: Float64Field.(Int) -> Double): Point<Float64> =
        ArrayRealVector(DoubleArray(size) { Float64Field.initializer(it) }).wrap()

-    override fun Matrix<Double>.plus(other: Matrix<Double>): CMMatrix =
+    override fun Matrix<Float64>.plus(other: Matrix<Float64>): CMMatrix =
        toCM().origin.add(other.toCM().origin).wrap()

-    override fun Point<Double>.plus(other: Point<Double>): CMVector =
+    override fun Point<Float64>.plus(other: Point<Float64>): CMVector =
        toCM().origin.add(other.toCM().origin).wrap()

-    override fun Point<Double>.minus(other: Point<Double>): CMVector =
+    override fun Point<Float64>.minus(other: Point<Float64>): CMVector =
        toCM().origin.subtract(other.toCM().origin).wrap()

-    override fun Matrix<Double>.dot(other: Matrix<Double>): CMMatrix =
+    override fun Matrix<Float64>.dot(other: Matrix<Float64>): CMMatrix =
        toCM().origin.multiply(other.toCM().origin).wrap()

-    override fun Matrix<Double>.dot(vector: Point<Double>): CMVector =
+    override fun Matrix<Float64>.dot(vector: Point<Float64>): CMVector =
        toCM().origin.preMultiply(vector.toCM().origin).wrap()

-    override operator fun Matrix<Double>.minus(other: Matrix<Double>): CMMatrix =
+    override operator fun Matrix<Float64>.minus(other: Matrix<Float64>): CMMatrix =
        toCM().origin.subtract(other.toCM().origin).wrap()

-    override operator fun Matrix<Double>.times(value: Double): CMMatrix =
+    override operator fun Matrix<Float64>.times(value: Double): CMMatrix =
        toCM().origin.scalarMultiply(value).wrap()

-    override fun Double.times(m: Matrix<Double>): CMMatrix =
+    override fun Double.times(m: Matrix<Float64>): CMMatrix =
        m * this

-    override fun Point<Double>.times(value: Double): CMVector =
+    override fun Point<Float64>.times(value: Double): CMVector =
        toCM().origin.mapMultiply(value).wrap()

-    override fun Double.times(v: Point<Double>): CMVector =
+    override fun Double.times(v: Point<Float64>): CMVector =
        v * this

-    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Double>, attribute: A): V? {
+    @OptIn(UnstableKMathAPI::class)
+    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Float64>, attribute: A): V? {

        val origin = structure.toCM().origin

@ -125,7 +127,7 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {

            Cholesky -> object : CholeskyDecomposition<Float64> {
                val cmCholesky by lazy { org.apache.commons.math3.linear.CholeskyDecomposition(origin) }
-                override val l: Matrix<Double> get() = cmCholesky.l.wrap()
+                override val l: Matrix<Float64> get() = cmCholesky.l.wrap()
            }

            QR -> object : QRDecomposition<Float64> {
@ -144,6 +146,13 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {

            }

+            EIG -> object : EigenDecomposition<Float64> {
+                val cmEigen by lazy { org.apache.commons.math3.linear.EigenDecomposition(origin) }
+
+                override val v: Matrix<Float64> get() = cmEigen.v.wrap()
+                override val d: Matrix<Float64> get() = cmEigen.d.wrap()
+            }
+
            else -> null
        }
        @Suppress("UNCHECKED_CAST")
--- a/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/linear/CMSolver.kt
+++ b/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/linear/CMSolver.kt
@ -9,6 +9,7 @@ import org.apache.commons.math3.linear.*
 import space.kscience.kmath.linear.LinearSolver
 import space.kscience.kmath.linear.Matrix
 import space.kscience.kmath.linear.Point
+import space.kscience.kmath.structures.Float64

 public enum class CMDecomposition {
    LUP,
@ -19,7 +20,7 @@ public enum class CMDecomposition {
 }

 private fun CMLinearSpace.solver(
-    a: Matrix<Double>,
+    a: Matrix<Float64>,
    decomposition: CMDecomposition = CMDecomposition.LUP,
 ): DecompositionSolver = when (decomposition) {
    CMDecomposition.LUP -> LUDecomposition(a.toCM().origin).solver
@ -30,31 +31,31 @@ private fun CMLinearSpace.solver(
 }

 public fun CMLinearSpace.solve(
-    a: Matrix<Double>,
-    b: Matrix<Double>,
+    a: Matrix<Float64>,
+    b: Matrix<Float64>,
    decomposition: CMDecomposition = CMDecomposition.LUP,
 ): CMMatrix = solver(a, decomposition).solve(b.toCM().origin).wrap()

 public fun CMLinearSpace.solve(
-    a: Matrix<Double>,
-    b: Point<Double>,
+    a: Matrix<Float64>,
+    b: Point<Float64>,
    decomposition: CMDecomposition = CMDecomposition.LUP,
 ): CMVector = solver(a, decomposition).solve(b.toCM().origin).toPoint()

 public fun CMLinearSpace.inverse(
-    a: Matrix<Double>,
+    a: Matrix<Float64>,
    decomposition: CMDecomposition = CMDecomposition.LUP,
 ): CMMatrix = solver(a, decomposition).inverse.wrap()


-public fun CMLinearSpace.solver(decomposition: CMDecomposition): LinearSolver<Double> = object : LinearSolver<Double> {
-    override fun solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> =
+public fun CMLinearSpace.solver(decomposition: CMDecomposition): LinearSolver<Float64> = object : LinearSolver<Float64> {
+    override fun solve(a: Matrix<Float64>, b: Matrix<Float64>): Matrix<Float64> =
        solver(a, decomposition).solve(b.toCM().origin).wrap()

-    override fun solve(a: Matrix<Double>, b: Point<Double>): Point<Double> =
+    override fun solve(a: Matrix<Float64>, b: Point<Float64>): Point<Float64> =
        solver(a, decomposition).solve(b.toCM().origin).toPoint()

-    override fun inverse(matrix: Matrix<Double>): Matrix<Double> = solver(matrix, decomposition).inverse.wrap()
+    override fun inverse(matrix: Matrix<Float64>): Matrix<Float64> = solver(matrix, decomposition).inverse.wrap()
 }

-public fun CMLinearSpace.lupSolver(): LinearSolver<Double> = solver((CMDecomposition.LUP))
+public fun CMLinearSpace.lupSolver(): LinearSolver<Float64> = solver((CMDecomposition.LUP))
--- a/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/optimization/CMOptimizer.kt
+++ b/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/optimization/CMOptimizer.kt
@ -22,6 +22,7 @@ import space.kscience.kmath.expressions.SymbolIndexer
 import space.kscience.kmath.expressions.derivative
 import space.kscience.kmath.expressions.withSymbols
 import space.kscience.kmath.optimization.*
+import space.kscience.kmath.structures.Float64
 import kotlin.collections.set
 import kotlin.reflect.KClass

@ -33,7 +34,7 @@ public object CMOptimizerEngine : OptimizationAttribute<() -> MultivariateOptimi
 /**
 * Specify a Commons-maths optimization engine
 */
-public fun AttributesBuilder<FunctionOptimization<Double>>.cmEngine(optimizerBuilder: () -> MultivariateOptimizer) {
+public fun AttributesBuilder<FunctionOptimization<Float64>>.cmEngine(optimizerBuilder: () -> MultivariateOptimizer) {
    set(CMOptimizerEngine, optimizerBuilder)
 }

@ -42,18 +43,18 @@ public object CMOptimizerData : SetAttribute<SymbolIndexer.() -> OptimizationDat
 /**
 * Specify Commons-maths optimization data.
 */
-public fun AttributesBuilder<FunctionOptimization<Double>>.cmOptimizationData(data: SymbolIndexer.() -> OptimizationData) {
+public fun AttributesBuilder<FunctionOptimization<Float64>>.cmOptimizationData(data: SymbolIndexer.() -> OptimizationData) {
    CMOptimizerData add data
 }

-public fun AttributesBuilder<FunctionOptimization<Double>>.simplexSteps(vararg steps: Pair<Symbol, Double>) {
+public fun AttributesBuilder<FunctionOptimization<Float64>>.simplexSteps(vararg steps: Pair<Symbol, Double>) {
    //TODO use convergence checker from features
    cmEngine { SimplexOptimizer(CMOptimizer.defaultConvergenceChecker) }
    cmOptimizationData { NelderMeadSimplex(mapOf(*steps).toDoubleArray()) }
 }

@OptIn(UnstableKMathAPI::class)
-public object CMOptimizer : Optimizer<Double, FunctionOptimization<Double>> {
+public object CMOptimizer : Optimizer<Double, FunctionOptimization<Float64>> {

    public const val DEFAULT_RELATIVE_TOLERANCE: Double = 1e-4
    public const val DEFAULT_ABSOLUTE_TOLERANCE: Double = 1e-4
@ -67,12 +68,12 @@ public object CMOptimizer : Optimizer<Double, FunctionOptimization<Double>> {


    override suspend fun optimize(
-        problem: FunctionOptimization<Double>,
-    ): FunctionOptimization<Double> {
+        problem: FunctionOptimization<Float64>,
+    ): FunctionOptimization<Float64> {
        val startPoint = problem.startPoint

        val parameters = problem.attributes[OptimizationParameters]
-            ?: problem.attributes[OptimizationStartPoint<Double>()]?.keys
+            ?: problem.attributes[OptimizationStartPoint<Float64>()]?.keys
            ?: startPoint.keys


--- a/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/transform/Transformations.kt
+++ b/kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/transform/Transformations.kt
@ -77,7 +77,7 @@ public fun Flow<Buffer<Complex>>.fft(
 }

@JvmName("realFFT")
-public fun Flow<Buffer<Double>>.fft(
+public fun Flow<Buffer<Float64>>.fft(
    normalization: DftNormalization = DftNormalization.STANDARD,
    direction: TransformType = TransformType.FORWARD,
 ): Flow<Buffer<Complex>> {
@ -89,7 +89,7 @@ public fun Flow<Buffer<Double>>.fft(
 * Process a continuous flow of real numbers in FFT splitting it in chunks of [bufferSize].
 */
@JvmName("realFFT")
-public fun Flow<Double>.fft(
+public fun Flow<Float64>.fft(
    bufferSize: Int = Int.MAX_VALUE,
    normalization: DftNormalization = DftNormalization.STANDARD,
    direction: TransformType = TransformType.FORWARD,
@ -99,4 +99,4 @@ public fun Flow<Double>.fft(
 * Map a complex flow into real flow by taking real part of each number
 */
@FlowPreview
-public fun Flow<Complex>.real(): Flow<Double> = map { it.re }
+public fun Flow<Complex>.real(): Flow<Float64> = map { it.re }
--- a/kmath-commons/src/jvmTest/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt
+++ b/kmath-commons/src/jvmTest/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt
@ -18,6 +18,7 @@ import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.optimization.*
 import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.stat.chiSquaredExpression
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer
 import space.kscience.kmath.structures.asBuffer
 import kotlin.test.Test
@ -70,7 +71,7 @@ internal class OptimizeTest {
            bindSymbol(a) * arg.pow(2) + bindSymbol(b) * arg + cWithDefault
        }

-        val result: FunctionOptimization<Double> = chi2.optimizeWith(
+        val result: FunctionOptimization<Float64> = chi2.optimizeWith(
            CMOptimizer,
            mapOf(a to 1.5, b to 0.9, c to 1.0),
        ) {
--- a/kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/Quaternion.kt
+++ b/kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/Quaternion.kt
@ -11,6 +11,7 @@ import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.memory.*
 import space.kscience.kmath.operations.*
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.MutableBuffer
 import space.kscience.kmath.structures.MutableBufferFactory
 import kotlin.math.*
@ -28,7 +29,7 @@ public class Quaternion(
    public val x: Double,
    public val y: Double,
    public val z: Double,
-) : Buffer<Double> {
+) : Buffer<Float64> {
    init {
        require(!w.isNaN()) { "w-component of quaternion is not-a-number" }
        require(!x.isNaN()) { "x-component of quaternion is not-a-number" }
@ -51,7 +52,7 @@ public class Quaternion(
        else -> error("Index $index out of bounds [0,3]")
    }

-    override fun iterator(): Iterator<Double> = listOf(w, x, y, z).iterator()
+    override fun iterator(): Iterator<Float64> = listOf(w, x, y, z).iterator()

    override fun equals(other: Any?): Boolean {
        if (this === other) return true
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/Domain1D.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/Domain1D.kt
@ -7,6 +7,7 @@ package space.kscience.kmath.domains

 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.linear.Point
+import space.kscience.kmath.structures.Float64

@UnstableKMathAPI
 public abstract class Domain1D<T : Comparable<T>>(public val range: ClosedRange<T>) : Domain<T> {
@ -22,8 +23,8 @@ public abstract class Domain1D<T : Comparable<T>>(public val range: ClosedRange<

@UnstableKMathAPI
 public class DoubleDomain1D(
-    @Suppress("CanBeParameter") public val doubleRange: ClosedFloatingPointRange<Double>,
-) : Domain1D<Double>(doubleRange), Float64Domain {
+    @Suppress("CanBeParameter") public val doubleRange: ClosedFloatingPointRange<Float64>,
+) : Domain1D<Float64>(doubleRange), Float64Domain {
    override fun getLowerBound(num: Int): Double {
        require(num == 0)
        return range.start
@ -55,5 +56,5 @@ public class DoubleDomain1D(
 }

@UnstableKMathAPI
-public val Domain1D<Double>.center: Double
+public val Domain1D<Float64>.center: Double
    get() = (range.endInclusive + range.start) / 2
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/Float64Domain.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/Float64Domain.kt
@ -5,6 +5,7 @@
 package space.kscience.kmath.domains

 import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.structures.Float64

 /**
 * n-dimensional volume
@ -12,7 +13,7 @@ import space.kscience.kmath.UnstableKMathAPI
 * @author Alexander Nozik
 */
@UnstableKMathAPI
-public interface Float64Domain : Domain<Double> {
+public interface Float64Domain : Domain<Float64> {
    /**
     * Global lower edge
     * @param num axis number
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/HyperSquareDomain.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/HyperSquareDomain.kt
@ -7,6 +7,7 @@ package space.kscience.kmath.domains
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.linear.Point
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer
 import space.kscience.kmath.structures.indices

@ -15,7 +16,7 @@ import space.kscience.kmath.structures.indices
 * and a [Buffer] of upper boundaries. Upper should be greater or equals than lower.
 */
@UnstableKMathAPI
-public class HyperSquareDomain(public val lower: Buffer<Double>, public val upper: Buffer<Double>) : Float64Domain {
+public class HyperSquareDomain(public val lower: Buffer<Float64>, public val upper: Buffer<Float64>) : Float64Domain {
    init {
        require(lower.size == upper.size) {
            "Domain borders size mismatch. Lower borders size is ${lower.size}, but upper borders size is ${upper.size}."
@ -29,7 +30,7 @@ public class HyperSquareDomain(public val lower: Buffer<Double>, public val uppe

    public val center: Float64Buffer get() = Float64Buffer(dimension) { (lower[it] + upper[it]) / 2.0 }

-    override operator fun contains(point: Point<Double>): Boolean = point.indices.all { i ->
+    override operator fun contains(point: Point<Float64>): Boolean = point.indices.all { i ->
        point[i] in lower[i]..upper[i]
    }

--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/UnconstrainedDomain.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/UnconstrainedDomain.kt
@ -6,10 +6,11 @@ package space.kscience.kmath.domains

 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.linear.Point
+import space.kscience.kmath.structures.Float64

@UnstableKMathAPI
 public class UnconstrainedDomain(override val dimension: Int) : Float64Domain {
-    override operator fun contains(point: Point<Double>): Boolean = true
+    override operator fun contains(point: Point<Float64>): Boolean = true

    override fun getLowerBound(num: Int): Double = Double.NEGATIVE_INFINITY

--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/Expression.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/Expression.kt
@ -13,6 +13,7 @@ import space.kscience.kmath.operations.Algebra
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.IntRing
 import space.kscience.kmath.operations.LongRing
+import space.kscience.kmath.structures.Float64
 import kotlin.jvm.JvmName
 import kotlin.properties.ReadOnlyProperty

@ -47,9 +48,9 @@ public inline fun <reified T> Expression(noinline block: (Map<Symbol, T>) -> T):
 * Specialization of [Expression] for [Double] allowing better performance because of using array.
 */
@UnstableKMathAPI
-public interface DoubleExpression : Expression<Double> {
+public interface DoubleExpression : Expression<Float64> {

-    override val type: SafeType<Double> get() = DoubleField.type
+    override val type: SafeType<Float64> get() = DoubleField.type

    /**
     * The indexer of this expression's arguments that should be used to build array for [invoke].
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/FunctionalExpressionAlgebra.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/FunctionalExpressionAlgebra.kt
@ -6,6 +6,7 @@
 package space.kscience.kmath.expressions

 import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.MutableBufferFactory
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract
@ -197,5 +198,5 @@ public inline fun <T, A : ExtendedField<T>> A.expressionInExtendedField(
 ): Expression<T> = FunctionalExpressionExtendedField(this).block()

 public inline fun Float64Field.expression(
-    block: FunctionalExpressionExtendedField<Double, Float64Field>.() -> Expression<Double>,
-): Expression<Double> = FunctionalExpressionExtendedField(this).block()
+    block: FunctionalExpressionExtendedField<Double, Float64Field>.() -> Expression<Float64>,
+): Expression<Float64> = FunctionalExpressionExtendedField(this).block()
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/EigenDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/EigenDecomposition.kt
@ -0,0 +1,32 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.linear
+
+import space.kscience.attributes.PolymorphicAttribute
+import space.kscience.attributes.safeTypeOf
+import space.kscience.kmath.UnstableKMathAPI
+
+@UnstableKMathAPI
+public interface EigenDecomposition<T> {
+    /**
+     * Eigenvector matrix.
+     */
+    public val v: Matrix<T>
+
+    /**
+     * A diagonal matrix of eigenvalues. Must have [IsDiagonal]
+     */
+    public val d: Matrix<T>
+}
+
+@UnstableKMathAPI
+public class EigenDecompositionAttribute<T> :
+    PolymorphicAttribute<EigenDecomposition<T>>(safeTypeOf()),
+    MatrixAttribute<EigenDecomposition<T>>
+
+@UnstableKMathAPI
+public val <T> MatrixScope<T>.EIG: EigenDecompositionAttribute<T>
+    get() = EigenDecompositionAttribute()
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/Float64LinearSpace.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/Float64LinearSpace.kt
@ -6,12 +6,19 @@
 package space.kscience.kmath.linear

 import space.kscience.kmath.PerformancePitfall
-import space.kscience.kmath.nd.*
+import space.kscience.kmath.nd.Floa64FieldOpsND
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.asND
 import space.kscience.kmath.operations.Float64BufferOps
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer
+import kotlin.collections.component1
+import kotlin.collections.component2
+import kotlin.collections.map

 public object Float64LinearSpace : LinearSpace<Double, Float64Field> {

@ -21,36 +28,36 @@ public object Float64LinearSpace : LinearSpace<Double, Float64Field> {
        rows: Int,
        columns: Int,
        initializer: Float64Field.(i: Int, j: Int) -> Double,
-    ): Matrix<Double> = Floa64FieldOpsND.structureND(ShapeND(rows, columns)) { (i, j) ->
+    ): Matrix<Float64> = Floa64FieldOpsND.structureND(ShapeND(rows, columns)) { (i, j) ->
        Float64Field.initializer(i, j)
    }.as2D()

    override fun buildVector(size: Int, initializer: Float64Field.(Int) -> Double): Float64Buffer =
        Float64Buffer(size) { Float64Field.initializer(it) }

-    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = Floa64FieldOpsND {
+    override fun Matrix<Float64>.unaryMinus(): Matrix<Float64> = Floa64FieldOpsND {
        asND().map { -it }.as2D()
    }

-    override fun Matrix<Double>.plus(other: Matrix<Double>): Matrix<Double> = Floa64FieldOpsND {
+    override fun Matrix<Float64>.plus(other: Matrix<Float64>): Matrix<Float64> = Floa64FieldOpsND {
        require(shape == other.shape) { "Shape mismatch on Matrix::plus. Expected $shape but found ${other.shape}" }
        asND().plus(other.asND()).as2D()
    }

-    override fun Matrix<Double>.minus(other: Matrix<Double>): Matrix<Double> = Floa64FieldOpsND {
+    override fun Matrix<Float64>.minus(other: Matrix<Float64>): Matrix<Float64> = Floa64FieldOpsND {
        require(shape == other.shape) { "Shape mismatch on Matrix::minus. Expected $shape but found ${other.shape}" }
        asND().minus(other.asND()).as2D()
    }

    // Create a continuous in-memory representation of this vector for better memory layout handling
-    private fun Buffer<Double>.linearize() = if (this is Float64Buffer) {
+    private fun Buffer<Float64>.linearize() = if (this is Float64Buffer) {
        this.array
    } else {
        DoubleArray(size) { get(it) }
    }

    @OptIn(PerformancePitfall::class)
-    override fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> {
+    override fun Matrix<Float64>.dot(other: Matrix<Float64>): Matrix<Float64> {
        require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
        val rows = this@dot.rows.map { it.linearize() }
        val columns = other.columns.map { it.linearize() }
@ -67,7 +74,7 @@ public object Float64LinearSpace : LinearSpace<Double, Float64Field> {
    }

    @OptIn(PerformancePitfall::class)
-    override fun Matrix<Double>.dot(vector: Point<Double>): Float64Buffer {
+    override fun Matrix<Float64>.dot(vector: Point<Float64>): Float64Buffer {
        require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
        val rows = this@dot.rows.map { it.linearize() }
        val indices = 0 until this.colNum
@ -82,27 +89,27 @@ public object Float64LinearSpace : LinearSpace<Double, Float64Field> {

    }

-    override fun Matrix<Double>.times(value: Double): Matrix<Double> = Floa64FieldOpsND {
+    override fun Matrix<Float64>.times(value: Double): Matrix<Float64> = Floa64FieldOpsND {
        asND().map { it * value }.as2D()
    }

-    public override fun Point<Double>.plus(other: Point<Double>): Float64Buffer = Float64BufferOps.run {
+    public override fun Point<Float64>.plus(other: Point<Float64>): Float64Buffer = Float64BufferOps.run {
        this@plus + other
    }

-    public override fun Point<Double>.minus(other: Point<Double>): Float64Buffer = Float64BufferOps.run {
+    public override fun Point<Float64>.minus(other: Point<Float64>): Float64Buffer = Float64BufferOps.run {
        this@minus - other
    }

-    public override fun Point<Double>.times(value: Double): Float64Buffer = Float64BufferOps.run {
+    public override fun Point<Float64>.times(value: Double): Float64Buffer = Float64BufferOps.run {
        scale(this@times, value)
    }

-    public operator fun Point<Double>.div(value: Double): Float64Buffer = Float64BufferOps.run {
+    public operator fun Point<Float64>.div(value: Double): Float64Buffer = Float64BufferOps.run {
        scale(this@div, 1.0 / value)
    }

-    public override fun Double.times(v: Point<Double>): Float64Buffer = v * this
+    public override fun Double.times(v: Point<Float64>): Float64Buffer = v * this


 }
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt
@ -163,9 +163,9 @@ public fun <T : Comparable<T>> Field<T>.lup(


 public fun Field<Float64>.lup(
-    matrix: Matrix<Double>,
+    matrix: Matrix<Float64>,
    singularityThreshold: Double = 1e-11,
-): GenericLupDecomposition<Double> = lup(matrix) { it < singularityThreshold }
+): GenericLupDecomposition<Float64> = lup(matrix) { it < singularityThreshold }

 private fun <T> Field<T>.solve(
    lup: LupDecomposition<T>,
@ -235,5 +235,5 @@ public fun <T : Comparable<T>> LinearSpace<T, Field<T>>.lupSolver(
    override fun inverse(matrix: Matrix<T>): Matrix<T> = solve(matrix, one(matrix.rowNum, matrix.colNum))
 }

-public fun LinearSpace<Double, Float64Field>.lupSolver(singularityThreshold: Double = 1e-11): LinearSolver<Double> =
+public fun LinearSpace<Double, Float64Field>.lupSolver(singularityThreshold: Double = 1e-11): LinearSolver<Float64> =
    lupSolver { it < singularityThreshold }
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt
@ -67,13 +67,14 @@ public fun <T : Any, A : Ring<T>> MatrixBuilder<T, A>.symmetric(
 ): Matrix<T> {
    require(columns == rows) { "In order to build symmetric matrix, number of rows $rows should be equal to number of columns $columns" }
    return with(BufferAccessor2D<T?>(rows, rows, MutableBufferFactory(type))) {
-        val cache = factory(rows * rows) { null }
+        val cache = HashMap<IntArray, T>()
        linearSpace.buildMatrix(rows, rows) { i, j ->
-            val cached = cache[i, j]
+            val index = intArrayOf(i, j)
+            val cached = cache[index]
            if (cached == null) {
                val value = if (i <= j) builder(i, j) else builder(j, i)
-                cache[i, j] = value
-                cache[j, i] = value
+                cache[index] = value
+                cache[index] = value
                value
            } else {
                cached
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixAttributes.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixAttributes.kt
@ -124,7 +124,7 @@ public val <T> MatrixScope<T>.QR: QRDecompositionAttribute<T>

 public interface CholeskyDecomposition<T> {
    /**
-     * The triangular matrix in this decomposition. It may have either [UpperTriangular] or [LowerTriangular].
+     * The lower triangular matrix in this decomposition. It should have [LowerTriangular].
     */
    public val l: Matrix<T>
 }
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/cumulative.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/cumulative.kt
@ -8,6 +8,7 @@ package space.kscience.kmath.misc
 import space.kscience.kmath.operations.Ring
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import kotlin.jvm.JvmName

 /**
@ -47,7 +48,7 @@ public fun <T> Iterable<T>.cumulativeSum(ring: Ring<T>): Iterable<T> =
    ring { cumulative(zero) { element: T, sum: T -> sum + element } }

@JvmName("cumulativeSumOfDouble")
-public fun Iterable<Double>.cumulativeSum(): Iterable<Double> = cumulative(0.0) { element, sum -> sum + element }
+public fun Iterable<Float64>.cumulativeSum(): Iterable<Float64> = cumulative(0.0) { element, sum -> sum + element }

@JvmName("cumulativeSumOfInt")
 public fun Iterable<Int>.cumulativeSum(): Iterable<Int> = cumulative(0) { element, sum -> sum + element }
@ -59,7 +60,7 @@ public fun <T> Sequence<T>.cumulativeSum(ring: Ring<T>): Sequence<T> =
    ring { cumulative(zero) { element: T, sum: T -> sum + element } }

@JvmName("cumulativeSumOfDouble")
-public fun Sequence<Double>.cumulativeSum(): Sequence<Double> = cumulative(0.0) { element, sum -> sum + element }
+public fun Sequence<Float64>.cumulativeSum(): Sequence<Float64> = cumulative(0.0) { element, sum -> sum + element }

@JvmName("cumulativeSumOfInt")
 public fun Sequence<Int>.cumulativeSum(): Sequence<Int> = cumulative(0) { element, sum -> sum + element }
@ -71,7 +72,7 @@ public fun <T> List<T>.cumulativeSum(group: Ring<T>): List<T> =
    group { cumulative(zero) { element: T, sum: T -> sum + element } }

@JvmName("cumulativeSumOfDouble")
-public fun List<Double>.cumulativeSum(): List<Double> = cumulative(0.0) { element, sum -> sum + element }
+public fun List<Float64>.cumulativeSum(): List<Float64> = cumulative(0.0) { element, sum -> sum + element }

@JvmName("cumulativeSumOfInt")
 public fun List<Int>.cumulativeSum(): List<Int> = cumulative(0) { element, sum -> sum + element }
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Float64FieldND.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Float64FieldND.kt
@ -8,6 +8,7 @@ package space.kscience.kmath.nd
 import space.kscience.kmath.PerformancePitfall
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract
@ -20,7 +21,7 @@ import kotlin.math.pow as kpow
 public class Float64BufferND(
    indexes: ShapeIndexer,
    override val buffer: Float64Buffer,
-) : MutableBufferND<Double>(indexes, buffer), MutableStructureNDOfDouble {
+) : MutableBufferND<Float64>(indexes, buffer), MutableStructureNDOfDouble {

    override fun getDouble(index: IntArray): Double = buffer[indices.offset(index)]

@ -30,11 +31,11 @@ public class Float64BufferND(
 }


-public sealed class Floa64FieldOpsND : BufferedFieldOpsND<Double, Float64Field>(Float64Field.bufferAlgebra),
-    ScaleOperations<StructureND<Double>>, ExtendedFieldOps<StructureND<Double>> {
+public sealed class Floa64FieldOpsND : BufferedFieldOpsND<Float64, Float64Field>(Float64Field.bufferAlgebra),
+    ScaleOperations<StructureND<Float64>>, ExtendedFieldOps<StructureND<Float64>> {

    @OptIn(PerformancePitfall::class)
-    override fun StructureND<Double>.toBufferND(): Float64BufferND = when (this) {
+    override fun StructureND<Float64>.toBufferND(): Float64BufferND = when (this) {
        is Float64BufferND -> this
        else -> {
            val indexer = indexerBuilder(shape)
@ -64,16 +65,16 @@ public sealed class Floa64FieldOpsND : BufferedFieldOpsND<Double, Float64Field>(
    }

    @OptIn(PerformancePitfall::class)
-    override fun StructureND<Double>.map(transform: Float64Field.(Double) -> Double): BufferND<Double> =
+    override fun StructureND<Float64>.map(transform: Float64Field.(Double) -> Double): BufferND<Float64> =
        mapInline(toBufferND()) { Float64Field.transform(it) }


    @OptIn(PerformancePitfall::class)
    override fun zip(
-        left: StructureND<Double>,
-        right: StructureND<Double>,
+        left: StructureND<Float64>,
+        right: StructureND<Float64>,
        transform: Float64Field.(Double, Double) -> Double,
-    ): BufferND<Double> = zipInline(left.toBufferND(), right.toBufferND()) { l, r -> Float64Field.transform(l, r) }
+    ): BufferND<Float64> = zipInline(left.toBufferND(), right.toBufferND()) { l, r -> Float64Field.transform(l, r) }

    override fun mutableStructureND(shape: ShapeND, initializer: Float64Field.(IntArray) -> Double): Float64BufferND {
        val indexer = indexerBuilder(shape)
@ -85,102 +86,102 @@ public sealed class Floa64FieldOpsND : BufferedFieldOpsND<Double, Float64Field>(
        )
    }

-    override fun add(left: StructureND<Double>, right: StructureND<Double>): Float64BufferND =
+    override fun add(left: StructureND<Float64>, right: StructureND<Float64>): Float64BufferND =
        zipInline(left.toBufferND(), right.toBufferND()) { l, r -> l + r }

-    override fun multiply(left: StructureND<Double>, right: StructureND<Double>): Float64BufferND =
+    override fun multiply(left: StructureND<Float64>, right: StructureND<Float64>): Float64BufferND =
        zipInline(left.toBufferND(), right.toBufferND()) { l, r -> l * r }

-    override fun StructureND<Double>.unaryMinus(): Float64BufferND = mapInline(toBufferND()) { -it }
+    override fun StructureND<Float64>.unaryMinus(): Float64BufferND = mapInline(toBufferND()) { -it }

-    override fun StructureND<Double>.div(arg: StructureND<Double>): Float64BufferND =
+    override fun StructureND<Float64>.div(arg: StructureND<Float64>): Float64BufferND =
        zipInline(toBufferND(), arg.toBufferND()) { l, r -> l / r }

-    override fun divide(left: StructureND<Double>, right: StructureND<Double>): Float64BufferND =
+    override fun divide(left: StructureND<Float64>, right: StructureND<Float64>): Float64BufferND =
        zipInline(left.toBufferND(), right.toBufferND()) { l: Double, r: Double -> l / r }

-    override fun StructureND<Double>.div(arg: Double): Float64BufferND =
+    override fun StructureND<Float64>.div(arg: Double): Float64BufferND =
        mapInline(toBufferND()) { it / arg }

-    override fun Double.div(arg: StructureND<Double>): Float64BufferND =
+    override fun Double.div(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { this / it }

-    override fun StructureND<Double>.unaryPlus(): Float64BufferND = toBufferND()
+    override fun StructureND<Float64>.unaryPlus(): Float64BufferND = toBufferND()

-    override fun StructureND<Double>.plus(arg: StructureND<Double>): Float64BufferND =
+    override fun StructureND<Float64>.plus(arg: StructureND<Float64>): Float64BufferND =
        zipInline(toBufferND(), arg.toBufferND()) { l: Double, r: Double -> l + r }

-    override fun StructureND<Double>.minus(arg: StructureND<Double>): Float64BufferND =
+    override fun StructureND<Float64>.minus(arg: StructureND<Float64>): Float64BufferND =
        zipInline(toBufferND(), arg.toBufferND()) { l: Double, r: Double -> l - r }

-    override fun StructureND<Double>.times(arg: StructureND<Double>): Float64BufferND =
+    override fun StructureND<Float64>.times(arg: StructureND<Float64>): Float64BufferND =
        zipInline(toBufferND(), arg.toBufferND()) { l: Double, r: Double -> l * r }

-    override fun StructureND<Double>.times(k: Number): Float64BufferND =
+    override fun StructureND<Float64>.times(k: Number): Float64BufferND =
        mapInline(toBufferND()) { it * k.toDouble() }

-    override fun StructureND<Double>.div(k: Number): Float64BufferND =
+    override fun StructureND<Float64>.div(k: Number): Float64BufferND =
        mapInline(toBufferND()) { it / k.toDouble() }

-    override fun Number.times(arg: StructureND<Double>): Float64BufferND = arg * this
+    override fun Number.times(arg: StructureND<Float64>): Float64BufferND = arg * this

-    override fun StructureND<Double>.plus(arg: Double): Float64BufferND = mapInline(toBufferND()) { it + arg }
+    override fun StructureND<Float64>.plus(arg: Double): Float64BufferND = mapInline(toBufferND()) { it + arg }

-    override fun StructureND<Double>.minus(arg: Double): StructureND<Double> = mapInline(toBufferND()) { it - arg }
+    override fun StructureND<Float64>.minus(arg: Double): StructureND<Float64> = mapInline(toBufferND()) { it - arg }

-    override fun Double.plus(arg: StructureND<Double>): StructureND<Double> = arg + this
+    override fun Double.plus(arg: StructureND<Float64>): StructureND<Float64> = arg + this

-    override fun Double.minus(arg: StructureND<Double>): StructureND<Double> = mapInline(arg.toBufferND()) { this - it }
+    override fun Double.minus(arg: StructureND<Float64>): StructureND<Float64> = mapInline(arg.toBufferND()) { this - it }

-    override fun scale(a: StructureND<Double>, value: Double): Float64BufferND =
+    override fun scale(a: StructureND<Float64>, value: Double): Float64BufferND =
        mapInline(a.toBufferND()) { it * value }

-    override fun exp(arg: StructureND<Double>): Float64BufferND =
+    override fun exp(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.exp(it) }

-    override fun ln(arg: StructureND<Double>): Float64BufferND =
+    override fun ln(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.ln(it) }

-    override fun sin(arg: StructureND<Double>): Float64BufferND =
+    override fun sin(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.sin(it) }

-    override fun cos(arg: StructureND<Double>): Float64BufferND =
+    override fun cos(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.cos(it) }

-    override fun tan(arg: StructureND<Double>): Float64BufferND =
+    override fun tan(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.tan(it) }

-    override fun asin(arg: StructureND<Double>): Float64BufferND =
+    override fun asin(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.asin(it) }

-    override fun acos(arg: StructureND<Double>): Float64BufferND =
+    override fun acos(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.acos(it) }

-    override fun atan(arg: StructureND<Double>): Float64BufferND =
+    override fun atan(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.atan(it) }

-    override fun sinh(arg: StructureND<Double>): Float64BufferND =
+    override fun sinh(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.sinh(it) }

-    override fun cosh(arg: StructureND<Double>): Float64BufferND =
+    override fun cosh(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.cosh(it) }

-    override fun tanh(arg: StructureND<Double>): Float64BufferND =
+    override fun tanh(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.tanh(it) }

-    override fun asinh(arg: StructureND<Double>): Float64BufferND =
+    override fun asinh(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.asinh(it) }

-    override fun acosh(arg: StructureND<Double>): Float64BufferND =
+    override fun acosh(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.acosh(it) }

-    override fun atanh(arg: StructureND<Double>): Float64BufferND =
+    override fun atanh(arg: StructureND<Float64>): Float64BufferND =
        mapInline(arg.toBufferND()) { kotlin.math.atanh(it) }

    override fun power(
-        arg: StructureND<Double>,
+        arg: StructureND<Float64>,
        pow: Number,
-    ): StructureND<Double> = if (pow is Int) {
+    ): StructureND<Float64> = if (pow is Int) {
        mapInline(arg.toBufferND()) { it.pow(pow) }
    } else {
        mapInline(arg.toBufferND()) { it.pow(pow.toDouble()) }
@ -191,18 +192,18 @@ public sealed class Floa64FieldOpsND : BufferedFieldOpsND<Double, Float64Field>(

@OptIn(UnstableKMathAPI::class)
 public class Float64FieldND(override val shape: ShapeND) :
-    Floa64FieldOpsND(), FieldND<Double, Float64Field>, NumbersAddOps<StructureND<Double>>,
-    ExtendedField<StructureND<Double>> {
+    Floa64FieldOpsND(), FieldND<Double, Float64Field>, NumbersAddOps<StructureND<Float64>>,
+    ExtendedField<StructureND<Float64>> {

-    override fun power(arg: StructureND<Double>, pow: UInt): Float64BufferND = mapInline(arg.toBufferND()) {
+    override fun power(arg: StructureND<Float64>, pow: UInt): Float64BufferND = mapInline(arg.toBufferND()) {
        it.kpow(pow.toInt())
    }

-    override fun power(arg: StructureND<Double>, pow: Int): Float64BufferND = mapInline(arg.toBufferND()) {
+    override fun power(arg: StructureND<Float64>, pow: Int): Float64BufferND = mapInline(arg.toBufferND()) {
        it.kpow(pow)
    }

-    override fun power(arg: StructureND<Double>, pow: Number): Float64BufferND = if (pow.isInteger()) {
+    override fun power(arg: StructureND<Float64>, pow: Number): Float64BufferND = if (pow.isInteger()) {
        power(arg, pow.toInt())
    } else {
        val dpow = pow.toDouble()
@ -212,17 +213,17 @@ public class Float64FieldND(override val shape: ShapeND) :
        }
    }

-    override fun sinh(arg: StructureND<Double>): Float64BufferND = super<Floa64FieldOpsND>.sinh(arg)
+    override fun sinh(arg: StructureND<Float64>): Float64BufferND = super<Floa64FieldOpsND>.sinh(arg)

-    override fun cosh(arg: StructureND<Double>): Float64BufferND = super<Floa64FieldOpsND>.cosh(arg)
+    override fun cosh(arg: StructureND<Float64>): Float64BufferND = super<Floa64FieldOpsND>.cosh(arg)

-    override fun tanh(arg: StructureND<Double>): Float64BufferND = super<Floa64FieldOpsND>.tan(arg)
+    override fun tanh(arg: StructureND<Float64>): Float64BufferND = super<Floa64FieldOpsND>.tan(arg)

-    override fun asinh(arg: StructureND<Double>): Float64BufferND = super<Floa64FieldOpsND>.asinh(arg)
+    override fun asinh(arg: StructureND<Float64>): Float64BufferND = super<Floa64FieldOpsND>.asinh(arg)

-    override fun acosh(arg: StructureND<Double>): Float64BufferND = super<Floa64FieldOpsND>.acosh(arg)
+    override fun acosh(arg: StructureND<Float64>): Float64BufferND = super<Floa64FieldOpsND>.acosh(arg)

-    override fun atanh(arg: StructureND<Double>): Float64BufferND = super<Floa64FieldOpsND>.atanh(arg)
+    override fun atanh(arg: StructureND<Float64>): Float64BufferND = super<Floa64FieldOpsND>.atanh(arg)

    override fun number(value: Number): Float64BufferND {
        val d = value.toDouble() // minimize conversions
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/ShapeND.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/ShapeND.kt
@ -6,7 +6,6 @@
 package space.kscience.kmath.nd

 import space.kscience.kmath.UnsafeKMathAPI
-import kotlin.jvm.JvmInline

 /**
 * A read-only ND shape
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt
@ -14,6 +14,7 @@ import space.kscience.kmath.linear.LinearSpace
 import space.kscience.kmath.operations.Ring
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import kotlin.math.abs

 public interface StructureAttribute<T> : Attribute<T>
@ -76,8 +77,8 @@ public interface StructureND<out T> : AttributeContainer, WithShape {

        @PerformancePitfall
        public fun contentEquals(
-            st1: StructureND<Double>,
-            st2: StructureND<Double>,
+            st1: StructureND<Float64>,
+            st2: StructureND<Float64>,
            tolerance: Double = 1e-11,
        ): Boolean {
            if (st1 === st2) return true
@ -210,7 +211,7 @@ public fun <T : Comparable<T>> LinearSpace<T, Ring<T>>.contentEquals(
@PerformancePitfall
 public operator fun <T> StructureND<T>.get(vararg index: Int): T = get(index)

-public operator fun StructureND<Double>.get(vararg index: Int): Double = getDouble(index)
+public operator fun StructureND<Float64>.get(vararg index: Int): Double = getDouble(index)

 public operator fun StructureND<Int>.get(vararg index: Int): Int = getInt(index)

--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/VirtualStructureND.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/VirtualStructureND.kt
@ -7,6 +7,7 @@ package space.kscience.kmath.nd

 import space.kscience.kmath.PerformancePitfall
 import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.structures.Float64

 public open class VirtualStructureND<T>(
    override val shape: ShapeND,
@ -24,7 +25,7 @@ public open class VirtualStructureND<T>(
 public class VirtualDoubleStructureND(
    shape: ShapeND,
    producer: (IntArray) -> Double,
-) : VirtualStructureND<Double>(shape, producer)
+) : VirtualStructureND<Float64>(shape, producer)

@UnstableKMathAPI
 public class VirtualIntStructureND(
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/primitiveStructureND.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/primitiveStructureND.kt
@ -6,8 +6,9 @@
 package space.kscience.kmath.nd

 import space.kscience.kmath.PerformancePitfall
+import space.kscience.kmath.structures.Float64

-public interface StructureNDOfDouble : StructureND<Double> {
+public interface StructureNDOfDouble : StructureND<Float64> {

    /**
     * Guaranteed non-blocking access to content
@ -19,10 +20,10 @@ public interface StructureNDOfDouble : StructureND<Double> {
 * Optimized method to access primitive without boxing if possible
 */
@OptIn(PerformancePitfall::class)
-public fun StructureND<Double>.getDouble(index: IntArray): Double =
+public fun StructureND<Float64>.getDouble(index: IntArray): Double =
    if (this is StructureNDOfDouble) getDouble(index) else get(index)

-public interface MutableStructureNDOfDouble : StructureNDOfDouble, MutableStructureND<Double> {
+public interface MutableStructureNDOfDouble : StructureNDOfDouble, MutableStructureND<Float64> {

    /**
     * Guaranteed non-blocking access to content
@ -31,7 +32,7 @@ public interface MutableStructureNDOfDouble : StructureNDOfDouble, MutableStruct
 }

@OptIn(PerformancePitfall::class)
-public fun MutableStructureND<Double>.getDouble(index: IntArray): Double =
+public fun MutableStructureND<Float64>.getDouble(index: IntArray): Double =
    if (this is StructureNDOfDouble) getDouble(index) else get(index)


--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt
@ -6,6 +6,7 @@
 package space.kscience.kmath.operations

 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.MutableBuffer
 import space.kscience.kmath.structures.MutableBufferFactory

@ -182,7 +183,7 @@ public fun <T, A : Field<T>> BufferFieldOps<T, A>.withSize(size: Int): BufferFie

 //Double buffer specialization

-public fun BufferField<Double, *>.buffer(vararg elements: Number): Buffer<Double> {
+public fun BufferField<Double, *>.buffer(vararg elements: Number): Buffer<Float64> {
    require(elements.size == size) { "Expected $size elements but found ${elements.size}" }
    return elementBufferFactory(size) { elements[it].toDouble() }
 }
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Float64BufferField.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Float64BufferField.kt
@ -6,6 +6,7 @@
 package space.kscience.kmath.operations

 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer

 /**
@ -13,23 +14,23 @@ import space.kscience.kmath.structures.Float64Buffer
 *
 * @property size the size of buffers to operate on.
 */
-public class Float64BufferField(public val size: Int) : ExtendedField<Buffer<Double>>, Float64BufferOps() {
-    override val zero: Buffer<Double> by lazy { Float64Buffer(size) { 0.0 } }
-    override val one: Buffer<Double> by lazy { Float64Buffer(size) { 1.0 } }
+public class Float64BufferField(public val size: Int) : ExtendedField<Buffer<Float64>>, Float64BufferOps() {
+    override val zero: Buffer<Float64> by lazy { Float64Buffer(size) { 0.0 } }
+    override val one: Buffer<Float64> by lazy { Float64Buffer(size) { 1.0 } }

-    override fun sinh(arg: Buffer<Double>): Float64Buffer = super<Float64BufferOps>.sinh(arg)
+    override fun sinh(arg: Buffer<Float64>): Float64Buffer = super<Float64BufferOps>.sinh(arg)

-    override fun cosh(arg: Buffer<Double>): Float64Buffer = super<Float64BufferOps>.cosh(arg)
+    override fun cosh(arg: Buffer<Float64>): Float64Buffer = super<Float64BufferOps>.cosh(arg)

-    override fun tanh(arg: Buffer<Double>): Float64Buffer = super<Float64BufferOps>.tanh(arg)
+    override fun tanh(arg: Buffer<Float64>): Float64Buffer = super<Float64BufferOps>.tanh(arg)

-    override fun asinh(arg: Buffer<Double>): Float64Buffer = super<Float64BufferOps>.asinh(arg)
+    override fun asinh(arg: Buffer<Float64>): Float64Buffer = super<Float64BufferOps>.asinh(arg)

-    override fun acosh(arg: Buffer<Double>): Float64Buffer = super<Float64BufferOps>.acosh(arg)
+    override fun acosh(arg: Buffer<Float64>): Float64Buffer = super<Float64BufferOps>.acosh(arg)

-    override fun atanh(arg: Buffer<Double>): Float64Buffer = super<Float64BufferOps>.atanh(arg)
+    override fun atanh(arg: Buffer<Float64>): Float64Buffer = super<Float64BufferOps>.atanh(arg)

-    override fun power(arg: Buffer<Double>, pow: Number): Float64Buffer = if (pow.isInteger()) {
+    override fun power(arg: Buffer<Float64>, pow: Number): Float64Buffer = if (pow.isInteger()) {
        arg.map { it.pow(pow.toInt()) }
    } else {
        arg.map {
@ -38,6 +39,6 @@ public class Float64BufferField(public val size: Int) : ExtendedField<Buffer<Dou
        }
    }

-    override fun unaryOperationFunction(operation: String): (arg: Buffer<Double>) -> Buffer<Double> =
+    override fun unaryOperationFunction(operation: String): (arg: Buffer<Float64>) -> Buffer<Float64> =
        super<ExtendedField>.unaryOperationFunction(operation)
 }
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Float64BufferOps.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Float64BufferOps.kt
@ -14,43 +14,43 @@ import kotlin.math.sqrt
 /**
 * [ExtendedFieldOps] over [Float64Buffer].
 */
-public abstract class Float64BufferOps : BufferAlgebra<Double, Float64Field>, ExtendedFieldOps<Buffer<Double>>,
-    Norm<Buffer<Double>, Double> {
+public abstract class Float64BufferOps : BufferAlgebra<Double, Float64Field>, ExtendedFieldOps<Buffer<Float64>>,
+    Norm<Buffer<Float64>, Double> {

    override val elementAlgebra: Float64Field get() = Float64Field

-    override val elementBufferFactory: MutableBufferFactory<Double> get() = elementAlgebra.bufferFactory
+    override val elementBufferFactory: MutableBufferFactory<Float64> get() = elementAlgebra.bufferFactory

    @Suppress("OVERRIDE_BY_INLINE")
    @OptIn(UnstableKMathAPI::class)
-    final override inline fun Buffer<Double>.map(block: Float64Field.(Double) -> Double): Float64Buffer =
+    final override inline fun Buffer<Float64>.map(block: Float64Field.(Double) -> Double): Float64Buffer =
        DoubleArray(size) { Float64Field.block(getDouble(it)) }.asBuffer()


    @OptIn(UnstableKMathAPI::class)
    @Suppress("OVERRIDE_BY_INLINE")
-    final override inline fun Buffer<Double>.mapIndexed(block: Float64Field.(index: Int, arg: Double) -> Double): Float64Buffer =
+    final override inline fun Buffer<Float64>.mapIndexed(block: Float64Field.(index: Int, arg: Double) -> Double): Float64Buffer =
        Float64Buffer(size) { Float64Field.block(it, getDouble(it)) }

    @OptIn(UnstableKMathAPI::class)
    @Suppress("OVERRIDE_BY_INLINE")
-    final override inline fun Buffer<Double>.zip(
-        other: Buffer<Double>,
+    final override inline fun Buffer<Float64>.zip(
+        other: Buffer<Float64>,
        block: Float64Field.(left: Double, right: Double) -> Double,
    ): Float64Buffer {
        require(size == other.size) { "Incompatible buffer sizes. left: ${size}, right: ${other.size}" }
        return Float64Buffer(size) { Float64Field.block(getDouble(it), other.getDouble(it)) }
    }

-    override fun unaryOperationFunction(operation: String): (arg: Buffer<Double>) -> Buffer<Double> =
+    override fun unaryOperationFunction(operation: String): (arg: Buffer<Float64>) -> Buffer<Float64> =
        super<ExtendedFieldOps>.unaryOperationFunction(operation)

-    override fun binaryOperationFunction(operation: String): (left: Buffer<Double>, right: Buffer<Double>) -> Buffer<Double> =
+    override fun binaryOperationFunction(operation: String): (left: Buffer<Float64>, right: Buffer<Float64>) -> Buffer<Float64> =
        super<ExtendedFieldOps>.binaryOperationFunction(operation)

-    override fun Buffer<Double>.unaryMinus(): Float64Buffer = map { -it }
+    override fun Buffer<Float64>.unaryMinus(): Float64Buffer = map { -it }

-    override fun add(left: Buffer<Double>, right: Buffer<Double>): Float64Buffer {
+    override fun add(left: Buffer<Float64>, right: Buffer<Float64>): Float64Buffer {
        require(right.size == left.size) {
            "The size of the first buffer ${left.size} should be the same as for second one: ${right.size} "
        }
@ -62,9 +62,9 @@ public abstract class Float64BufferOps : BufferAlgebra<Double, Float64Field>, Ex
        } else Float64Buffer(DoubleArray(left.size) { left[it] + right[it] })
    }

-    override fun Buffer<Double>.plus(arg: Buffer<Double>): Float64Buffer = add(this, arg)
+    override fun Buffer<Float64>.plus(arg: Buffer<Float64>): Float64Buffer = add(this, arg)

-    override fun Buffer<Double>.minus(arg: Buffer<Double>): Float64Buffer {
+    override fun Buffer<Float64>.minus(arg: Buffer<Float64>): Float64Buffer {
        require(arg.size == this.size) {
            "The size of the first buffer ${this.size} should be the same as for second one: ${arg.size} "
        }
@ -77,7 +77,7 @@ public abstract class Float64BufferOps : BufferAlgebra<Double, Float64Field>, Ex
    }

    //
-//    override fun multiply(a: Buffer<Double>, k: Number): RealBuffer {
+//    override fun multiply(a: Buffer<Float64>, k: Number): RealBuffer {
 //        val kValue = k.toDouble()
 //
 //        return if (a is RealBuffer) {
@ -86,7 +86,7 @@ public abstract class Float64BufferOps : BufferAlgebra<Double, Float64Field>, Ex
 //        } else RealBuffer(DoubleArray(a.size) { a[it] * kValue })
 //    }
 //
-//    override fun divide(a: Buffer<Double>, k: Number): RealBuffer {
+//    override fun divide(a: Buffer<Float64>, k: Number): RealBuffer {
 //        val kValue = k.toDouble()
 //
 //        return if (a is RealBuffer) {
@ -96,7 +96,7 @@ public abstract class Float64BufferOps : BufferAlgebra<Double, Float64Field>, Ex
 //    }

    @UnstableKMathAPI
-    override fun multiply(left: Buffer<Double>, right: Buffer<Double>): Float64Buffer {
+    override fun multiply(left: Buffer<Float64>, right: Buffer<Float64>): Float64Buffer {
        require(right.size == left.size) {
            "The size of the first buffer ${left.size} should be the same as for second one: ${right.size} "
        }
@ -108,7 +108,7 @@ public abstract class Float64BufferOps : BufferAlgebra<Double, Float64Field>, Ex
        } else Float64Buffer(DoubleArray(left.size) { left[it] * right[it] })
    }

-    override fun divide(left: Buffer<Double>, right: Buffer<Double>): Float64Buffer {
+    override fun divide(left: Buffer<Float64>, right: Buffer<Float64>): Float64Buffer {
        require(right.size == left.size) {
            "The size of the first buffer ${left.size} should be the same as for second one: ${right.size} "
        }
@ -120,39 +120,39 @@ public abstract class Float64BufferOps : BufferAlgebra<Double, Float64Field>, Ex
        } else Float64Buffer(DoubleArray(left.size) { left[it] / right[it] })
    }

-    override fun sin(arg: Buffer<Double>): Float64Buffer = arg.map { sin(it) }
+    override fun sin(arg: Buffer<Float64>): Float64Buffer = arg.map { sin(it) }

-    override fun cos(arg: Buffer<Double>): Float64Buffer = arg.map { cos(it) }
+    override fun cos(arg: Buffer<Float64>): Float64Buffer = arg.map { cos(it) }

-    override fun tan(arg: Buffer<Double>): Float64Buffer = arg.map { tan(it) }
+    override fun tan(arg: Buffer<Float64>): Float64Buffer = arg.map { tan(it) }

-    override fun asin(arg: Buffer<Double>): Float64Buffer = arg.map { asin(it) }
+    override fun asin(arg: Buffer<Float64>): Float64Buffer = arg.map { asin(it) }

-    override fun acos(arg: Buffer<Double>): Float64Buffer = arg.map { acos(it) }
+    override fun acos(arg: Buffer<Float64>): Float64Buffer = arg.map { acos(it) }

-    override fun atan(arg: Buffer<Double>): Float64Buffer = arg.map { atan(it) }
+    override fun atan(arg: Buffer<Float64>): Float64Buffer = arg.map { atan(it) }

-    override fun sinh(arg: Buffer<Double>): Float64Buffer = arg.map { sinh(it) }
+    override fun sinh(arg: Buffer<Float64>): Float64Buffer = arg.map { sinh(it) }

-    override fun cosh(arg: Buffer<Double>): Float64Buffer = arg.map { cosh(it) }
+    override fun cosh(arg: Buffer<Float64>): Float64Buffer = arg.map { cosh(it) }

-    override fun tanh(arg: Buffer<Double>): Float64Buffer = arg.map { tanh(it) }
+    override fun tanh(arg: Buffer<Float64>): Float64Buffer = arg.map { tanh(it) }

-    override fun asinh(arg: Buffer<Double>): Float64Buffer = arg.map { asinh(it) }
+    override fun asinh(arg: Buffer<Float64>): Float64Buffer = arg.map { asinh(it) }

-    override fun acosh(arg: Buffer<Double>): Float64Buffer = arg.map { acosh(it) }
+    override fun acosh(arg: Buffer<Float64>): Float64Buffer = arg.map { acosh(it) }

-    override fun atanh(arg: Buffer<Double>): Float64Buffer = arg.map { atanh(it) }
+    override fun atanh(arg: Buffer<Float64>): Float64Buffer = arg.map { atanh(it) }

-    override fun exp(arg: Buffer<Double>): Float64Buffer = arg.map { exp(it) }
+    override fun exp(arg: Buffer<Float64>): Float64Buffer = arg.map { exp(it) }

-    override fun ln(arg: Buffer<Double>): Float64Buffer = arg.map { ln(it) }
+    override fun ln(arg: Buffer<Float64>): Float64Buffer = arg.map { ln(it) }

-    override fun norm(arg: Buffer<Double>): Double = Float64L2Norm.norm(arg)
+    override fun norm(arg: Buffer<Float64>): Double = Float64L2Norm.norm(arg)

-    override fun scale(a: Buffer<Double>, value: Double): Float64Buffer = a.map { it * value }
+    override fun scale(a: Buffer<Float64>, value: Double): Float64Buffer = a.map { it * value }

-    override fun power(arg: Buffer<Double>, pow: Number): Buffer<Double> = if (pow is Int) {
+    override fun power(arg: Buffer<Float64>, pow: Number): Buffer<Float64> = if (pow is Int) {
        arg.map { it.pow(pow) }
    } else {
        arg.map { it.pow(pow.toDouble()) }
@ -161,11 +161,11 @@ public abstract class Float64BufferOps : BufferAlgebra<Double, Float64Field>, Ex
    public companion object : Float64BufferOps()
 }

-public object Float64L2Norm : Norm<Point<Double>, Double> {
-    override fun norm(arg: Point<Double>): Double = sqrt(arg.fold(0.0) { acc: Double, d: Double -> acc + d.pow(2) })
+public object Float64L2Norm : Norm<Point<Float64>, Double> {
+    override fun norm(arg: Point<Float64>): Double = sqrt(arg.fold(0.0) { acc: Double, d: Double -> acc + d.pow(2) })
 }

-public fun Float64BufferOps.sum(buffer: Buffer<Double>): Double = buffer.reduce(Double::plus)
+public fun Float64BufferOps.sum(buffer: Buffer<Float64>): Double = buffer.reduce(Double::plus)

 /**
 * Sum of elements using given [conversion]
@ -173,7 +173,7 @@ public fun Float64BufferOps.sum(buffer: Buffer<Double>): Double = buffer.reduce(
 public inline fun <T> Float64BufferOps.sumOf(buffer: Buffer<T>, conversion: (T) -> Double): Double =
    buffer.fold(0.0) { acc, value -> acc + conversion(value) }

-public fun Float64BufferOps.average(buffer: Buffer<Double>): Double = sum(buffer) / buffer.size
+public fun Float64BufferOps.average(buffer: Buffer<Float64>): Double = sum(buffer) / buffer.size

 /**
 * Average of elements using given [conversion]
@ -181,14 +181,14 @@ public fun Float64BufferOps.average(buffer: Buffer<Double>): Double = sum(buffer
 public inline fun <T> Float64BufferOps.averageOf(buffer: Buffer<T>, conversion: (T) -> Double): Double =
    sumOf(buffer, conversion) / buffer.size

-public fun Float64BufferOps.dispersion(buffer: Buffer<Double>): Double {
+public fun Float64BufferOps.dispersion(buffer: Buffer<Float64>): Double {
    val av = average(buffer)
    return buffer.fold(0.0) { acc, value -> acc + (value - av).pow(2) } / buffer.size
 }

-public fun Float64BufferOps.std(buffer: Buffer<Double>): Double = sqrt(dispersion(buffer))
+public fun Float64BufferOps.std(buffer: Buffer<Float64>): Double = sqrt(dispersion(buffer))

-public fun Float64BufferOps.covariance(x: Buffer<Double>, y: Buffer<Double>): Double {
+public fun Float64BufferOps.covariance(x: Buffer<Float64>, y: Buffer<Float64>): Double {
    require(x.size == y.size) { "Expected buffers of the same size, but x.size == ${x.size} and y.size == ${y.size}" }
    val xMean = average(x)
    val yMean = average(y)
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt
@ -4,6 +4,7 @@
 */
 package space.kscience.kmath.operations

+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.MutableBufferFactory
 import kotlin.math.pow as kpow

@ -66,8 +67,8 @@ public interface ExtendedField<T> : ExtendedFieldOps<T>, Field<T>, NumericAlgebr
 * A field for [Double] without boxing. Does not produce an appropriate field element.
 */
@Suppress("EXTENSION_SHADOWED_BY_MEMBER")
-public object Float64Field : ExtendedField<Double>, Norm<Double, Double>, ScaleOperations<Double> {
-    override val bufferFactory: MutableBufferFactory<Double> = MutableBufferFactory()
+public object Float64Field : ExtendedField<Float64>, Norm<Double, Double>, ScaleOperations<Float64> {
+    override val bufferFactory: MutableBufferFactory<Float64> = MutableBufferFactory()

    override val zero: Double get() = 0.0
    override val one: Double get() = 1.0
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt
@ -115,7 +115,7 @@ public fun <T> Buffer(
    size: Int,
    initializer: (Int) -> T,
 ): Buffer<T> = when (type.kType) {
-    typeOf<Double>() -> MutableBuffer.double(size) { initializer(it) as Double } as Buffer<T>
+    typeOf<Float64>() -> MutableBuffer.double(size) { initializer(it) as Double } as Buffer<T>
    typeOf<Short>() -> MutableBuffer.short(size) { initializer(it) as Short } as Buffer<T>
    typeOf<Int>() -> MutableBuffer.int(size) { initializer(it) as Int } as Buffer<T>
    typeOf<Long>() -> MutableBuffer.long(size) { initializer(it) as Long } as Buffer<T>
@ -134,7 +134,7 @@ public inline fun <reified T> Buffer(size: Int, initializer: (Int) -> T): Buffer
    //code duplication here because we want to inline initializers
    val type = safeTypeOf<T>()
    return when (type.kType) {
-        typeOf<Double>() -> MutableBuffer.double(size) { initializer(it) as Double } as Buffer<T>
+        typeOf<Float64>() -> MutableBuffer.double(size) { initializer(it) as Double } as Buffer<T>
        typeOf<Short>() -> MutableBuffer.short(size) { initializer(it) as Short } as Buffer<T>
        typeOf<Int>() -> MutableBuffer.int(size) { initializer(it) as Int } as Buffer<T>
        typeOf<Long>() -> MutableBuffer.long(size) { initializer(it) as Long } as Buffer<T>
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Float64Buffer.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Float64Buffer.kt
@ -14,7 +14,7 @@ import kotlin.jvm.JvmInline
 * @property array the underlying array.
 */
@JvmInline
-public value class Float64Buffer(public val array: DoubleArray) : PrimitiveBuffer<Double> {
+public value class Float64Buffer(public val array: DoubleArray) : PrimitiveBuffer<Float64> {

    override val size: Int get() = array.size

@ -54,7 +54,7 @@ public fun Float64Buffer(vararg doubles: Double): Float64Buffer = Float64Buffer(
 /**
 * Returns a new [DoubleArray] containing all the elements of this [Buffer].
 */
-public fun Buffer<Double>.toDoubleArray(): DoubleArray = when (this) {
+public fun Buffer<Float64>.toDoubleArray(): DoubleArray = when (this) {
    is Float64Buffer -> array
    else -> DoubleArray(size, ::get)
 }
@ -62,7 +62,7 @@ public fun Buffer<Double>.toDoubleArray(): DoubleArray = when (this) {
 /**
 * Represent this buffer as [Float64Buffer]. Does not guarantee that changes in the original buffer are reflected on this buffer.
 */
-public fun Buffer<Double>.toFloat64Buffer(): Float64Buffer = when (this) {
+public fun Buffer<Float64>.toFloat64Buffer(): Float64Buffer = when (this) {
    is Float64Buffer -> this
    else -> DoubleArray(size, ::get).asBuffer()
 }
@ -79,5 +79,5 @@ public fun DoubleArray.asBuffer(): Float64Buffer = Float64Buffer(this)
 public fun interface Float64BufferTransform : BufferTransform<Double, Double> {
    public fun transform(arg: Float64Buffer): Float64Buffer

-    override fun transform(arg: Buffer<Double>): Float64Buffer = arg.toFloat64Buffer()
+    override fun transform(arg: Buffer<Float64>): Float64Buffer = arg.toFloat64Buffer()
 }
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MutableBuffer.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MutableBuffer.kt
@ -99,7 +99,7 @@ public inline fun <T> MutableBuffer(
    typeOf<Int32>() -> MutableBuffer.int(size) { initializer(it) as Int32 } as MutableBuffer<T>
    typeOf<Int64>() -> MutableBuffer.long(size) { initializer(it) as Int64 } as MutableBuffer<T>
    typeOf<Float>() -> MutableBuffer.float(size) { initializer(it) as Float } as MutableBuffer<T>
-    typeOf<Double>() -> MutableBuffer.double(size) { initializer(it) as Double } as MutableBuffer<T>
+    typeOf<Float64>() -> MutableBuffer.double(size) { initializer(it) as Double } as MutableBuffer<T>
    //TODO add unsigned types
    else -> MutableListBuffer(MutableList(size, initializer))
 }
--- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferPrimitiveAccess.kt
+++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferPrimitiveAccess.kt
@ -11,7 +11,7 @@ import space.kscience.kmath.UnstableKMathAPI
 * Non-boxing access to primitive [Double]
 */
@UnstableKMathAPI
-public fun Buffer<Double>.getDouble(index: Int): Double = if (this is BufferView) {
+public fun Buffer<Float64>.getDouble(index: Int): Double = if (this is BufferView) {
    val originIndex = originIndex(index)
    if (originIndex >= 0) {
        origin.getDouble(originIndex)
--- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/ExpressionFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/ExpressionFieldTest.kt
@ -6,6 +6,7 @@
 package space.kscience.kmath.expressions

 import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
 import kotlin.test.Test
 import kotlin.test.assertEquals
 import kotlin.test.assertFails
@ -37,7 +38,7 @@ class ExpressionFieldTest {

    @Test
    fun valueExpression() {
-        val expressionBuilder: FunctionalExpressionField<Double, *>.() -> Expression<Double> = {
+        val expressionBuilder: FunctionalExpressionField<Double, *>.() -> Expression<Float64> = {
            val x by binding
            x * x + 2 * x + one
        }
--- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/SimpleAutoDiffTest.kt
+++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/SimpleAutoDiffTest.kt
@ -8,6 +8,7 @@ package space.kscience.kmath.expressions
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.bindSymbol
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.asBuffer
 import kotlin.math.E
 import kotlin.math.PI
@ -21,18 +22,18 @@ internal class SimpleAutoDiffTest {

    fun dx(
        xBinding: Pair<Symbol, Double>,
-        body: SimpleAutoDiffField<Double, Float64Field>.(x: AutoDiffValue<Double>) -> AutoDiffValue<Double>,
-    ): DerivationResult<Double> = Float64Field.simpleAutoDiff(xBinding) { body(bindSymbol(xBinding.first)) }
+        body: SimpleAutoDiffField<Double, Float64Field>.(x: AutoDiffValue<Float64>) -> AutoDiffValue<Float64>,
+    ): DerivationResult<Float64> = Float64Field.simpleAutoDiff(xBinding) { body(bindSymbol(xBinding.first)) }

    fun dxy(
        xBinding: Pair<Symbol, Double>,
        yBinding: Pair<Symbol, Double>,
-        body: SimpleAutoDiffField<Double, Float64Field>.(x: AutoDiffValue<Double>, y: AutoDiffValue<Double>) -> AutoDiffValue<Double>,
-    ): DerivationResult<Double> = Float64Field.simpleAutoDiff(xBinding, yBinding) {
+        body: SimpleAutoDiffField<Double, Float64Field>.(x: AutoDiffValue<Float64>, y: AutoDiffValue<Float64>) -> AutoDiffValue<Float64>,
+    ): DerivationResult<Float64> = Float64Field.simpleAutoDiff(xBinding, yBinding) {
        body(bindSymbol(xBinding.first), bindSymbol(yBinding.first))
    }

-    fun diff(block: SimpleAutoDiffField<Double, Float64Field>.() -> AutoDiffValue<Double>): SimpleAutoDiffExpression<Double, Float64Field> {
+    fun diff(block: SimpleAutoDiffField<Double, Float64Field>.() -> AutoDiffValue<Float64>): SimpleAutoDiffExpression<Double, Float64Field> {
        return SimpleAutoDiffExpression(Float64Field, block)
    }

--- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixTest.kt
+++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixTest.kt
@ -10,6 +10,7 @@ import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.structures.Float64
 import kotlin.test.Test
 import kotlin.test.assertEquals
 import kotlin.test.assertTrue
@ -38,7 +39,7 @@ class MatrixTest {

    @Test
    fun testMatrixExtension() = Double.algebra.linearSpace.run {
-        val transitionMatrix: Matrix<Double> = VirtualMatrix(6, 6) { row, col ->
+        val transitionMatrix: Matrix<Float64> = VirtualMatrix(6, 6) { row, col ->
            when {
                col == 0 -> .50
                row + 1 == col -> .50
@ -47,7 +48,7 @@ class MatrixTest {
            }
        }

-        infix fun Matrix<Double>.pow(power: Int): Matrix<Double> {
+        infix fun Matrix<Float64>.pow(power: Int): Matrix<Float64> {
            var res = this
            repeat(power - 1) {
                res = res dot this@pow
--- a/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/linear/Float64ParallelLinearSpace.kt
+++ b/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/linear/Float64ParallelLinearSpace.kt
@ -11,6 +11,7 @@ import space.kscience.kmath.operations.Float64BufferOps
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer
 import space.kscience.kmath.structures.asBuffer
 import java.util.stream.IntStream
@ -25,7 +26,7 @@ public object Float64ParallelLinearSpace : LinearSpace<Double, Float64Field> {
        rows: Int,
        columns: Int,
        initializer: Float64Field.(i: Int, j: Int) -> Double,
-    ): Matrix<Double> {
+    ): Matrix<Float64> {
        val shape = ShapeND(rows, columns)
        val indexer = BufferAlgebraND.defaultIndexerBuilder(shape)

@ -43,29 +44,29 @@ public object Float64ParallelLinearSpace : LinearSpace<Double, Float64Field> {
    override fun buildVector(size: Int, initializer: Float64Field.(Int) -> Double): Float64Buffer =
        IntStream.range(0, size).parallel().mapToDouble { Float64Field.initializer(it) }.toArray().asBuffer()

-    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = Floa64FieldOpsND {
+    override fun Matrix<Float64>.unaryMinus(): Matrix<Float64> = Floa64FieldOpsND {
        asND().map { -it }.as2D()
    }

-    override fun Matrix<Double>.plus(other: Matrix<Double>): Matrix<Double> = Floa64FieldOpsND {
+    override fun Matrix<Float64>.plus(other: Matrix<Float64>): Matrix<Float64> = Floa64FieldOpsND {
        require(shape == other.shape) { "Shape mismatch on Matrix::plus. Expected $shape but found ${other.shape}" }
        asND().plus(other.asND()).as2D()
    }

-    override fun Matrix<Double>.minus(other: Matrix<Double>): Matrix<Double> = Floa64FieldOpsND {
+    override fun Matrix<Float64>.minus(other: Matrix<Float64>): Matrix<Float64> = Floa64FieldOpsND {
        require(shape == other.shape) { "Shape mismatch on Matrix::minus. Expected $shape but found ${other.shape}" }
        asND().minus(other.asND()).as2D()
    }

    // Create a continuous in-memory representation of this vector for better memory layout handling
-    private fun Buffer<Double>.linearize() = if (this is Float64Buffer) {
+    private fun Buffer<Float64>.linearize() = if (this is Float64Buffer) {
        this.array
    } else {
        DoubleArray(size) { get(it) }
    }

    @OptIn(PerformancePitfall::class)
-    override fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> {
+    override fun Matrix<Float64>.dot(other: Matrix<Float64>): Matrix<Float64> {
        require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
        val rows = this@dot.rows.map { it.linearize() }
        val columns = other.columns.map { it.linearize() }
@ -82,7 +83,7 @@ public object Float64ParallelLinearSpace : LinearSpace<Double, Float64Field> {
    }

    @OptIn(PerformancePitfall::class)
-    override fun Matrix<Double>.dot(vector: Point<Double>): Float64Buffer {
+    override fun Matrix<Float64>.dot(vector: Point<Float64>): Float64Buffer {
        require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
        val rows = this@dot.rows.map { it.linearize() }
        val indices = 0 until this.colNum
@ -97,27 +98,27 @@ public object Float64ParallelLinearSpace : LinearSpace<Double, Float64Field> {

    }

-    override fun Matrix<Double>.times(value: Double): Matrix<Double> = Floa64FieldOpsND {
+    override fun Matrix<Float64>.times(value: Double): Matrix<Float64> = Floa64FieldOpsND {
        asND().map { it * value }.as2D()
    }

-    public override fun Point<Double>.plus(other: Point<Double>): Float64Buffer = Float64BufferOps.run {
+    public override fun Point<Float64>.plus(other: Point<Float64>): Float64Buffer = Float64BufferOps.run {
        this@plus + other
    }

-    public override fun Point<Double>.minus(other: Point<Double>): Float64Buffer = Float64BufferOps.run {
+    public override fun Point<Float64>.minus(other: Point<Float64>): Float64Buffer = Float64BufferOps.run {
        this@minus - other
    }

-    public override fun Point<Double>.times(value: Double): Float64Buffer = Float64BufferOps.run {
+    public override fun Point<Float64>.times(value: Double): Float64Buffer = Float64BufferOps.run {
        scale(this@times, value)
    }

-    public operator fun Point<Double>.div(value: Double): Float64Buffer = Float64BufferOps.run {
+    public operator fun Point<Float64>.div(value: Double): Float64Buffer = Float64BufferOps.run {
        scale(this@div, 1.0 / value)
    }

-    public override fun Double.times(v: Point<Double>): Float64Buffer = v * this
+    public override fun Double.times(v: Point<Float64>): Float64Buffer = v * this


 }
--- a/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/structures/parallelMutableBuffer.kt
+++ b/kmath-core/src/jvmMain/kotlin/space/kscience/kmath/structures/parallelMutableBuffer.kt
@ -30,7 +30,7 @@ public fun <T> MutableBuffer.Companion.parallel(
        .asBuffer() as MutableBuffer<T>

    typeOf<Float>() -> Float32Buffer(size) { initializer(it) as Float } as MutableBuffer<T>
-    typeOf<Double>() -> IntStream.range(0, size).parallel().mapToDouble { initializer(it) as Float64 }.toArray()
+    typeOf<Float64>() -> IntStream.range(0, size).parallel().mapToDouble { initializer(it) as Float64 }.toArray()
        .asBuffer() as MutableBuffer<T>
    //TODO add unsigned types
    else -> IntStream.range(0, size).parallel().mapToObj { initializer(it) }.collect(Collectors.toList<T>())
--- a/kmath-core/src/jvmTest/kotlin/space/kscience/kmath/linear/ParallelMatrixTest.kt
+++ b/kmath-core/src/jvmTest/kotlin/space/kscience/kmath/linear/ParallelMatrixTest.kt
@ -9,6 +9,7 @@ import space.kscience.kmath.PerformancePitfall
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
 import kotlin.test.Test
 import kotlin.test.assertEquals
 import kotlin.test.assertTrue
@ -37,7 +38,7 @@ class ParallelMatrixTest {

    @Test
    fun testMatrixExtension() = Float64Field.linearSpace.parallel {
-        val transitionMatrix: Matrix<Double> = VirtualMatrix(6, 6) { row, col ->
+        val transitionMatrix: Matrix<Float64> = VirtualMatrix(6, 6) { row, col ->
            when {
                col == 0 -> .50
                row + 1 == col -> .50
@ -46,7 +47,7 @@ class ParallelMatrixTest {
            }
        }

-        infix fun Matrix<Double>.pow(power: Int): Matrix<Double> {
+        infix fun Matrix<Float64>.pow(power: Int): Matrix<Float64> {
            var res = this
            repeat(power - 1) {
                res = res dot this@pow
--- a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingChain.kt
+++ b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingChain.kt
@ -40,13 +40,13 @@ public interface BlockingBufferChain<out T> : BlockingChain<T>, BufferChain<T> {
 }


-public suspend inline fun <reified T : Any> Chain<T>.nextBuffer(size: Int): Buffer<T> = if (this is BufferChain) {
+public suspend inline fun <reified T> Chain<T>.nextBuffer(size: Int): Buffer<T> = if (this is BufferChain) {
    nextBuffer(size)
 } else {
    Buffer(size) { next() }
 }

-public inline fun <reified T : Any> BlockingChain<T>.nextBufferBlocking(
+public inline fun <reified T> BlockingChain<T>.nextBufferBlocking(
    size: Int,
 ): Buffer<T> = if (this is BlockingBufferChain) {
    nextBufferBlocking(size)
--- a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingDoubleChain.kt
+++ b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingDoubleChain.kt
@ -5,12 +5,13 @@

 package space.kscience.kmath.chains

+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer

 /**
 * Chunked, specialized chain for double values, which supports blocking [nextBlocking] operation
 */
-public interface BlockingDoubleChain : BlockingBufferChain<Double> {
+public interface BlockingDoubleChain : BlockingBufferChain<Float64> {

    /**
     * Returns an [DoubleArray] chunk of [size] values of [next].
--- a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/Chain.kt
+++ b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/Chain.kt
@ -23,7 +23,7 @@ public interface Chain<out T> : Flow<T> {
    public suspend fun next(): T

    /**
-     * Create a copy of current chain state. Consuming resulting chain does not affect initial chain.
+     * Create a copy of the current chain state. Consuming the resulting chain does not affect the initial chain.
     */
    public suspend fun fork(): Chain<T>

@ -47,7 +47,7 @@ public class SimpleChain<out R>(private val gen: suspend () -> R) : Chain<R> {
 /**
 * A stateless Markov chain
 */
-public class MarkovChain<out R : Any>(private val seed: suspend () -> R, private val gen: suspend (R) -> R) : Chain<R> {
+public class MarkovChain<out R>(private val seed: suspend () -> R, private val gen: suspend (R) -> R) : Chain<R> {
    private val mutex: Mutex = Mutex()
    private var value: R? = null

@ -71,8 +71,8 @@ public class MarkovChain<out R : Any>(private val seed: suspend () -> R, private
 */
 public class StatefulChain<S, out R>(
    private val state: S,
-    private val seed: S.() -> R,
-    private val forkState: ((S) -> S),
+    private val seed: suspend S.() -> R,
+    private val forkState: suspend ((S) -> S),
    private val gen: suspend S.(R) -> R,
 ) : Chain<R> {
    private val mutex: Mutex = Mutex()
@ -97,6 +97,11 @@ public class ConstantChain<out T>(public val value: T) : Chain<T> {
    override suspend fun fork(): Chain<T> = this
 }

+/**
+ * Discard a fixed number of samples
+ */
+public suspend inline fun <reified T> Chain<T>.discard(number: Int): Chain<T> = apply { nextBuffer<T>(number) }
+
 /**
 * Map the chain result using suspended transformation. Initial chain result can no longer be safely consumed
 * since mapped chain consumes tokens. Accepts regular transformation function.
--- a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/BufferFlow.kt
+++ b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/BufferFlow.kt
@ -17,6 +17,7 @@ import space.kscience.kmath.chains.BlockingDoubleChain
 import space.kscience.kmath.operations.Group
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.BufferFactory
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer

 /**
@ -56,7 +57,7 @@ public fun <T> Flow<T>.chunked(bufferSize: Int, bufferFactory: BufferFactory<T>)
 /**
 * Specialized flow chunker for real buffer
 */
-public fun Flow<Double>.chunked(bufferSize: Int): Flow<Float64Buffer> = flow {
+public fun Flow<Float64>.chunked(bufferSize: Int): Flow<Float64Buffer> = flow {
    require(bufferSize > 0) { "Resulting chunk size must be more than zero" }

    if (this@chunked is BlockingDoubleChain) {
--- a/kmath-ejml/build.gradle.kts
+++ b/kmath-ejml/build.gradle.kts
@ -1,15 +1,20 @@
 plugins {
-    id("space.kscience.gradle.jvm")
+    id("space.kscience.gradle.mpp")
 }

 val ejmlVerision = "0.43.1"

-dependencies {
-    api("org.ejml:ejml-ddense:$ejmlVerision")
-    api("org.ejml:ejml-fdense:$ejmlVerision")
-    api("org.ejml:ejml-dsparse:$ejmlVerision")
-    api("org.ejml:ejml-fsparse:$ejmlVerision")
-    api(projects.kmathCore)
+kscience {
+    jvm()
+    jvmMain {
+        api(projects.kmathCore)
+        api(projects.kmathComplex)
+        api("org.ejml:ejml-all:$ejmlVerision")
+    }
+
+    jvmTest {
+        implementation(projects.testUtils)
+    }
 }

 readme {
@ -30,12 +35,4 @@ readme {
        id = "ejml-linear-space",
        ref = "src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt"
    ) { "LinearSpace implementations." }
-}
-
-//kotlin.sourceSets.main {
-//    val codegen by tasks.creating {
-//        ejmlCodegen(kotlin.srcDirs.first().absolutePath + "/space/kscience/kmath/ejml/_generated.kt")
-//    }
-//
-//    kotlin.srcDirs(files().builtBy(codegen))
-//}
+}
--- a/kmath-ejml/src/jvmMain/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt
+++ b/kmath-ejml/src/jvmMain/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt
@ -6,9 +6,13 @@
 package space.kscience.kmath.ejml

 import space.kscience.kmath.UnstableKMathAPI
-import space.kscience.kmath.linear.*
-import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.linear.Inverted
+import space.kscience.kmath.linear.LinearSpace
+import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.linear.Point
+import space.kscience.kmath.nd.Structure2D
 import space.kscience.kmath.operations.Ring
+import space.kscience.kmath.structures.Float64

 /**
 * [LinearSpace] implementation specialized for a certain EJML type.
@ -38,6 +42,6 @@ public abstract class EjmlLinearSpace<T : Any, out A : Ring<T>, out M : org.ejml
    public abstract override fun buildVector(size: Int, initializer: A.(Int) -> T): EjmlVector<T, M>

    @UnstableKMathAPI
-    public fun EjmlMatrix<T, *>.inverted(): Matrix<Double> =
-        computeAttribute(this, Float64Field.linearSpace.Inverted)!!
+    public fun Structure2D<T>.inverted(): Matrix<Float64> =
+        computeAttribute(this, Inverted()) ?: error("Can't invert matrix")
 }
--- a/kmath-ejml/src/jvmMain/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt
+++ b/kmath-ejml/src/jvmMain/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt
--- a/kmath-ejml/src/jvmMain/kotlin/space/kscience/kmath/ejml/EjmlVector.kt
+++ b/kmath-ejml/src/jvmMain/kotlin/space/kscience/kmath/ejml/EjmlVector.kt
--- a/kmath-ejml/src/jvmMain/kotlin/space/kscience/kmath/ejml/implementations.kt
+++ b/kmath-ejml/src/jvmMain/kotlin/space/kscience/kmath/ejml/implementations.kt
@ -12,6 +12,8 @@ import org.ejml.dense.row.CommonOps_DDRM
 import org.ejml.dense.row.CommonOps_FDRM
 import org.ejml.dense.row.factory.DecompositionFactory_DDRM
 import org.ejml.dense.row.factory.DecompositionFactory_FDRM
+import org.ejml.interfaces.decomposition.EigenDecomposition_F32
+import org.ejml.interfaces.decomposition.EigenDecomposition_F64
 import org.ejml.sparse.FillReducing
 import org.ejml.sparse.csc.CommonOps_DSCC
 import org.ejml.sparse.csc.CommonOps_FSCC
@ -19,6 +21,8 @@ import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC
 import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC
 import space.kscience.attributes.SafeType
 import space.kscience.attributes.safeTypeOf
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.complex.Complex
 import space.kscience.kmath.linear.*
 import space.kscience.kmath.linear.Matrix
 import space.kscience.kmath.nd.Structure2D
@ -27,9 +31,15 @@ import space.kscience.kmath.operations.Float32Field
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.Float32
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.IntBuffer
 import space.kscience.kmath.structures.asBuffer

+/**
+ * Copy EJML [Complex_F64] into KMath [Complex]
+ */
+public fun Complex_F64.toKMathComplex(): Complex = Complex(real, imaginary)
+
 /**
 * [EjmlVector] specialization for [Double].
 */
@ -79,16 +89,16 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
     */
    override val elementAlgebra: Float64Field get() = Float64Field

-    override val type: SafeType<Double> get() = safeTypeOf()
+    override val type: SafeType<Float64> get() = safeTypeOf()

    @Suppress("UNCHECKED_CAST")
-    override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixRMaj> = when {
+    override fun Matrix<Float64>.toEjml(): EjmlDoubleMatrix<DMatrixRMaj> = when {
        this is EjmlDoubleMatrix<*> && origin is DMatrixRMaj -> this as EjmlDoubleMatrix<DMatrixRMaj>
        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
    }

    @Suppress("UNCHECKED_CAST")
-    override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixRMaj> = when {
+    override fun Point<Float64>.toEjml(): EjmlDoubleVector<DMatrixRMaj> = when {
        this is EjmlDoubleVector<*> && origin is DMatrixRMaj -> this as EjmlDoubleVector<DMatrixRMaj>
        else -> EjmlDoubleVector(DMatrixRMaj(size, 1).also {
            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
@ -115,21 +125,21 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
    private fun <T : DMatrix> T.wrapVector() = EjmlDoubleVector(this)

-    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
+    override fun Matrix<Float64>.unaryMinus(): Matrix<Float64> = this * elementAlgebra { -one }

-    override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+    override fun Matrix<Float64>.dot(other: Matrix<Float64>): EjmlDoubleMatrix<DMatrixRMaj> {
        val out = DMatrixRMaj(1, 1)
        CommonOps_DDRM.mult(toEjml().origin, other.toEjml().origin, out)
        return out.wrapMatrix()
    }

-    override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+    override fun Matrix<Float64>.dot(vector: Point<Float64>): EjmlDoubleVector<DMatrixRMaj> {
        val out = DMatrixRMaj(1, 1)
        CommonOps_DDRM.mult(toEjml().origin, vector.toEjml().origin, out)
        return out.wrapVector()
    }

-    override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+    override operator fun Matrix<Float64>.minus(other: Matrix<Float64>): EjmlDoubleMatrix<DMatrixRMaj> {
        val out = DMatrixRMaj(1, 1)

        CommonOps_DDRM.add(
@ -143,19 +153,19 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
        return out.wrapMatrix()
    }

-    override operator fun Matrix<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixRMaj> {
+    override operator fun Matrix<Float64>.times(value: Double): EjmlDoubleMatrix<DMatrixRMaj> {
        val res = DMatrixRMaj(1, 1)
        CommonOps_DDRM.scale(value, toEjml().origin, res)
        return res.wrapMatrix()
    }

-    override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixRMaj> {
+    override fun Point<Float64>.unaryMinus(): EjmlDoubleVector<DMatrixRMaj> {
        val res = DMatrixRMaj(1, 1)
        CommonOps_DDRM.changeSign(toEjml().origin, res)
        return res.wrapVector()
    }

-    override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+    override fun Matrix<Float64>.plus(other: Matrix<Float64>): EjmlDoubleMatrix<DMatrixRMaj> {
        val out = DMatrixRMaj(1, 1)

        CommonOps_DDRM.add(
@ -169,7 +179,7 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
        return out.wrapMatrix()
    }

-    override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+    override fun Point<Float64>.plus(other: Point<Float64>): EjmlDoubleVector<DMatrixRMaj> {
        val out = DMatrixRMaj(1, 1)

        CommonOps_DDRM.add(
@ -183,7 +193,7 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
        return out.wrapVector()
    }

-    override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+    override fun Point<Float64>.minus(other: Point<Float64>): EjmlDoubleVector<DMatrixRMaj> {
        val out = DMatrixRMaj(1, 1)

        CommonOps_DDRM.add(
@ -197,18 +207,19 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
        return out.wrapVector()
    }

-    override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> = m * this
+    override fun Double.times(m: Matrix<Float64>): EjmlDoubleMatrix<DMatrixRMaj> = m * this

-    override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixRMaj> {
+    override fun Point<Float64>.times(value: Double): EjmlDoubleVector<DMatrixRMaj> {
        val res = DMatrixRMaj(1, 1)
        CommonOps_DDRM.scale(value, toEjml().origin, res)
        return res.wrapVector()
    }

-    override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixRMaj> = v * this
+    override fun Double.times(v: Point<Float64>): EjmlDoubleVector<DMatrixRMaj> = v * this

-    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Double>, attribute: A): V? {
-        val origin = structure.toEjml().origin
+    @OptIn(UnstableKMathAPI::class)
+    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Float64>, attribute: A): V? {
+        val origin: DMatrixRMaj = structure.toEjml().origin

        val raw: Any? = when (attribute) {
            Inverted -> {
@ -218,28 +229,28 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
            }

            Determinant -> CommonOps_DDRM.det(origin)
-            SVD -> object : SingularValueDecomposition<Double> {
+            SVD -> object : SingularValueDecomposition<Float64> {
                val ejmlSvd by lazy {
                    DecompositionFactory_DDRM
                        .svd(origin.numRows, origin.numCols, true, true, false)
                        .apply { decompose(origin.copy()) }
                }
-                override val u: Matrix<Double> get() = ejmlSvd.getU(null, false).wrapMatrix()
+                override val u: Matrix<Float64> get() = ejmlSvd.getU(null, false).wrapMatrix()

-                override val s: Matrix<Double> get() = ejmlSvd.getW(null).wrapMatrix()
-                override val v: Matrix<Double> get() = ejmlSvd.getV(null, false).wrapMatrix()
-                override val singularValues: Point<Double> get() = ejmlSvd.singularValues.asBuffer()
+                override val s: Matrix<Float64> get() = ejmlSvd.getW(null).wrapMatrix()
+                override val v: Matrix<Float64> get() = ejmlSvd.getV(null, false).wrapMatrix()
+                override val singularValues: Point<Float64> get() = ejmlSvd.singularValues.asBuffer()

            }

-            QR -> object : QRDecomposition<Double> {
+            QR -> object : QRDecomposition<Float64> {
                val ejmlQr by lazy { DecompositionFactory_DDRM.qr().apply { decompose(origin.copy()) } }
-                override val q: Matrix<Double> get() = ejmlQr.getQ(null, false).wrapMatrix()
-                override val r: Matrix<Double> get() = ejmlQr.getR(null, false).wrapMatrix()
+                override val q: Matrix<Float64> get() = ejmlQr.getQ(null, false).wrapMatrix()
+                override val r: Matrix<Float64> get() = ejmlQr.getR(null, false).wrapMatrix()
            }

-            Cholesky -> object : CholeskyDecomposition<Double> {
-                override val l: Matrix<Double> by lazy {
+            Cholesky -> object : CholeskyDecomposition<Float64> {
+                override val l: Matrix<Float64> by lazy {
                    val cholesky =
                        DecompositionFactory_DDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }

@ -247,20 +258,49 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
                }
            }

-            LUP -> object : LupDecomposition<Double> {
+            LUP -> object : LupDecomposition<Float64> {
                private val lup by lazy {
                    DecompositionFactory_DDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
                }

-                override val l: Matrix<Double>
+                override val l: Matrix<Float64>
                    get() = lup.getLower(null).wrapMatrix().withAttribute(LowerTriangular)

-
-                override val u: Matrix<Double>
+                override val u: Matrix<Float64>
                    get() = lup.getUpper(null).wrapMatrix().withAttribute(UpperTriangular)
+
                override val pivot: IntBuffer get() = lup.getRowPivotV(null).asBuffer()
            }

+            EIG -> {
+                check(origin.numCols == origin.numRows) { "Eigenvalue decomposition requires symmetric matrix" }
+                object : EigenDecomposition<Float64> {
+                    val cmEigen: EigenDecomposition_F64<DMatrixRMaj> by lazy {
+                        DecompositionFactory_DDRM.eig(origin.numRows, true).apply { decompose(origin) }
+                    }
+
+                    override val v: Matrix<Float64> by lazy {
+                        val eigenvectors = List(origin.numRows) { cmEigen.getEigenVector(it) }.filterNotNull()
+                        buildMatrix(eigenvectors.size, origin.numCols) { row, column ->
+                            eigenvectors[row][column]
+                        }
+                    }
+
+                    override val d: Matrix<Float64> by lazy {
+                        val eigenvalues = List(origin.numRows) { cmEigen.getEigenvalue(it) }
+
+                        buildMatrix(origin.numRows, origin.numCols) { row, column ->
+                            when (row) {
+                                column -> eigenvalues[row].real
+                                column - 1 -> eigenvalues[row].imaginary
+                                column + 1 -> -eigenvalues[row].imaginary
+                                else -> 0.0
+                            }
+                        }
+                    }
+                }
+            }
+
            else -> null
        }

@ -275,7 +315,7 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
     * @param b n by p matrix.
     * @return the solution for *x* that is n by p.
     */
-    public fun solve(a: Matrix<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+    public fun solve(a: Matrix<Float64>, b: Matrix<Float64>): EjmlDoubleMatrix<DMatrixRMaj> {
        val res = DMatrixRMaj(1, 1)
        CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res)
        return res.wrapMatrix()
@ -288,7 +328,7 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
     * @param b n by p vector.
     * @return the solution for *x* that is n by p.
     */
-    public fun solve(a: Matrix<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+    public fun solve(a: Matrix<Float64>, b: Point<Float64>): EjmlDoubleVector<DMatrixRMaj> {
        val res = DMatrixRMaj(1, 1)
        CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res)
        return EjmlDoubleVector(res)
@ -434,6 +474,7 @@ public object EjmlLinearSpaceFDRM : EjmlLinearSpace<Float, Float32Field, FMatrix

    override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixRMaj> = v * this

+    @OptIn(UnstableKMathAPI::class)
    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Float32>, attribute: A): V? {
        val origin = structure.toEjml().origin

@ -488,6 +529,35 @@ public object EjmlLinearSpaceFDRM : EjmlLinearSpace<Float, Float32Field, FMatrix
                override val pivot: IntBuffer get() = lup.getRowPivotV(null).asBuffer()
            }

+            EIG -> {
+                check(origin.numCols == origin.numRows) { "Eigenvalue decomposition requires symmetric matrix" }
+                object : EigenDecomposition<Float32> {
+                    val cmEigen: EigenDecomposition_F32<FMatrixRMaj> by lazy {
+                        DecompositionFactory_FDRM.eig(origin.numRows, true).apply { decompose(origin) }
+                    }
+
+                    override val v by lazy {
+                        val eigenvectors: List<FMatrixRMaj> = List(origin.numRows) { cmEigen.getEigenVector(it) }
+                        buildMatrix(origin.numRows, origin.numCols) { row, column ->
+                            eigenvectors[row][column]
+                        }
+                    }
+
+                    override val d: Matrix<Float32> by lazy {
+                        val eigenvalues = List(origin.numRows) { cmEigen.getEigenvalue(it) }
+
+                        buildMatrix(origin.numRows, origin.numCols) { row, column ->
+                            when (row) {
+                                column -> eigenvalues[row].real
+                                column - 1 -> eigenvalues[row].imaginary
+                                column + 1 -> -eigenvalues[row].imaginary
+                                else -> 0f
+                            }
+                        }
+                    }
+                }
+            }
+
            else -> null
        }

@ -533,16 +603,16 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
     */
    override val elementAlgebra: Float64Field get() = Float64Field

-    override val type: SafeType<Double> get() = safeTypeOf()
+    override val type: SafeType<Float64> get() = safeTypeOf()

    @Suppress("UNCHECKED_CAST")
-    override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixSparseCSC> = when {
+    override fun Matrix<Float64>.toEjml(): EjmlDoubleMatrix<DMatrixSparseCSC> = when {
        this is EjmlDoubleMatrix<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleMatrix<DMatrixSparseCSC>
        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
    }

    @Suppress("UNCHECKED_CAST")
-    override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixSparseCSC> = when {
+    override fun Point<Float64>.toEjml(): EjmlDoubleVector<DMatrixSparseCSC> = when {
        this is EjmlDoubleVector<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleVector<DMatrixSparseCSC>
        else -> EjmlDoubleVector(DMatrixSparseCSC(size, 1).also {
            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
@ -569,21 +639,21 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
    private fun <T : DMatrix> T.wrapVector() = EjmlDoubleVector(this)

-    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
+    override fun Matrix<Float64>.unaryMinus(): Matrix<Float64> = this * elementAlgebra { -one }

-    override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+    override fun Matrix<Float64>.dot(other: Matrix<Float64>): EjmlDoubleMatrix<DMatrixSparseCSC> {
        val out = DMatrixSparseCSC(1, 1)
        CommonOps_DSCC.mult(toEjml().origin, other.toEjml().origin, out)
        return out.wrapMatrix()
    }

-    override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+    override fun Matrix<Float64>.dot(vector: Point<Float64>): EjmlDoubleVector<DMatrixSparseCSC> {
        val out = DMatrixSparseCSC(1, 1)
        CommonOps_DSCC.mult(toEjml().origin, vector.toEjml().origin, out)
        return out.wrapVector()
    }

-    override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+    override operator fun Matrix<Float64>.minus(other: Matrix<Float64>): EjmlDoubleMatrix<DMatrixSparseCSC> {
        val out = DMatrixSparseCSC(1, 1)

        CommonOps_DSCC.add(
@ -599,19 +669,19 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
        return out.wrapMatrix()
    }

-    override operator fun Matrix<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixSparseCSC> {
+    override operator fun Matrix<Float64>.times(value: Double): EjmlDoubleMatrix<DMatrixSparseCSC> {
        val res = DMatrixSparseCSC(1, 1)
        CommonOps_DSCC.scale(value, toEjml().origin, res)
        return res.wrapMatrix()
    }

-    override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixSparseCSC> {
+    override fun Point<Float64>.unaryMinus(): EjmlDoubleVector<DMatrixSparseCSC> {
        val res = DMatrixSparseCSC(1, 1)
        CommonOps_DSCC.changeSign(toEjml().origin, res)
        return res.wrapVector()
    }

-    override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+    override fun Matrix<Float64>.plus(other: Matrix<Float64>): EjmlDoubleMatrix<DMatrixSparseCSC> {
        val out = DMatrixSparseCSC(1, 1)

        CommonOps_DSCC.add(
@ -627,7 +697,7 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
        return out.wrapMatrix()
    }

-    override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+    override fun Point<Float64>.plus(other: Point<Float64>): EjmlDoubleVector<DMatrixSparseCSC> {
        val out = DMatrixSparseCSC(1, 1)

        CommonOps_DSCC.add(
@ -643,7 +713,7 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
        return out.wrapVector()
    }

-    override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+    override fun Point<Float64>.minus(other: Point<Float64>): EjmlDoubleVector<DMatrixSparseCSC> {
        val out = DMatrixSparseCSC(1, 1)

        CommonOps_DSCC.add(
@ -659,17 +729,17 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
        return out.wrapVector()
    }

-    override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> = m * this
+    override fun Double.times(m: Matrix<Float64>): EjmlDoubleMatrix<DMatrixSparseCSC> = m * this

-    override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixSparseCSC> {
+    override fun Point<Float64>.times(value: Double): EjmlDoubleVector<DMatrixSparseCSC> {
        val res = DMatrixSparseCSC(1, 1)
        CommonOps_DSCC.scale(value, toEjml().origin, res)
        return res.wrapVector()
    }

-    override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> = v * this
+    override fun Double.times(v: Point<Float64>): EjmlDoubleVector<DMatrixSparseCSC> = v * this

-    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Double>, attribute: A): V? {
+    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Float64>, attribute: A): V? {
        val origin = structure.toEjml().origin

        val raw: Any? = when (attribute) {
@ -681,16 +751,16 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri

            Determinant -> CommonOps_DSCC.det(origin)

-            QR -> object : QRDecomposition<Double> {
+            QR -> object : QRDecomposition<Float64> {
                val ejmlQr by lazy {
                    DecompositionFactory_DSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
                }
-                override val q: Matrix<Double> get() = ejmlQr.getQ(null, false).wrapMatrix()
-                override val r: Matrix<Double> get() = ejmlQr.getR(null, false).wrapMatrix()
+                override val q: Matrix<Float64> get() = ejmlQr.getQ(null, false).wrapMatrix()
+                override val r: Matrix<Float64> get() = ejmlQr.getR(null, false).wrapMatrix()
            }

-            Cholesky -> object : CholeskyDecomposition<Double> {
-                override val l: Matrix<Double> by lazy {
+            Cholesky -> object : CholeskyDecomposition<Float64> {
+                override val l: Matrix<Float64> by lazy {
                    val cholesky =
                        DecompositionFactory_DSCC.cholesky().apply { decompose(origin.copy()) }

@ -698,16 +768,16 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
                }
            }

-            LUP -> object : LupDecomposition<Double> {
+            LUP -> object : LupDecomposition<Float64> {
                private val lup by lazy {
                    DecompositionFactory_DSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
                }

-                override val l: Matrix<Double>
+                override val l: Matrix<Float64>
                    get() = lup.getLower(null).wrapMatrix().withAttribute(LowerTriangular)


-                override val u: Matrix<Double>
+                override val u: Matrix<Float64>
                    get() = lup.getUpper(null).wrapMatrix().withAttribute(UpperTriangular)
                override val pivot: IntBuffer get() = lup.getRowPivotV(null).asBuffer()
            }
@ -726,7 +796,7 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
     * @param b n by p matrix.
     * @return the solution for *x* that is n by p.
     */
-    public fun solve(a: Matrix<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+    public fun solve(a: Matrix<Float64>, b: Matrix<Float64>): EjmlDoubleMatrix<DMatrixSparseCSC> {
        val res = DMatrixSparseCSC(1, 1)
        CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res)
        return res.wrapMatrix()
@ -739,7 +809,7 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
     * @param b n by p vector.
     * @return the solution for *x* that is n by p.
     */
-    public fun solve(a: Matrix<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+    public fun solve(a: Matrix<Float64>, b: Point<Float64>): EjmlDoubleVector<DMatrixSparseCSC> {
        val res = DMatrixSparseCSC(1, 1)
        CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res)
        return EjmlDoubleVector(res)
--- a/kmath-ejml/src/jvmTest/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt
+++ b/kmath-ejml/src/jvmTest/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt
@ -17,12 +17,16 @@ import space.kscience.kmath.linear.*
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.nd.toArray
 import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.testutils.assertStructureEquals
 import kotlin.random.Random
 import kotlin.random.asJavaRandom
 import kotlin.test.*

 internal fun <T : Any> assertMatrixEquals(expected: StructureND<T>, actual: StructureND<T>) {
-    assertTrue { StructureND.contentEquals(expected, actual) }
+    expected.elements().forEach { (index, value) ->
+        assertEquals(value, actual[index], "Structure element with index ${index.toList()} should be equal to $value but is ${actual[index]}")
+    }
 }

@OptIn(UnstableKMathAPI::class)
@ -63,7 +67,7 @@ internal class EjmlMatrixTest {
        val w = EjmlDoubleMatrix(m)
        val det: Double = w.getOrComputeAttribute(Determinant) ?: fail()
        assertEquals(CommonOps_DDRM.det(m), det)
-        val lup: LupDecomposition<Double> = w.getOrComputeAttribute(LUP) ?: fail()
+        val lup: LupDecomposition<Float64> = w.getOrComputeAttribute(LUP) ?: fail()

        val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows, m.numCols)
            .also { it.decompose(m.copy()) }
@ -104,4 +108,15 @@ internal class EjmlMatrixTest {

        assertTrue { StructureND.contentEquals(one(dim, dim), res, 1e-3) }
    }
+
+    @Test
+    fun eigenValueDecomposition() = EjmlLinearSpaceDDRM {
+        val dim = 46
+        val u = buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+        val matrix = buildMatrix(dim, dim) { row, col ->
+            if (row >= col) u[row, col] else u[col, row]
+        }
+        val eigen = matrix.getOrComputeAttribute(EIG) ?: fail()
+        assertStructureEquals(matrix, eigen.v dot eigen.d dot eigen.v.transposed())
+    }
 }
--- a/kmath-ejml/src/jvmTest/kotlin/space/kscience/kmath/ejml/EjmlVectorTest.kt
+++ b/kmath-ejml/src/jvmTest/kotlin/space/kscience/kmath/ejml/EjmlVectorTest.kt
--- a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt
+++ b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt
@ -9,12 +9,13 @@ import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.linear.Point
 import space.kscience.kmath.operations.Float64L2Norm
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.MutableBuffer.Companion.double
 import space.kscience.kmath.structures.asBuffer
 import space.kscience.kmath.structures.indices
 import kotlin.math.pow

-public typealias DoubleVector = Point<Double>
+public typealias DoubleVector = Point<Float64>

 public fun DoubleVector(vararg doubles: Double): DoubleVector = doubles.asBuffer()

@ -41,7 +42,7 @@ public operator fun DoubleVector.plus(number: Number): DoubleVector = map { it +

 public operator fun Number.plus(vector: DoubleVector): DoubleVector = vector + this

-public operator fun DoubleVector.unaryMinus(): Buffer<Double> = map { -it }
+public operator fun DoubleVector.unaryMinus(): Buffer<Float64> = map { -it }

 public operator fun DoubleVector.minus(other: DoubleVector): DoubleVector {
    require(size == other.size) { "Vector size $size expected but ${other.size} found" }
--- a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt
+++ b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt
@ -15,6 +15,7 @@ import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.algebra
 import space.kscience.kmath.operations.asIterable
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer
 import kotlin.math.pow

@ -30,7 +31,7 @@ import kotlin.math.pow
 *  Functions that help create a real (Double) matrix
 */

-public typealias RealMatrix = Matrix<Double>
+public typealias RealMatrix = Matrix<Float64>

 public fun realMatrix(rowNum: Int, colNum: Int, initializer: Float64Field.(i: Int, j: Int) -> Double): RealMatrix =
    Double.algebra.linearSpace.buildMatrix(rowNum, colNum, initializer)
@ -114,7 +115,7 @@ public operator fun RealMatrix.minus(other: RealMatrix): RealMatrix =
 *  Operations on columns
 */

-public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer<Double>) -> Double): RealMatrix =
+public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer<Float64>) -> Double): RealMatrix =
    Double.algebra.linearSpace.buildMatrix(rowNum, colNum + 1) { row, col ->
        if (col < colNum)
            get(row, col)
--- a/Show More
+++ b/Show More