Merge branch 'dev' into feature/tensorflow

# Conflicts:
#	kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt

CHANGELOG.md

@@ -17,6 +17,7 @@
 - `BigInt` operation performance improvement and fixes by @zhelenskiy (#328)
 - Integration between `MST` and Symja `IExpr`
 - Complex power
+- Separate methods for UInt, Int and Number powers. NaN safety.
 
 ### Changed
 - Exponential operations merged with hyperbolic functions

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt

@@ -13,6 +13,7 @@ import space.kscience.kmath.commons.linear.CMLinearSpace
 import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
 import space.kscience.kmath.linear.invoke
 import space.kscience.kmath.linear.linearSpace
+import space.kscience.kmath.multik.multikAlgebra
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.structures.Buffer
 import kotlin.random.Random
@@ -58,6 +59,16 @@ internal class DotBenchmark {
         blackhole.consume(matrix1 dot matrix2)
     }
 
+//    @Benchmark
+//    fun tensorDot(blackhole: Blackhole) = with(Double.tensorAlgebra) {
+//        blackhole.consume(matrix1 dot matrix2)
+//    }
+
+    @Benchmark
+    fun multikDot(blackhole: Blackhole) = with(Double.multikAlgebra) {
+        blackhole.consume(matrix1 dot matrix2)
+    }
+
     @Benchmark
     fun bufferedDot(blackhole: Blackhole) = with(DoubleField.linearSpace(Buffer.Companion::auto)) {
         blackhole.consume(matrix1 dot matrix2)

examples/src/main/kotlin/space/kscience/kmath/tensors/neuralNetwork.kt

@@ -163,7 +163,7 @@ class NeuralNetwork(private val layers: List<Layer>) {
             for ((xBatch, yBatch) in iterBatch(xTrain, yTrain)) {
                 train(xBatch, yBatch)
             }
-            println("Accuracy:${accuracy(yTrain, predict(xTrain).argMax(1, true))}")
+            println("Accuracy:${accuracy(yTrain, predict(xTrain).argMax(1, true).asDouble())}")
         }
     }
 
@@ -230,7 +230,7 @@ fun main() = BroadcastDoubleTensorAlgebra {
     val prediction = model.predict(xTest)
 
     // process raw prediction via argMax
-    val predictionLabels = prediction.argMax(1, true)
+    val predictionLabels = prediction.argMax(1, true).asDouble()
 
     // find out accuracy
     val acc = accuracy(yTest, predictionLabels)

kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/ComplexFieldND.kt

@@ -18,7 +18,7 @@ import kotlin.contracts.contract
  */
 @OptIn(UnstableKMathAPI::class)
 public sealed class ComplexFieldOpsND : BufferedFieldOpsND<Complex, ComplexField>(ComplexField.bufferAlgebra),
-    ScaleOperations<StructureND<Complex>>, ExtendedFieldOps<StructureND<Complex>> {
+    ScaleOperations<StructureND<Complex>>, ExtendedFieldOps<StructureND<Complex>>, PowerOperations<StructureND<Complex>> {
 
     override fun StructureND<Complex>.toBufferND(): BufferND<Complex> = when (this) {
         is BufferND -> this
@@ -33,9 +33,6 @@ public sealed class ComplexFieldOpsND : BufferedFieldOpsND<Complex, ComplexField
     override fun scale(a: StructureND<Complex>, value: Double): BufferND<Complex> =
         mapInline(a.toBufferND()) { it * value }
 
-    override fun power(arg: StructureND<Complex>, pow: Number): BufferND<Complex> =
-        mapInline(arg.toBufferND()) { power(it, pow) }
-
     override fun exp(arg: StructureND<Complex>): BufferND<Complex> = mapInline(arg.toBufferND()) { exp(it) }
     override fun ln(arg: StructureND<Complex>): BufferND<Complex> = mapInline(arg.toBufferND()) { ln(it) }
 
@@ -53,6 +50,9 @@ public sealed class ComplexFieldOpsND : BufferedFieldOpsND<Complex, ComplexField
     override fun acosh(arg: StructureND<Complex>): BufferND<Complex> = mapInline(arg.toBufferND()) { acosh(it) }
     override fun atanh(arg: StructureND<Complex>): BufferND<Complex> = mapInline(arg.toBufferND()) { atanh(it) }
 
+    override fun power(arg: StructureND<Complex>, pow: Number): StructureND<Complex> =
+        mapInline(arg.toBufferND()) { power(it,pow) }
+
     public companion object : ComplexFieldOpsND()
 }
 
@@ -63,7 +63,8 @@ public val ComplexField.bufferAlgebra: BufferFieldOps<Complex, ComplexField>
 
 @OptIn(UnstableKMathAPI::class)
 public class ComplexFieldND(override val shape: Shape) :
-    ComplexFieldOpsND(), FieldND<Complex, ComplexField>, NumbersAddOps<StructureND<Complex>> {
+    ComplexFieldOpsND(), FieldND<Complex, ComplexField>,
+    NumbersAddOps<StructureND<Complex>> {
 
     override fun number(value: Number): BufferND<Complex> {
         val d = value.toDouble() // minimize conversions

kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt

@@ -252,7 +252,7 @@ public class SimpleAutoDiffExpression<T : Any, F : Field<T>>(
  * Generate [AutoDiffProcessor] for [SimpleAutoDiffExpression]
  */
 public fun <T : Any, F : Field<T>> simpleAutoDiff(
-    field: F
+    field: F,
 ): AutoDiffProcessor<T, AutoDiffValue<T>, SimpleAutoDiffField<T, F>> =
     AutoDiffProcessor { function ->
         SimpleAutoDiffExpression(field, function)
@@ -272,8 +272,8 @@ public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.sqrt(x: Aut
 public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.pow(
     x: AutoDiffValue<T>,
     y: Double,
-): AutoDiffValue<T> = derive(const { power(x.value, y) }) { z ->
+): AutoDiffValue<T> = derive(const { x.value.pow(y)}) { z ->
-    x.d += z.d * y * power(x.value, y - 1)
+    x.d += z.d * y * x.value.pow(y - 1)
 }
 
 public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.pow(

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt

@@ -35,7 +35,7 @@ public interface WithShape {
  * @param T the type of ND-structure element.
  * @param C the type of the element context.
  */
-public interface AlgebraND<T, out C : Algebra<T>> {
+public interface AlgebraND<T, out C : Algebra<T>>: Algebra<StructureND<T>> {
     /**
      * The algebra over elements of ND structure.
      */

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt

@@ -53,21 +53,28 @@ public interface BufferAlgebraND<T, out A : Algebra<T>> : AlgebraND<T, A> {
 
 public inline fun <T, A : Algebra<T>> BufferAlgebraND<T, A>.mapInline(
     arg: BufferND<T>,
-    crossinline transform: A.(T) -> T
+    crossinline transform: A.(T) -> T,
 ): BufferND<T> {
     val indexes = arg.indices
-    return BufferND(indexes, bufferAlgebra.mapInline(arg.buffer, transform))
+    val buffer = arg.buffer
+    return BufferND(
+        indexes,
+        bufferAlgebra.run {
+            bufferFactory(buffer.size) { elementAlgebra.transform(buffer[it]) }
+        }
+    )
 }
 
 internal inline fun <T, A : Algebra<T>> BufferAlgebraND<T, A>.mapIndexedInline(
     arg: BufferND<T>,
-    crossinline transform: A.(index: IntArray, arg: T) -> T
+    crossinline transform: A.(index: IntArray, arg: T) -> T,
 ): BufferND<T> {
     val indexes = arg.indices
+    val buffer = arg.buffer
     return BufferND(
         indexes,
-        bufferAlgebra.mapIndexedInline(arg.buffer) { offset, value ->
+        bufferAlgebra.run {
-            transform(indexes.index(offset), value)
+            bufferFactory(buffer.size) { elementAlgebra.transform(indexes.index(it), buffer[it]) }
         }
     )
 }
@@ -75,35 +82,42 @@ internal inline fun <T, A : Algebra<T>> BufferAlgebraND<T, A>.mapIndexedInline(
 internal inline fun <T, A : Algebra<T>> BufferAlgebraND<T, A>.zipInline(
     l: BufferND<T>,
     r: BufferND<T>,
-    crossinline block: A.(l: T, r: T) -> T
+    crossinline block: A.(l: T, r: T) -> T,
 ): BufferND<T> {
     require(l.indices == r.indices) { "Zip requires the same shapes, but found ${l.shape} on the left and ${r.shape} on the right" }
     val indexes = l.indices
-    return BufferND(indexes, bufferAlgebra.zipInline(l.buffer, r.buffer, block))
+    val lbuffer = l.buffer
+    val rbuffer = r.buffer
+    return BufferND(
+        indexes,
+        bufferAlgebra.run {
+            bufferFactory(lbuffer.size) { elementAlgebra.block(lbuffer[it], rbuffer[it]) }
+        }
+    )
 }
 
 @OptIn(PerformancePitfall::class)
 public open class BufferedGroupNDOps<T, out A : Group<T>>(
     override val bufferAlgebra: BufferAlgebra<T, A>,
-    override val indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder
+    override val indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
 ) : GroupOpsND<T, A>, BufferAlgebraND<T, A> {
     override fun StructureND<T>.unaryMinus(): StructureND<T> = map { -it }
 }
 
 public open class BufferedRingOpsND<T, out A : Ring<T>>(
     bufferAlgebra: BufferAlgebra<T, A>,
-    indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder
+    indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
 ) : BufferedGroupNDOps<T, A>(bufferAlgebra, indexerBuilder), RingOpsND<T, A>
 
 public open class BufferedFieldOpsND<T, out A : Field<T>>(
     bufferAlgebra: BufferAlgebra<T, A>,
-    indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder
+    indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
 ) : BufferedRingOpsND<T, A>(bufferAlgebra, indexerBuilder), FieldOpsND<T, A> {
 
     public constructor(
         elementAlgebra: A,
         bufferFactory: BufferFactory<T>,
-        indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder
+        indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
     ) : this(BufferFieldOps(elementAlgebra, bufferFactory), indexerBuilder)
 
     @OptIn(PerformancePitfall::class)
@@ -117,11 +131,11 @@ public val <T, A : Field<T>> BufferAlgebra<T, A>.nd: BufferedFieldOpsND<T, A> ge
 
 public fun <T, A : Algebra<T>> BufferAlgebraND<T, A>.structureND(
     vararg shape: Int,
-    initializer: A.(IntArray) -> T
+    initializer: A.(IntArray) -> T,
 ): BufferND<T> = structureND(shape, initializer)
 
 public fun <T, EA : Algebra<T>, A> A.structureND(
-    initializer: EA.(IntArray) -> T
+    initializer: EA.(IntArray) -> T,
 ): BufferND<T> where A : BufferAlgebraND<T, EA>, A : WithShape = structureND(shape, initializer)
 
 //// group factories

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt

@@ -5,7 +5,6 @@
 
 package space.kscience.kmath.nd
 
-import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.BufferFactory
 import space.kscience.kmath.structures.MutableBuffer
@@ -27,11 +26,6 @@ public open class BufferND<out T>(
 
     override val shape: IntArray get() = indices.shape
 
-    @PerformancePitfall
-    override fun elements(): Sequence<Pair<IntArray, T>> = indices.asSequence().map {
-        it to this[it]
-    }
-
     override fun toString(): String = StructureND.toString(this)
 }
 

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/DoubleFieldND.kt

@@ -11,7 +11,7 @@ import space.kscience.kmath.operations.*
 import space.kscience.kmath.structures.DoubleBuffer
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract
-import kotlin.math.pow
+import kotlin.math.pow as kpow
 
 public class DoubleBufferND(
     indexes: ShapeIndexer,
@@ -30,9 +30,9 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
         }
     }
 
-    private inline fun mapInline(
+    protected inline fun mapInline(
         arg: DoubleBufferND,
-        transform: (Double) -> Double
+        transform: (Double) -> Double,
     ): DoubleBufferND {
         val indexes = arg.indices
         val array = arg.buffer.array
@@ -42,7 +42,7 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
     private inline fun zipInline(
         l: DoubleBufferND,
         r: DoubleBufferND,
-        block: (l: Double, r: Double) -> Double
+        block: (l: Double, r: Double) -> Double,
     ): DoubleBufferND {
         require(l.indices == r.indices) { "Zip requires the same shapes, but found ${l.shape} on the left and ${r.shape} on the right" }
         val indexes = l.indices
@@ -60,7 +60,7 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
     override fun zip(
         left: StructureND<Double>,
         right: StructureND<Double>,
-        transform: DoubleField.(Double, Double) -> Double
+        transform: DoubleField.(Double, Double) -> Double,
     ): BufferND<Double> = zipInline(left.toBufferND(), right.toBufferND()) { l, r -> DoubleField.transform(l, r) }
 
     override fun structureND(shape: Shape, initializer: DoubleField.(IntArray) -> Double): DoubleBufferND {
@@ -81,8 +81,8 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
 
     override fun StructureND<Double>.unaryMinus(): DoubleBufferND = mapInline(toBufferND()) { -it }
 
-    override fun StructureND<Double>.div(other: StructureND<Double>): DoubleBufferND =
+    override fun StructureND<Double>.div(arg: StructureND<Double>): DoubleBufferND =
-        zipInline(toBufferND(), other.toBufferND()) { l, r -> l / r }
+        zipInline(toBufferND(), arg.toBufferND()) { l, r -> l / r }
 
     override fun divide(left: StructureND<Double>, right: StructureND<Double>): DoubleBufferND =
         zipInline(left.toBufferND(), right.toBufferND()) { l: Double, r: Double -> l / r }
@@ -101,8 +101,8 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
     override fun StructureND<Double>.minus(arg: StructureND<Double>): DoubleBufferND =
         zipInline(toBufferND(), arg.toBufferND()) { l: Double, r: Double -> l - r }
 
-    override fun StructureND<Double>.times(other: StructureND<Double>): DoubleBufferND =
+    override fun StructureND<Double>.times(arg: StructureND<Double>): DoubleBufferND =
-        zipInline(toBufferND(), other.toBufferND()) { l: Double, r: Double -> l * r }
+        zipInline(toBufferND(), arg.toBufferND()) { l: Double, r: Double -> l * r }
 
     override fun StructureND<Double>.times(k: Number): DoubleBufferND =
         mapInline(toBufferND()) { it * k.toDouble() }
@@ -123,9 +123,6 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
     override fun scale(a: StructureND<Double>, value: Double): DoubleBufferND =
         mapInline(a.toBufferND()) { it * value }
 
-    override fun power(arg: StructureND<Double>, pow: Number): DoubleBufferND =
-        mapInline(arg.toBufferND()) { it.pow(pow.toDouble()) }
-
     override fun exp(arg: StructureND<Double>): DoubleBufferND =
         mapInline(arg.toBufferND()) { kotlin.math.exp(it) }
 
@@ -173,7 +170,38 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
 
 @OptIn(UnstableKMathAPI::class)
 public class DoubleFieldND(override val shape: Shape) :
-    DoubleFieldOpsND(), FieldND<Double, DoubleField>, NumbersAddOps<StructureND<Double>> {
+    DoubleFieldOpsND(), FieldND<Double, DoubleField>, NumbersAddOps<StructureND<Double>>,
+    ExtendedField<StructureND<Double>> {
+
+    override fun power(arg: StructureND<Double>, pow: UInt): DoubleBufferND = mapInline(arg.toBufferND()) {
+        it.kpow(pow.toInt())
+    }
+
+    override fun power(arg: StructureND<Double>, pow: Int): DoubleBufferND = mapInline(arg.toBufferND()) {
+        it.kpow(pow)
+    }
+
+    override fun power(arg: StructureND<Double>, pow: Number): DoubleBufferND = if(pow.isInteger()){
+        power(arg, pow.toInt())
+    } else {
+        val dpow = pow.toDouble()
+        mapInline(arg.toBufferND()) {
+            if (it < 0) throw IllegalArgumentException("Can't raise negative $it to a fractional power")
+            else it.kpow(dpow)
+        }
+    }
+
+    override fun sinh(arg: StructureND<Double>): DoubleBufferND = super<DoubleFieldOpsND>.sinh(arg)
+
+    override fun cosh(arg: StructureND<Double>): DoubleBufferND = super<DoubleFieldOpsND>.cosh(arg)
+
+    override fun tanh(arg: StructureND<Double>): DoubleBufferND = super<DoubleFieldOpsND>.tan(arg)
+
+    override fun asinh(arg: StructureND<Double>): DoubleBufferND = super<DoubleFieldOpsND>.asinh(arg)
+
+    override fun acosh(arg: StructureND<Double>): DoubleBufferND = super<DoubleFieldOpsND>.acosh(arg)
+
+    override fun atanh(arg: StructureND<Double>): DoubleBufferND = super<DoubleFieldOpsND>.atanh(arg)
 
     override fun number(value: Number): DoubleBufferND {
         val d = value.toDouble() // minimize conversions

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt

@@ -6,6 +6,8 @@
 package space.kscience.kmath.operations
 
 import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.Ring.Companion.optimizedPower
 
 /**
  * Stub for DSL the [Algebra] is.
@@ -99,6 +101,14 @@ public interface Algebra<T> {
      */
     public fun binaryOperation(operation: String, left: T, right: T): T =
         binaryOperationFunction(operation)(left, right)
+
+    /**
+     * Export an algebra element, so it could be accessed even after algebra scope is closed.
+     * This method must be used on algebras where data is stored externally or any local algebra state is used.
+     * By default (if not overridden), exports the object itself.
+     */
+    @UnstableKMathAPI
+    public fun export(arg: T): T = arg
 }
 
 public fun <T> Algebra<T>.bindSymbolOrNull(symbol: Symbol): T? = bindSymbolOrNull(symbol.identity)
@@ -162,6 +172,7 @@ public interface GroupOps<T> : Algebra<T> {
      * @return the difference.
      */
     public operator fun T.minus(arg: T): T = add(this, -arg)
+
     // Dynamic dispatch of operations
     override fun unaryOperationFunction(operation: String): (arg: T) -> T = when (operation) {
         PLUS_OPERATION -> { arg -> +arg }
@@ -247,6 +258,40 @@ public interface Ring<T> : Group<T>, RingOps<T> {
      * The neutral element of multiplication
      */
     public val one: T
+
+    /**
+     * Raises [arg] to the integer power [pow].
+     */
+    public fun power(arg: T, pow: UInt): T = optimizedPower(arg, pow)
+
+    public companion object{
+        /**
+         * Raises [arg] to the non-negative integer power [exponent].
+         *
+         * Special case: 0 ^ 0 is 1.
+         *
+         * @receiver the algebra to provide multiplication.
+         * @param arg the base.
+         * @param exponent the exponent.
+         * @return the base raised to the power.
+         * @author Evgeniy Zhelenskiy
+         */
+        internal fun <T> Ring<T>.optimizedPower(arg: T, exponent: UInt): T = when {
+            arg == zero && exponent > 0U -> zero
+            arg == one -> arg
+            arg == -one -> powWithoutOptimization(arg, exponent % 2U)
+            else -> powWithoutOptimization(arg, exponent)
+        }
+
+        private fun <T> Ring<T>.powWithoutOptimization(base: T, exponent: UInt): T = when (exponent) {
+            0U -> one
+            1U -> base
+            else -> {
+                val pre = powWithoutOptimization(base, exponent shr 1).let { it * it }
+                if (exponent and 1U == 0U) pre else pre * base
+            }
+        }
+    }
 }
 
 /**
@@ -270,10 +315,10 @@ public interface FieldOps<T> : RingOps<T> {
      * Division of two elements.
      *
      * @receiver the dividend.
-     * @param other the divisor.
+     * @param arg the divisor.
      * @return the quotient.
      */
-    public operator fun T.div(other: T): T = divide(this, other)
+    public operator fun T.div(arg: T): T = divide(this, arg)
 
     override fun binaryOperationFunction(operation: String): (left: T, right: T) -> T = when (operation) {
         DIV_OPERATION -> ::divide
@@ -297,4 +342,24 @@ public interface FieldOps<T> : RingOps<T> {
  */
 public interface Field<T> : Ring<T>, FieldOps<T>, ScaleOperations<T>, NumericAlgebra<T> {
     override fun number(value: Number): T = scale(one, value.toDouble())
+
+    public fun power(arg: T, pow: Int): T = optimizedPower(arg, pow)
+
+    public companion object{
+        /**
+         * Raises [arg] to the integer power [exponent].
+         *
+         * Special case: 0 ^ 0 is 1.
+         *
+         * @receiver the algebra to provide multiplication and division.
+         * @param arg the base.
+         * @param exponent the exponent.
+         * @return the base raised to the power.
+         * @author Iaroslav Postovalov, Evgeniy Zhelenskiy
+         */
+        private fun <T> Field<T>.optimizedPower(arg: T, exponent: Int): T = when {
+            exponent < 0 -> one / (this as Ring<T>).optimizedPower(arg, if (exponent == Int.MIN_VALUE) Int.MAX_VALUE.toUInt().inc() else (-exponent).toUInt())
+            else -> (this as Ring<T>).optimizedPower(arg, exponent.toUInt())
+        }
+    }
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.operations
 
+import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.BufferFactory
 import space.kscience.kmath.structures.DoubleBuffer
@@ -34,11 +35,13 @@ public interface BufferAlgebra<T, out A : Algebra<T>> : Algebra<Buffer<T>> {
     public fun Buffer<T>.zip(other: Buffer<T>, block: A.(left: T, right: T) -> T): Buffer<T> =
         zipInline(this, other, block)
 
+    @UnstableKMathAPI
     override fun unaryOperationFunction(operation: String): (arg: Buffer<T>) -> Buffer<T> {
         val operationFunction = elementAlgebra.unaryOperationFunction(operation)
         return { arg -> bufferFactory(arg.size) { operationFunction(arg[it]) } }
     }
 
+    @UnstableKMathAPI
     override fun binaryOperationFunction(operation: String): (left: Buffer<T>, right: Buffer<T>) -> Buffer<T> {
         val operationFunction = elementAlgebra.binaryOperationFunction(operation)
         return { left, right ->
@@ -50,7 +53,7 @@ public interface BufferAlgebra<T, out A : Algebra<T>> : Algebra<Buffer<T>> {
 /**
  * Inline map
  */
-public inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.mapInline(
+private inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.mapInline(
     buffer: Buffer<T>,
     crossinline block: A.(T) -> T
 ): Buffer<T> = bufferFactory(buffer.size) { elementAlgebra.block(buffer[it]) }
@@ -58,7 +61,7 @@ public inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.mapInline(
 /**
  * Inline map
  */
-public inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.mapIndexedInline(
+private inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.mapIndexedInline(
     buffer: Buffer<T>,
     crossinline block: A.(index: Int, arg: T) -> T
 ): Buffer<T> = bufferFactory(buffer.size) { elementAlgebra.block(it, buffer[it]) }
@@ -66,7 +69,7 @@ public inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.mapIndexedInline(
 /**
  * Inline zip
  */
-public inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.zipInline(
+private inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.zipInline(
     l: Buffer<T>,
     r: Buffer<T>,
     crossinline block: A.(l: T, r: T) -> T
@@ -126,7 +129,7 @@ public fun <T, A : ExponentialOperations<T>> BufferAlgebra<T, A>.atanh(arg: Buff
     mapInline(arg) { atanh(it) }
 
 public fun <T, A : PowerOperations<T>> BufferAlgebra<T, A>.pow(arg: Buffer<T>, pow: Number): Buffer<T> =
-    mapInline(arg) { power(it, pow) }
+    mapInline(arg) {it.pow(pow) }
 
 
 public open class BufferRingOps<T, A: Ring<T>>(
@@ -138,9 +141,11 @@ public open class BufferRingOps<T, A: Ring<T>>(
     override fun multiply(left: Buffer<T>, right: Buffer<T>): Buffer<T> = zipInline(left, right) { l, r -> l * r }
     override fun Buffer<T>.unaryMinus(): Buffer<T> = map { -it }
 
+    @UnstableKMathAPI
     override fun unaryOperationFunction(operation: String): (arg: Buffer<T>) -> Buffer<T> =
         super<BufferAlgebra>.unaryOperationFunction(operation)
 
+    @UnstableKMathAPI
     override fun binaryOperationFunction(operation: String): (left: Buffer<T>, right: Buffer<T>) -> Buffer<T> =
         super<BufferAlgebra>.binaryOperationFunction(operation)
 }
@@ -160,6 +165,7 @@ public open class BufferFieldOps<T, A : Field<T>>(
     override fun scale(a: Buffer<T>, value: Double): Buffer<T> = a.map { scale(it, value) }
     override fun Buffer<T>.unaryMinus(): Buffer<T> = map { -it }
 
+    @UnstableKMathAPI
     override fun binaryOperationFunction(operation: String): (left: Buffer<T>, right: Buffer<T>) -> Buffer<T> =
         super<BufferRingOps>.binaryOperationFunction(operation)
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferField.kt

@@ -5,6 +5,8 @@
 
 package space.kscience.kmath.operations
 
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.DoubleField.pow
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
 
@@ -27,7 +29,20 @@ public class DoubleBufferField(public val size: Int) : ExtendedField<Buffer<Doub
 
     override fun acosh(arg: Buffer<Double>): DoubleBuffer = super<DoubleBufferOps>.acosh(arg)
 
-    override fun atanh(arg: Buffer<Double>): DoubleBuffer= super<DoubleBufferOps>.atanh(arg)
+    override fun atanh(arg: Buffer<Double>): DoubleBuffer = super<DoubleBufferOps>.atanh(arg)
+
+    override fun power(arg: Buffer<Double>, pow: Number): DoubleBuffer = if (pow.isInteger()) {
+        arg.mapInline { it.pow(pow.toInt()) }
+    } else {
+        arg.mapInline {
+            if(it<0) throw IllegalArgumentException("Negative argument $it could not be raised to the fractional power")
+            it.pow(pow.toDouble())
+        }
+    }
+
+    @UnstableKMathAPI
+    override fun unaryOperationFunction(operation: String): (arg: Buffer<Double>) -> Buffer<Double> =
+        super<ExtendedField>.unaryOperationFunction(operation)
 
     //    override fun number(value: Number): Buffer<Double> = DoubleBuffer(size) { value.toDouble() }
 //

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferOps.kt

@@ -6,21 +6,35 @@
 package space.kscience.kmath.operations
 
 import space.kscience.kmath.linear.Point
+import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.BufferFactory
 import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.asBuffer
 
 import kotlin.math.*
 
 /**
  * [ExtendedFieldOps] over [DoubleBuffer].
  */
-public abstract class DoubleBufferOps : ExtendedFieldOps<Buffer<Double>>, Norm<Buffer<Double>, Double> {
+public abstract class DoubleBufferOps : BufferAlgebra<Double, DoubleField>, ExtendedFieldOps<Buffer<Double>>,
+    Norm<Buffer<Double>, Double> {
 
-    override fun Buffer<Double>.unaryMinus(): DoubleBuffer = if (this is DoubleBuffer) {
+    override val elementAlgebra: DoubleField get() = DoubleField
-        DoubleBuffer(size) { -array[it] }
+    override val bufferFactory: BufferFactory<Double> get() = ::DoubleBuffer
-    } else {
+
-        DoubleBuffer(size) { -get(it) }
+    override fun Buffer<Double>.map(block: DoubleField.(Double) -> Double): DoubleBuffer =
-    }
+        mapInline { DoubleField.block(it) }
+
+    @UnstableKMathAPI
+    override fun unaryOperationFunction(operation: String): (arg: Buffer<Double>) -> Buffer<Double> =
+        super<ExtendedFieldOps>.unaryOperationFunction(operation)
+
+    @UnstableKMathAPI
+    override fun binaryOperationFunction(operation: String): (left: Buffer<Double>, right: Buffer<Double>) -> Buffer<Double> =
+        super<ExtendedFieldOps>.binaryOperationFunction(operation)
+
+    override fun Buffer<Double>.unaryMinus(): DoubleBuffer = mapInline { -it }
 
     override fun add(left: Buffer<Double>, right: Buffer<Double>): DoubleBuffer {
         require(right.size == left.size) {
@@ -34,18 +48,18 @@ public abstract class DoubleBufferOps : ExtendedFieldOps<Buffer<Double>>, Norm<B
         } else DoubleBuffer(DoubleArray(left.size) { left[it] + right[it] })
     }
 
-    override fun Buffer<Double>.plus(other: Buffer<Double>): DoubleBuffer = add(this, other)
+    override fun Buffer<Double>.plus(arg: Buffer<Double>): DoubleBuffer = add(this, arg)
 
-    override fun Buffer<Double>.minus(other: Buffer<Double>): DoubleBuffer {
+    override fun Buffer<Double>.minus(arg: Buffer<Double>): DoubleBuffer {
-        require(other.size == this.size) {
+        require(arg.size == this.size) {
-            "The size of the first buffer ${this.size} should be the same as for second one: ${other.size} "
+            "The size of the first buffer ${this.size} should be the same as for second one: ${arg.size} "
         }
 
-        return if (this is DoubleBuffer && other is DoubleBuffer) {
+        return if (this is DoubleBuffer && arg is DoubleBuffer) {
             val aArray = this.array
-            val bArray = other.array
+            val bArray = arg.array
             DoubleBuffer(DoubleArray(this.size) { aArray[it] - bArray[it] })
-        } else DoubleBuffer(DoubleArray(this.size) { this[it] - other[it] })
+        } else DoubleBuffer(DoubleArray(this.size) { this[it] - arg[it] })
     }
 
     //
@@ -76,8 +90,7 @@ public abstract class DoubleBufferOps : ExtendedFieldOps<Buffer<Double>>, Norm<B
             val aArray = left.array
             val bArray = right.array
             DoubleBuffer(DoubleArray(left.size) { aArray[it] * bArray[it] })
-        } else
+        } else DoubleBuffer(DoubleArray(left.size) { left[it] * right[it] })
-            DoubleBuffer(DoubleArray(left.size) { left[it] * right[it] })
     }
 
     override fun divide(left: Buffer<Double>, right: Buffer<Double>): DoubleBuffer {
@@ -92,101 +105,46 @@ public abstract class DoubleBufferOps : ExtendedFieldOps<Buffer<Double>>, Norm<B
         } else DoubleBuffer(DoubleArray(left.size) { left[it] / right[it] })
     }
 
-    override fun sin(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun sin(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::sin)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { sin(array[it]) })
-    } else DoubleBuffer(DoubleArray(arg.size) { sin(arg[it]) })
 
-    override fun cos(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun cos(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::cos)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { cos(array[it]) })
-    } else DoubleBuffer(DoubleArray(arg.size) { cos(arg[it]) })
 
-    override fun tan(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun tan(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::tan)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { tan(array[it]) })
-    } else DoubleBuffer(DoubleArray(arg.size) { tan(arg[it]) })
 
-    override fun asin(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun asin(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::asin)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { asin(array[it]) })
-    } else
-        DoubleBuffer(DoubleArray(arg.size) { asin(arg[it]) })
 
-    override fun acos(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun acos(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::acos)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { acos(array[it]) })
-    } else
-        DoubleBuffer(DoubleArray(arg.size) { acos(arg[it]) })
 
-    override fun atan(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun atan(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::atan)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { atan(array[it]) })
-    } else
-        DoubleBuffer(DoubleArray(arg.size) { atan(arg[it]) })
 
-    override fun sinh(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun sinh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::sinh)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { sinh(array[it]) })
-    } else
-        DoubleBuffer(DoubleArray(arg.size) { sinh(arg[it]) })
 
-    override fun cosh(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun cosh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::cosh)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { cosh(array[it]) })
-    } else
-        DoubleBuffer(DoubleArray(arg.size) { cosh(arg[it]) })
 
-    override fun tanh(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun tanh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::tanh)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { tanh(array[it]) })
-    } else
-        DoubleBuffer(DoubleArray(arg.size) { tanh(arg[it]) })
 
-    override fun asinh(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun asinh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::asinh)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { asinh(array[it]) })
-    } else
-        DoubleBuffer(DoubleArray(arg.size) { asinh(arg[it]) })
 
-    override fun acosh(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun acosh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::acosh)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { acosh(array[it]) })
-    } else
-        DoubleBuffer(DoubleArray(arg.size) { acosh(arg[it]) })
 
-    override fun atanh(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun atanh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::atanh)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { atanh(array[it]) })
-    } else
-        DoubleBuffer(DoubleArray(arg.size) { atanh(arg[it]) })
 
-    override fun power(arg: Buffer<Double>, pow: Number): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun exp(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::exp)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { array[it].pow(pow.toDouble()) })
-    } else
-        DoubleBuffer(DoubleArray(arg.size) { arg[it].pow(pow.toDouble()) })
 
-    override fun exp(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
+    override fun ln(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::ln)
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { exp(array[it]) })
-    } else DoubleBuffer(DoubleArray(arg.size) { exp(arg[it]) })
-
-    override fun ln(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
-        val array = arg.array
-        DoubleBuffer(DoubleArray(arg.size) { ln(array[it]) })
-    } else {
-        DoubleBuffer(DoubleArray(arg.size) { ln(arg[it]) })
-    }
 
     override fun norm(arg: Buffer<Double>): Double = DoubleL2Norm.norm(arg)
 
-    override fun scale(a: Buffer<Double>, value: Double): DoubleBuffer = if (a is DoubleBuffer) {
+    override fun scale(a: Buffer<Double>, value: Double): DoubleBuffer = a.mapInline { it * value }
-        val aArray = a.array
-        DoubleBuffer(DoubleArray(a.size) { aArray[it] * value })
-    } else DoubleBuffer(DoubleArray(a.size) { a[it] * value })
 
-    public companion object : DoubleBufferOps()
+    public companion object : DoubleBufferOps() {
+        public inline fun Buffer<Double>.mapInline(block: (Double) -> Double): DoubleBuffer =
+            if (this is DoubleBuffer) {
+                DoubleArray(size) { block(array[it]) }.asBuffer()
+            } else {
+                DoubleArray(size) { block(get(it)) }.asBuffer()
+            }
+    }
 }
 
 public object DoubleL2Norm : Norm<Point<Double>, Double> {

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/OptionalOperations.kt

@@ -74,14 +74,21 @@ public interface TrigonometricOperations<T> : Algebra<T> {
     }
 }
 
+/**
+ * Check if number is an integer from platform point of view
+ */
+public expect fun Number.isInteger(): Boolean
+
 /**
  * A context extension to include power operations based on exponentiation.
  *
  * @param T the type of element of this structure.
  */
-public interface PowerOperations<T> : Algebra<T> {
+public interface PowerOperations<T> : FieldOps<T> {
+
     /**
-     * Raises [arg] to the power [pow].
+     * Raises [arg] to a power if possible (negative number could not be raised to a fractional power).
+     * Throws [IllegalArgumentException] if not possible.
      */
     public fun power(arg: T, pow: Number): T
 
@@ -108,6 +115,7 @@ public interface PowerOperations<T> : Algebra<T> {
     }
 }
 
+
 /**
  * A container for operations related to `exp` and `ln` functions.
  *

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/algebraExtensions.kt

@@ -96,46 +96,3 @@ public fun <T, S> Iterable<T>.averageWith(space: S): T where S : Ring<T>, S : Sc
 public fun <T, S> Sequence<T>.averageWith(space: S): T where S : Ring<T>, S : ScaleOperations<T> =
     space.average(this)
 
-/**
- * Raises [arg] to the non-negative integer power [exponent].
- *
- * Special case: 0 ^ 0 is 1.
- *
- * @receiver the algebra to provide multiplication.
- * @param arg the base.
- * @param exponent the exponent.
- * @return the base raised to the power.
- * @author Evgeniy Zhelenskiy
- */
-public fun <T> Ring<T>.power(arg: T, exponent: UInt): T = when {
-    arg == zero && exponent > 0U -> zero
-    arg == one -> arg
-    arg == -one -> powWithoutOptimization(arg, exponent % 2U)
-    else -> powWithoutOptimization(arg, exponent)
-}
-
-private fun <T> Ring<T>.powWithoutOptimization(base: T, exponent: UInt): T = when (exponent) {
-    0U -> one
-    1U -> base
-    else -> {
-        val pre = powWithoutOptimization(base, exponent shr 1).let { it * it }
-        if (exponent and 1U == 0U) pre else pre * base
-    }
-}
-
-
-/**
- * Raises [arg] to the integer power [exponent].
- *
- * Special case: 0 ^ 0 is 1.
- *
- * @receiver the algebra to provide multiplication and division.
- * @param arg the base.
- * @param exponent the exponent.
- * @return the base raised to the power.
- * @author Iaroslav Postovalov, Evgeniy Zhelenskiy
- */
-public fun <T> Field<T>.power(arg: T, exponent: Int): T = when {
-    exponent < 0 -> one / (this as Ring<T>).power(arg, if (exponent == Int.MIN_VALUE) Int.MAX_VALUE.toUInt().inc() else (-exponent).toUInt())
-    else -> (this as Ring<T>).power(arg, exponent.toUInt())
-}

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt

@@ -13,9 +13,8 @@ import kotlin.math.pow as kpow
 public interface ExtendedFieldOps<T> :
     FieldOps<T>,
     TrigonometricOperations<T>,
-    PowerOperations<T>,
     ExponentialOperations<T>,
-    ScaleOperations<T>  {
+    ScaleOperations<T> {
     override fun tan(arg: T): T = sin(arg) / cos(arg)
     override fun tanh(arg: T): T = sinh(arg) / cosh(arg)
 
@@ -26,7 +25,6 @@ public interface ExtendedFieldOps<T> :
         TrigonometricOperations.ACOS_OPERATION -> ::acos
         TrigonometricOperations.ASIN_OPERATION -> ::asin
         TrigonometricOperations.ATAN_OPERATION -> ::atan
-        PowerOperations.SQRT_OPERATION -> ::sqrt
         ExponentialOperations.EXP_OPERATION -> ::exp
         ExponentialOperations.LN_OPERATION -> ::ln
         ExponentialOperations.COSH_OPERATION -> ::cosh
@@ -42,7 +40,7 @@ public interface ExtendedFieldOps<T> :
 /**
  * Advanced Number-like field that implements basic operations.
  */
-public interface ExtendedField<T> : ExtendedFieldOps<T>, Field<T>, NumericAlgebra<T>{
+public interface ExtendedField<T> : ExtendedFieldOps<T>, Field<T>, PowerOperations<T>, NumericAlgebra<T> {
     override fun sinh(arg: T): T = (exp(arg) - exp(-arg)) / 2.0
     override fun cosh(arg: T): T = (exp(arg) + exp(-arg)) / 2.0
     override fun tanh(arg: T): T = (exp(arg) - exp(-arg)) / (exp(-arg) + exp(arg))
@@ -50,6 +48,11 @@ public interface ExtendedField<T> : ExtendedFieldOps<T>, Field<T>, NumericAlgebr
     override fun acosh(arg: T): T = ln(arg + sqrt((arg - one) * (arg + one)))
     override fun atanh(arg: T): T = (ln(arg + one) - ln(one - arg)) / 2.0
 
+    override fun unaryOperationFunction(operation: String): (arg: T) -> T {
+        return if (operation == PowerOperations.SQRT_OPERATION) ::sqrt
+        else super<ExtendedFieldOps>.unaryOperationFunction(operation)
+    }
+
     override fun rightSideNumberOperationFunction(operation: String): (left: T, right: Number) -> T =
         when (operation) {
             PowerOperations.POW_OPERATION -> ::power
@@ -69,7 +72,7 @@ public object DoubleField : ExtendedField<Double>, Norm<Double, Double>, ScaleOp
 
     override fun binaryOperationFunction(operation: String): (left: Double, right: Double) -> Double =
         when (operation) {
-            PowerOperations.POW_OPERATION -> ::power
+            PowerOperations.POW_OPERATION -> { l, r -> l.kpow(r) }
             else -> super<ExtendedField>.binaryOperationFunction(operation)
         }
 
@@ -94,18 +97,23 @@ public object DoubleField : ExtendedField<Double>, Norm<Double, Double>, ScaleOp
     override inline fun acosh(arg: Double): Double = kotlin.math.acosh(arg)
     override inline fun atanh(arg: Double): Double = kotlin.math.atanh(arg)
 
-    override inline fun sqrt(arg: Double): Double = kotlin.math.sqrt(arg)
+    override fun sqrt(arg: Double): Double = kotlin.math.sqrt(arg)
-    override inline fun power(arg: Double, pow: Number): Double = arg.kpow(pow.toDouble())
+    override fun power(arg: Double, pow: Number): Double = when {
+        pow.isInteger() -> arg.kpow(pow.toInt())
+        arg < 0 -> throw IllegalArgumentException("Can't raise negative $arg to a fractional power $pow")
+        else -> arg.kpow(pow.toDouble())
+    }
+
     override inline fun exp(arg: Double): Double = kotlin.math.exp(arg)
     override inline fun ln(arg: Double): Double = kotlin.math.ln(arg)
 
     override inline fun norm(arg: Double): Double = abs(arg)
 
     override inline fun Double.unaryMinus(): Double = -this
-    override inline fun Double.plus(other: Double): Double = this + other
+    override inline fun Double.plus(arg: Double): Double = this + arg
-    override inline fun Double.minus(other: Double): Double = this - other
+    override inline fun Double.minus(arg: Double): Double = this - arg
-    override inline fun Double.times(other: Double): Double = this * other
+    override inline fun Double.times(arg: Double): Double = this * arg
-    override inline fun Double.div(other: Double): Double = this / other
+    override inline fun Double.div(arg: Double): Double = this / arg
 }
 
 public val Double.Companion.algebra: DoubleField get() = DoubleField
@@ -122,7 +130,7 @@ public object FloatField : ExtendedField<Float>, Norm<Float, Float> {
 
     override fun binaryOperationFunction(operation: String): (left: Float, right: Float) -> Float =
         when (operation) {
-            PowerOperations.POW_OPERATION -> ::power
+            PowerOperations.POW_OPERATION -> { l, r -> l.kpow(r) }
             else -> super.binaryOperationFunction(operation)
         }
 
@@ -149,16 +157,17 @@ public object FloatField : ExtendedField<Float>, Norm<Float, Float> {
 
     override inline fun sqrt(arg: Float): Float = kotlin.math.sqrt(arg)
     override inline fun power(arg: Float, pow: Number): Float = arg.kpow(pow.toFloat())
+
     override inline fun exp(arg: Float): Float = kotlin.math.exp(arg)
     override inline fun ln(arg: Float): Float = kotlin.math.ln(arg)
 
     override inline fun norm(arg: Float): Float = abs(arg)
 
     override inline fun Float.unaryMinus(): Float = -this
-    override inline fun Float.plus(other: Float): Float = this + other
+    override inline fun Float.plus(arg: Float): Float = this + arg
-    override inline fun Float.minus(other: Float): Float = this - other
+    override inline fun Float.minus(arg: Float): Float = this - arg
-    override inline fun Float.times(other: Float): Float = this * other
+    override inline fun Float.times(arg: Float): Float = this * arg
-    override inline fun Float.div(other: Float): Float = this / other
+    override inline fun Float.div(arg: Float): Float = this / arg
 }
 
 public val Float.Companion.algebra: FloatField get() = FloatField
@@ -180,9 +189,9 @@ public object IntRing : Ring<Int>, Norm<Int, Int>, NumericAlgebra<Int> {
     override inline fun norm(arg: Int): Int = abs(arg)
 
     override inline fun Int.unaryMinus(): Int = -this
-    override inline fun Int.plus(other: Int): Int = this + other
+    override inline fun Int.plus(arg: Int): Int = this + arg
-    override inline fun Int.minus(other: Int): Int = this - other
+    override inline fun Int.minus(arg: Int): Int = this - arg
-    override inline fun Int.times(other: Int): Int = this * other
+    override inline fun Int.times(arg: Int): Int = this * arg
 }
 
 public val Int.Companion.algebra: IntRing get() = IntRing
@@ -204,9 +213,9 @@ public object ShortRing : Ring<Short>, Norm<Short, Short>, NumericAlgebra<Short>
     override fun norm(arg: Short): Short = if (arg > 0) arg else (-arg).toShort()
 
     override inline fun Short.unaryMinus(): Short = (-this).toShort()
-    override inline fun Short.plus(other: Short): Short = (this + other).toShort()
+    override inline fun Short.plus(arg: Short): Short = (this + arg).toShort()
-    override inline fun Short.minus(other: Short): Short = (this - other).toShort()
+    override inline fun Short.minus(arg: Short): Short = (this - arg).toShort()
-    override inline fun Short.times(other: Short): Short = (this * other).toShort()
+    override inline fun Short.times(arg: Short): Short = (this * arg).toShort()
 }
 
 public val Short.Companion.algebra: ShortRing get() = ShortRing
@@ -230,7 +239,7 @@ public object ByteRing : Ring<Byte>, Norm<Byte, Byte>, NumericAlgebra<Byte> {
     override inline fun Byte.unaryMinus(): Byte = (-this).toByte()
     override inline fun Byte.plus(arg: Byte): Byte = (this + arg).toByte()
     override inline fun Byte.minus(arg: Byte): Byte = (this - arg).toByte()
-    override inline fun Byte.times(other: Byte): Byte = (this * other).toByte()
+    override inline fun Byte.times(arg: Byte): Byte = (this * arg).toByte()
 }
 
 public val Byte.Companion.algebra: ByteRing get() = ByteRing
@@ -252,9 +261,9 @@ public object LongRing : Ring<Long>, Norm<Long, Long>, NumericAlgebra<Long> {
     override fun norm(arg: Long): Long = abs(arg)
 
     override inline fun Long.unaryMinus(): Long = (-this)
-    override inline fun Long.plus(other: Long): Long = (this + other)
+    override inline fun Long.plus(arg: Long): Long = (this + arg)
-    override inline fun Long.minus(other: Long): Long = (this - other)
+    override inline fun Long.minus(arg: Long): Long = (this - arg)
-    override inline fun Long.times(other: Long): Long = (this * other)
+    override inline fun Long.times(arg: Long): Long = (this * arg)
 }
 
 public val Long.Companion.algebra: LongRing get() = LongRing

kmath-core/src/jsMain/kotlin/space/kscience/kmath/operations/isInteger.kt

@@ -0,0 +1,6 @@
+package space.kscience.kmath.operations
+
+/**
+ * Check if number is an integer
+ */
+public actual fun Number.isInteger(): Boolean = js("Number").isInteger(this) as Boolean

kmath-core/src/jvmMain/kotlin/space/kscience/kmath/operations/isInteger.kt

@@ -0,0 +1,6 @@
+package space.kscience.kmath.operations
+
+/**
+ * Check if number is an integer
+ */
+public actual fun Number.isInteger(): Boolean = (this is Int) || (this is Long) || (this is Short) || (this.toDouble() % 1 == 0.0)

kmath-core/src/nativeMain/kotlin/space/kscience/kmath/operations/isInteger.kt

@@ -0,0 +1,6 @@
+package space.kscience.kmath.operations
+
+/**
+ * Check if number is an integer
+ */
+public actual fun Number.isInteger(): Boolean = (this is Int) || (this is Long) || (this is Short) || (this.toDouble() % 1 == 0.0)

kmath-jafama/src/main/kotlin/space/kscience/kmath/jafama/KMathJafama.kt

@@ -59,8 +59,8 @@ public object JafamaDoubleField : ExtendedField<Double>, Norm<Double, Double>, S
     override inline fun Double.unaryMinus(): Double = -this
     override inline fun Double.plus(arg: Double): Double = this + arg
     override inline fun Double.minus(arg: Double): Double = this - arg
-    override inline fun Double.times(other: Double): Double = this * other
+    override inline fun Double.times(arg: Double): Double = this * arg
-    override inline fun Double.div(other: Double): Double = this / other
+    override inline fun Double.div(arg: Double): Double = this / arg
 }
 
 /**
@@ -108,8 +108,8 @@ public object StrictJafamaDoubleField : ExtendedField<Double>, Norm<Double, Doub
     override inline fun norm(arg: Double): Double = StrictFastMath.abs(arg)
 
     override inline fun Double.unaryMinus(): Double = -this
-    override inline fun Double.plus(other: Double): Double = this + other
+    override inline fun Double.plus(arg: Double): Double = this + arg
-    override inline fun Double.minus(other: Double): Double = this - other
+    override inline fun Double.minus(arg: Double): Double = this - arg
-    override inline fun Double.times(other: Double): Double = this * other
+    override inline fun Double.times(arg: Double): Double = this * arg
-    override inline fun Double.div(other: Double): Double = this / other
+    override inline fun Double.div(arg: Double): Double = this / arg
 }

kmath-kotlingrad/src/test/kotlin/space/kscience/kmath/kotlingrad/AdaptingTests.kt

@@ -62,6 +62,6 @@ internal class AdaptingTests {
             .parseMath()
             .compileToExpression(DoubleField)
 
-        assertEquals(actualDerivative(x to 0.1), expectedDerivative(x to 0.1))
+        assertEquals(actualDerivative(x to -0.1), expectedDerivative(x to -0.1))
     }
 }

kmath-multik/src/main/kotlin/space/kscience/kmath/multik/MultikDoubleAlgebra.kt

@@ -0,0 +1,49 @@
+package space.kscience.kmath.multik
+
+import org.jetbrains.kotlinx.multik.ndarray.data.DataType
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.ExponentialOperations
+import space.kscience.kmath.operations.TrigonometricOperations
+
+public object MultikDoubleAlgebra : MultikDivisionTensorAlgebra<Double, DoubleField>(),
+    TrigonometricOperations<StructureND<Double>>, ExponentialOperations<StructureND<Double>> {
+    override val elementAlgebra: DoubleField get() = DoubleField
+    override val type: DataType get() = DataType.DoubleDataType
+
+    override fun sin(arg: StructureND<Double>): MultikTensor<Double> = multikMath.mathEx.sin(arg.asMultik().array).wrap()
+
+    override fun cos(arg: StructureND<Double>): MultikTensor<Double> = multikMath.mathEx.cos(arg.asMultik().array).wrap()
+
+    override fun tan(arg: StructureND<Double>): MultikTensor<Double> = sin(arg) / cos(arg)
+
+    override fun asin(arg: StructureND<Double>): MultikTensor<Double>  = arg.map { asin(it) }
+
+    override fun acos(arg: StructureND<Double>): MultikTensor<Double> = arg.map { acos(it) }
+
+    override fun atan(arg: StructureND<Double>): MultikTensor<Double> = arg.map { atan(it) }
+
+    override fun exp(arg: StructureND<Double>): MultikTensor<Double> = multikMath.mathEx.exp(arg.asMultik().array).wrap()
+
+    override fun ln(arg: StructureND<Double>): MultikTensor<Double> = multikMath.mathEx.log(arg.asMultik().array).wrap()
+
+    override fun sinh(arg: StructureND<Double>): MultikTensor<Double> = (exp(arg) - exp(-arg)) / 2.0
+
+    override fun cosh(arg: StructureND<Double>): MultikTensor<Double> = (exp(arg) + exp(-arg)) / 2.0
+
+    override fun tanh(arg: StructureND<Double>): MultikTensor<Double> {
+        val expPlus = exp(arg)
+        val expMinus = exp(-arg)
+        return (expPlus - expMinus) / (expPlus + expMinus)
+    }
+
+    override fun asinh(arg: StructureND<Double>): MultikTensor<Double> = arg.map { asinh(it) }
+
+    override fun acosh(arg: StructureND<Double>): MultikTensor<Double> = arg.map { acosh(it) }
+
+    override fun atanh(arg: StructureND<Double>): MultikTensor<Double> = arg.map { atanh(it) }
+}
+
+public val Double.Companion.multikAlgebra: MultikTensorAlgebra<Double, DoubleField> get() = MultikDoubleAlgebra
+public val DoubleField.multikAlgebra: MultikTensorAlgebra<Double, DoubleField> get() = MultikDoubleAlgebra
+

kmath-multik/src/main/kotlin/space/kscience/kmath/multik/MultikTensorAlgebra.kt

@@ -7,11 +7,9 @@
 
 package space.kscience.kmath.multik
 
-import org.jetbrains.kotlinx.multik.api.Multik
+import org.jetbrains.kotlinx.multik.api.*
-import org.jetbrains.kotlinx.multik.api.linalg.dot
+import org.jetbrains.kotlinx.multik.api.linalg.LinAlg
-import org.jetbrains.kotlinx.multik.api.mk
+import org.jetbrains.kotlinx.multik.api.math.Math
-import org.jetbrains.kotlinx.multik.api.ndarrayOf
-import org.jetbrains.kotlinx.multik.api.zeros
 import org.jetbrains.kotlinx.multik.ndarray.data.*
 import org.jetbrains.kotlinx.multik.ndarray.operations.*
 import space.kscience.kmath.misc.PerformancePitfall
@@ -52,10 +50,16 @@ private fun <T, D : Dimension> MultiArray<T, D>.asD2Array(): D2Array<T> {
     else throw ClassCastException("Cannot cast MultiArray to NDArray.")
 }
 
-public abstract class MultikTensorAlgebra<T, A : Ring<T>> : TensorAlgebra<T, A> where T : Number, T : Comparable<T> {
+public abstract class MultikTensorAlgebra<T, A : Ring<T>> : TensorAlgebra<T, A>
+        where T : Number, T : Comparable<T> {
 
     public abstract val type: DataType
 
+    protected val multikMath: Math = mk.math
+    protected val multikLinAl: LinAlg = mk.linalg
+    protected val multikStat: Statistics = mk.stat
+
+
     override fun structureND(shape: Shape, initializer: A.(IntArray) -> T): MultikTensor<T> {
         val strides = DefaultStrides(shape)
         val memoryView = initMemoryView<T>(strides.linearSize, type)
@@ -65,6 +69,7 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>> : TensorAlgebra<T, A>
         return MultikTensor(NDArray(memoryView, shape = shape, dim = DN(shape.size)))
     }
 
+    @OptIn(PerformancePitfall::class)
     override fun StructureND<T>.map(transform: A.(T) -> T): MultikTensor<T> = if (this is MultikTensor) {
         val data = initMemoryView<T>(array.size, type)
         var count = 0
@@ -76,6 +81,7 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>> : TensorAlgebra<T, A>
         }
     }
 
+    @OptIn(PerformancePitfall::class)
     override fun StructureND<T>.mapIndexed(transform: A.(index: IntArray, T) -> T): MultikTensor<T> =
         if (this is MultikTensor) {
             val array = asMultik().array
@@ -96,6 +102,7 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>> : TensorAlgebra<T, A>
             }
         }
 
+    @OptIn(PerformancePitfall::class)
     override fun zip(left: StructureND<T>, right: StructureND<T>, transform: A.(T, T) -> T): MultikTensor<T> {
         require(left.shape.contentEquals(right.shape)) { "ND array shape mismatch" } //TODO replace by ShapeMismatchException
         val leftArray = left.asMultik().array
@@ -208,9 +215,9 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>> : TensorAlgebra<T, A>
     override fun StructureND<T>.unaryMinus(): MultikTensor<T> =
         asMultik().array.unaryMinus().wrap()
 
-    override fun StructureND<T>.get(i: Int): MultikTensor<T> = asMultik().array.mutableView(i).wrap()
+    override fun Tensor<T>.get(i: Int): MultikTensor<T> = asMultik().array.mutableView(i).wrap()
 
-    override fun StructureND<T>.transpose(i: Int, j: Int): MultikTensor<T> = asMultik().array.transpose(i, j).wrap()
+    override fun Tensor<T>.transpose(i: Int, j: Int): MultikTensor<T> = asMultik().array.transpose(i, j).wrap()
 
     override fun Tensor<T>.view(shape: IntArray): MultikTensor<T> {
         require(shape.all { it > 0 })
@@ -236,12 +243,12 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>> : TensorAlgebra<T, A>
     override fun StructureND<T>.dot(other: StructureND<T>): MultikTensor<T> =
         if (this.shape.size == 1 && other.shape.size == 1) {
             Multik.ndarrayOf(
-                asMultik().array.asD1Array() dot other.asMultik().array.asD1Array()
+                multikLinAl.linAlgEx.dotVV(asMultik().array.asD1Array(), other.asMultik().array.asD1Array())
-            ).asDNArray().wrap()
+            ).wrap()
         } else if (this.shape.size == 2 && other.shape.size == 2) {
-            (asMultik().array.asD2Array() dot other.asMultik().array.asD2Array()).asDNArray().wrap()
+            multikLinAl.linAlgEx.dotMM(asMultik().array.asD2Array(), other.asMultik().array.asD2Array()).wrap()
         } else if (this.shape.size == 2 && other.shape.size == 1) {
-            (asMultik().array.asD2Array() dot other.asMultik().array.asD1Array()).asDNArray().wrap()
+            multikLinAl.linAlgEx.dotMV(asMultik().array.asD2Array(), other.asMultik().array.asD1Array()).wrap()
         } else {
             TODO("Not implemented for broadcasting")
         }
@@ -270,7 +277,7 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>> : TensorAlgebra<T, A>
         TODO("Not yet implemented")
     }
 
-    override fun StructureND<T>.argMax(dim: Int, keepDim: Boolean): Tensor<T> {
+    override fun StructureND<T>.argMax(dim: Int, keepDim: Boolean): Tensor<Int> {
         TODO("Not yet implemented")
     }
 }
@@ -294,23 +301,15 @@ public abstract class MultikDivisionTensorAlgebra<T, A : Field<T>>
         }
     }
 
-    override fun Tensor<T>.divAssign(other: StructureND<T>) {
+    override fun Tensor<T>.divAssign(arg: StructureND<T>) {
         if (this is MultikTensor) {
-            array.divAssign(other.asMultik().array)
+            array.divAssign(arg.asMultik().array)
         } else {
-            mapInPlace { index, t -> elementAlgebra.divide(t, other[index]) }
+            mapInPlace { index, t -> elementAlgebra.divide(t, arg[index]) }
         }
     }
 }
 
-public object MultikDoubleAlgebra : MultikDivisionTensorAlgebra<Double, DoubleField>() {
-    override val elementAlgebra: DoubleField get() = DoubleField
-    override val type: DataType get() = DataType.DoubleDataType
-}
-
-public val Double.Companion.multikAlgebra: MultikTensorAlgebra<Double, DoubleField> get() = MultikDoubleAlgebra
-public val DoubleField.multikAlgebra: MultikTensorAlgebra<Double, DoubleField> get() = MultikDoubleAlgebra
-
 public object MultikFloatAlgebra : MultikDivisionTensorAlgebra<Float, FloatField>() {
     override val elementAlgebra: FloatField get() = FloatField
     override val type: DataType get() = DataType.FloatDataType

kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jArrayAlgebra.kt

@@ -155,7 +155,7 @@ public sealed interface Nd4jArrayField<T, out F : Field<T>> : FieldOpsND<T, F>,
  * Represents intersection of [ExtendedField] and [Field] over [Nd4jArrayStructure].
  */
 public sealed interface Nd4jArrayExtendedFieldOps<T, out F : ExtendedField<T>> :
-    ExtendedFieldOps<StructureND<T>>, Nd4jArrayField<T, F> {
+    ExtendedFieldOps<StructureND<T>>, Nd4jArrayField<T, F>, PowerOperations<StructureND<T>> {
 
     override fun sin(arg: StructureND<T>): StructureND<T> = Transforms.sin(arg.ndArray).wrap()
     override fun cos(arg: StructureND<T>): StructureND<T> = Transforms.cos(arg.ndArray).wrap()

kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jTensorAlgebra.kt

@@ -92,8 +92,8 @@ public sealed interface Nd4jTensorAlgebra<T : Number, A : Field<T>> : AnalyticTe
     }
 
     override fun StructureND<T>.unaryMinus(): Nd4jArrayStructure<T> = ndArray.neg().wrap()
-    override fun StructureND<T>.get(i: Int): Nd4jArrayStructure<T> = ndArray.slice(i.toLong()).wrap()
+    override fun Tensor<T>.get(i: Int): Nd4jArrayStructure<T> = ndArray.slice(i.toLong()).wrap()
-    override fun StructureND<T>.transpose(i: Int, j: Int): Nd4jArrayStructure<T> = ndArray.swapAxes(i, j).wrap()
+    override fun Tensor<T>.transpose(i: Int, j: Int): Nd4jArrayStructure<T> = ndArray.swapAxes(i, j).wrap()
     override fun StructureND<T>.dot(other: StructureND<T>): Nd4jArrayStructure<T> = ndArray.mmul(other.ndArray).wrap()
 
     override fun StructureND<T>.min(dim: Int, keepDim: Boolean): Nd4jArrayStructure<T> =
@@ -108,8 +108,8 @@ public sealed interface Nd4jTensorAlgebra<T : Number, A : Field<T>> : AnalyticTe
     override fun Tensor<T>.view(shape: IntArray): Nd4jArrayStructure<T> = ndArray.reshape(shape).wrap()
     override fun Tensor<T>.viewAs(other: StructureND<T>): Nd4jArrayStructure<T> = view(other.shape)
 
-    override fun StructureND<T>.argMax(dim: Int, keepDim: Boolean): Nd4jArrayStructure<T> =
+    override fun StructureND<T>.argMax(dim: Int, keepDim: Boolean): Tensor<Int> =
-        ndBase.get().argmax(ndArray, keepDim, dim).wrap()
+        ndBase.get().argmax(ndArray, keepDim, dim).asIntStructure()
 
     override fun StructureND<T>.mean(dim: Int, keepDim: Boolean): Nd4jArrayStructure<T> =
         ndArray.mean(keepDim, dim).wrap()
@@ -148,8 +148,8 @@ public sealed interface Nd4jTensorAlgebra<T : Number, A : Field<T>> : AnalyticTe
         ndArray.divi(value)
     }
 
-    override fun Tensor<T>.divAssign(other: StructureND<T>) {
+    override fun Tensor<T>.divAssign(arg: StructureND<T>) {
-        ndArray.divi(other.ndArray)
+        ndArray.divi(arg.ndArray)
     }
 
     override fun StructureND<T>.variance(dim: Int, keepDim: Boolean): Nd4jArrayStructure<T> =

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/AnalyticTensorAlgebra.kt

@@ -6,6 +6,7 @@
 package space.kscience.kmath.tensors.api
 
 import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.operations.ExtendedFieldOps
 import space.kscience.kmath.operations.Field
 
 
@@ -14,7 +15,8 @@ import space.kscience.kmath.operations.Field
  *
  * @param T the type of items closed under analytic functions in the tensors.
  */
-public interface AnalyticTensorAlgebra<T, A : Field<T>> : TensorPartialDivisionAlgebra<T, A> {
+public interface AnalyticTensorAlgebra<T, A : Field<T>> :
+    TensorPartialDivisionAlgebra<T, A>, ExtendedFieldOps<StructureND<T>> {
 
     /**
      * @return the mean of all elements in the input tensor.
@@ -121,4 +123,27 @@ public interface AnalyticTensorAlgebra<T, A : Field<T>> : TensorPartialDivisionA
     //For information: https://pytorch.org/docs/stable/generated/torch.floor.html#torch.floor
     public fun StructureND<T>.floor(): Tensor<T>
 
+    override fun sin(arg: StructureND<T>): StructureND<T> = arg.sin()
+
+    override fun cos(arg: StructureND<T>): StructureND<T> = arg.cos()
+
+    override fun asin(arg: StructureND<T>): StructureND<T> = arg.asin()
+
+    override fun acos(arg: StructureND<T>): StructureND<T> = arg.acos()
+
+    override fun atan(arg: StructureND<T>): StructureND<T> = arg.atan()
+
+    override fun exp(arg: StructureND<T>): StructureND<T> = arg.exp()
+
+    override fun ln(arg: StructureND<T>): StructureND<T> = arg.ln()
+
+    override fun sinh(arg: StructureND<T>): StructureND<T> = arg.sinh()
+
+    override fun cosh(arg: StructureND<T>): StructureND<T> = arg.cosh()
+
+    override fun asinh(arg: StructureND<T>): StructureND<T> = arg.asinh()
+
+    override fun acosh(arg: StructureND<T>): StructureND<T> = arg.acosh()
+
+    override fun atanh(arg: StructureND<T>): StructureND<T> = arg.atanh()
 }

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorAlgebra.kt

@@ -166,7 +166,7 @@ public interface TensorAlgebra<T, A : Ring<T>> : RingOpsND<T, A> {
      * @param i index of the extractable tensor
      * @return subtensor of the original tensor with index [i]
      */
-    public operator fun StructureND<T>.get(i: Int): Tensor<T>
+    public operator fun Tensor<T>.get(i: Int): Tensor<T>
 
     /**
      * Returns a tensor that is a transposed version of this tensor. The given dimensions [i] and [j] are swapped.
@@ -176,7 +176,7 @@ public interface TensorAlgebra<T, A : Ring<T>> : RingOpsND<T, A> {
      * @param j the second dimension to be transposed
      * @return transposed tensor
      */
-    public fun StructureND<T>.transpose(i: Int = -2, j: Int = -1): Tensor<T>
+    public fun Tensor<T>.transpose(i: Int = -2, j: Int = -1): Tensor<T>
 
     /**
      * Returns a new tensor with the same data as the self tensor but of a different shape.
@@ -324,7 +324,7 @@ public interface TensorAlgebra<T, A : Ring<T>> : RingOpsND<T, A> {
      * @param keepDim whether the output tensor has [dim] retained or not.
      * @return the index of maximum value of each row of the input tensor in the given dimension [dim].
      */
-    public fun StructureND<T>.argMax(dim: Int, keepDim: Boolean): Tensor<T>
+    public fun StructureND<T>.argMax(dim: Int, keepDim: Boolean): Tensor<Int>
 
     override fun add(left: StructureND<T>, right: StructureND<T>): Tensor<T> = left + right
 

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/api/TensorPartialDivisionAlgebra.kt

@@ -53,9 +53,9 @@ public interface TensorPartialDivisionAlgebra<T, A : Field<T>> : TensorAlgebra<T
     public operator fun Tensor<T>.divAssign(value: T)
 
     /**
-     * Each element of this tensor is divided by each element of the [other] tensor.
+     * Each element of this tensor is divided by each element of the [arg] tensor.
      *
-     * @param other tensor to be divided by.
+     * @param arg tensor to be divided by.
      */
-    public operator fun Tensor<T>.divAssign(other: StructureND<T>)
+    public operator fun Tensor<T>.divAssign(arg: StructureND<T>)
 }

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.tensors.core
 
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.tensors.api.Tensor
@@ -17,6 +18,8 @@ import space.kscience.kmath.tensors.core.internal.tensor
  * Basic linear algebra operations implemented with broadcasting.
  * For more information: https://pytorch.org/docs/stable/notes/broadcasting.html
  */
+
+@PerformancePitfall
 public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
 
     override fun StructureND<Double>.plus(arg: StructureND<Double>): DoubleTensor {
@@ -74,8 +77,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
         }
     }
 
-    override fun StructureND<Double>.div(other: StructureND<Double>): DoubleTensor {
+    override fun StructureND<Double>.div(arg: StructureND<Double>): DoubleTensor {
-        val broadcast = broadcastTensors(tensor, other.tensor)
+        val broadcast = broadcastTensors(tensor, arg.tensor)
         val newThis = broadcast[0]
         val newOther = broadcast[1]
         val resBuffer = DoubleArray(newThis.indices.linearSize) { i ->
@@ -85,8 +88,8 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
         return DoubleTensor(newThis.shape, resBuffer)
     }
 
-    override fun Tensor<Double>.divAssign(other: StructureND<Double>) {
+    override fun Tensor<Double>.divAssign(arg: StructureND<Double>) {
-        val newOther = broadcastTo(other.tensor, tensor.shape)
+        val newOther = broadcastTo(arg.tensor, tensor.shape)
         for (i in 0 until tensor.indices.linearSize) {
             tensor.mutableBuffer.array()[tensor.bufferStart + i] /=
                 newOther.mutableBuffer.array()[tensor.bufferStart + i]
@@ -99,5 +102,6 @@ public object BroadcastDoubleTensorAlgebra : DoubleTensorAlgebra() {
  * Compute a value using broadcast double tensor algebra
  */
 @UnstableKMathAPI
+@PerformancePitfall
 public fun <R> DoubleTensorAlgebra.withBroadcast(block: BroadcastDoubleTensorAlgebra.() -> R): R =
     BroadcastDoubleTensorAlgebra.block()

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt

@@ -9,7 +9,6 @@ import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.nd.Strides
 import space.kscience.kmath.structures.MutableBuffer
 import space.kscience.kmath.tensors.api.Tensor
-import space.kscience.kmath.tensors.core.internal.TensorLinearStructure
 
 /**
  * Represents [Tensor] over a [MutableBuffer] intended to be used through [DoubleTensor] and [IntTensor]

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt

@@ -3,13 +3,18 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
  */
 
+
+@file:OptIn(PerformancePitfall::class)
+
 package space.kscience.kmath.tensors.core
 
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.nd.MutableStructure2D
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.nd.as1D
 import space.kscience.kmath.nd.as2D
 import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.structures.MutableBuffer
 import space.kscience.kmath.structures.indices
 import space.kscience.kmath.tensors.api.AnalyticTensorAlgebra
 import space.kscience.kmath.tensors.api.LinearOpsTensorAlgebra
@@ -38,6 +43,7 @@ public open class DoubleTensorAlgebra :
      * @param transform the function to be applied to each element of the tensor.
      * @return the resulting tensor after applying the function.
      */
+    @PerformancePitfall
     @Suppress("OVERRIDE_BY_INLINE")
     final override inline fun StructureND<Double>.map(transform: DoubleField.(Double) -> Double): DoubleTensor {
         val tensor = this.tensor
@@ -51,6 +57,7 @@ public open class DoubleTensorAlgebra :
         )
     }
 
+    @PerformancePitfall
     @Suppress("OVERRIDE_BY_INLINE")
     final override inline fun StructureND<Double>.mapIndexed(transform: DoubleField.(index: IntArray, Double) -> Double): DoubleTensor {
         val tensor = this.tensor
@@ -64,12 +71,13 @@ public open class DoubleTensorAlgebra :
         )
     }
 
+    @PerformancePitfall
     override fun zip(
         left: StructureND<Double>,
         right: StructureND<Double>,
-        transform: DoubleField.(Double, Double) -> Double
+        transform: DoubleField.(Double, Double) -> Double,
     ): DoubleTensor {
-        require(left.shape.contentEquals(right.shape)){
+        require(left.shape.contentEquals(right.shape)) {
             "The shapes in zip are not equal: left - ${left.shape}, right - ${right.shape}"
         }
         val leftTensor = left.tensor
@@ -115,7 +123,7 @@ public open class DoubleTensorAlgebra :
         TensorLinearStructure(shape).asSequence().map { DoubleField.initializer(it) }.toMutableList().toDoubleArray()
     )
 
-    override operator fun StructureND<Double>.get(i: Int): DoubleTensor {
+    override operator fun Tensor<Double>.get(i: Int): DoubleTensor {
         val lastShape = tensor.shape.drop(1).toIntArray()
         val newShape = if (lastShape.isNotEmpty()) lastShape else intArrayOf(1)
         val newStart = newShape.reduce(Int::times) * i + tensor.bufferStart
@@ -314,11 +322,11 @@ public open class DoubleTensorAlgebra :
         return DoubleTensor(shape, resBuffer)
     }
 
-    override fun StructureND<Double>.div(other: StructureND<Double>): DoubleTensor {
+    override fun StructureND<Double>.div(arg: StructureND<Double>): DoubleTensor {
-        checkShapesCompatible(tensor, other)
+        checkShapesCompatible(tensor, arg)
         val resBuffer = DoubleArray(tensor.numElements) { i ->
-            tensor.mutableBuffer.array()[other.tensor.bufferStart + i] /
+            tensor.mutableBuffer.array()[arg.tensor.bufferStart + i] /
-                    other.tensor.mutableBuffer.array()[other.tensor.bufferStart + i]
+                    arg.tensor.mutableBuffer.array()[arg.tensor.bufferStart + i]
         }
         return DoubleTensor(tensor.shape, resBuffer)
     }
@@ -329,11 +337,11 @@ public open class DoubleTensorAlgebra :
         }
     }
 
-    override fun Tensor<Double>.divAssign(other: StructureND<Double>) {
+    override fun Tensor<Double>.divAssign(arg: StructureND<Double>) {
-        checkShapesCompatible(tensor, other)
+        checkShapesCompatible(tensor, arg)
         for (i in 0 until tensor.numElements) {
             tensor.mutableBuffer.array()[tensor.bufferStart + i] /=
-                other.tensor.mutableBuffer.array()[tensor.bufferStart + i]
+                arg.tensor.mutableBuffer.array()[tensor.bufferStart + i]
         }
     }
 
@@ -344,7 +352,7 @@ public open class DoubleTensorAlgebra :
         return DoubleTensor(tensor.shape, resBuffer)
     }
 
-    override fun StructureND<Double>.transpose(i: Int, j: Int): DoubleTensor {
+    override fun Tensor<Double>.transpose(i: Int, j: Int): DoubleTensor {
         val ii = tensor.minusIndex(i)
         val jj = tensor.minusIndex(j)
         checkTranspose(tensor.dimension, ii, jj)
@@ -376,6 +384,7 @@ public open class DoubleTensorAlgebra :
     override fun Tensor<Double>.viewAs(other: StructureND<Double>): DoubleTensor =
         tensor.view(other.shape)
 
+    @PerformancePitfall
     override infix fun StructureND<Double>.dot(other: StructureND<Double>): DoubleTensor {
         if (tensor.shape.size == 1 && other.shape.size == 1) {
             return DoubleTensor(intArrayOf(1), doubleArrayOf(tensor.times(other).tensor.mutableBuffer.array().sum()))
@@ -413,14 +422,11 @@ public open class DoubleTensorAlgebra :
 
         for ((res, ab) in resTensor.matrixSequence().zip(newThis.matrixSequence().zip(newOther.matrixSequence()))) {
             val (a, b) = ab
-            dotHelper(a.as2D(), b.as2D(), res.as2D(), l, m1, n)
+            dotTo(a.as2D(), b.as2D(), res.as2D(), l, m1, n)
         }
 
         return if (penultimateDim) {
-            resTensor.view(
+            resTensor.view(resTensor.shape.dropLast(2).toIntArray() + intArrayOf(resTensor.shape.last()))
-                resTensor.shape.dropLast(2).toIntArray() +
-                        intArrayOf(resTensor.shape.last())
-            )
         } else if (lastDim) {
             resTensor.view(resTensor.shape.dropLast(1).toIntArray())
         } else {
@@ -432,7 +438,7 @@ public open class DoubleTensorAlgebra :
         diagonalEntries: Tensor<Double>,
         offset: Int,
         dim1: Int,
-        dim2: Int
+        dim2: Int,
     ): DoubleTensor {
         val n = diagonalEntries.shape.size
         val d1 = minusIndexFrom(n + 1, dim1)
@@ -568,14 +574,14 @@ public open class DoubleTensorAlgebra :
      */
     public fun Tensor<Double>.rowsByIndices(indices: IntArray): DoubleTensor = stack(indices.map { this[it] })
 
-    internal inline fun StructureND<Double>.fold(foldFunction: (DoubleArray) -> Double): Double =
+    private inline fun StructureND<Double>.fold(foldFunction: (DoubleArray) -> Double): Double =
         foldFunction(tensor.copyArray())
 
-    internal inline fun StructureND<Double>.foldDim(
+    private inline fun <reified R : Any> StructureND<Double>.foldDim(
-        foldFunction: (DoubleArray) -> Double,
         dim: Int,
         keepDim: Boolean,
-    ): DoubleTensor {
+        foldFunction: (DoubleArray) -> R,
+    ): BufferedTensor<R> {
         check(dim < dimension) { "Dimension $dim out of range $dimension" }
         val resShape = if (keepDim) {
             shape.take(dim).toIntArray() + intArrayOf(1) + shape.takeLast(dimension - dim - 1).toIntArray()
@@ -583,80 +589,75 @@ public open class DoubleTensorAlgebra :
             shape.take(dim).toIntArray() + shape.takeLast(dimension - dim - 1).toIntArray()
         }
         val resNumElements = resShape.reduce(Int::times)
-        val resTensor = DoubleTensor(resShape, DoubleArray(resNumElements) { 0.0 }, 0)
+        val init = foldFunction(DoubleArray(1) { 0.0 })
-        for (index in resTensor.indices.asSequence()) {
+        val resTensor = BufferedTensor(resShape,
+            MutableBuffer.auto(resNumElements) { init }, 0)
+        for (index in resTensor.indices) {
             val prefix = index.take(dim).toIntArray()
             val suffix = index.takeLast(dimension - dim - 1).toIntArray()
             resTensor[index] = foldFunction(DoubleArray(shape[dim]) { i ->
                 tensor[prefix + intArrayOf(i) + suffix]
             })
         }
-
         return resTensor
     }
 
     override fun StructureND<Double>.sum(): Double = tensor.fold { it.sum() }
 
     override fun StructureND<Double>.sum(dim: Int, keepDim: Boolean): DoubleTensor =
-        foldDim({ x -> x.sum() }, dim, keepDim)
+        foldDim(dim, keepDim) { x -> x.sum() }.toDoubleTensor()
 
     override fun StructureND<Double>.min(): Double = this.fold { it.minOrNull()!! }
 
     override fun StructureND<Double>.min(dim: Int, keepDim: Boolean): DoubleTensor =
-        foldDim({ x -> x.minOrNull()!! }, dim, keepDim)
+        foldDim(dim, keepDim) { x -> x.minOrNull()!! }.toDoubleTensor()
 
     override fun StructureND<Double>.max(): Double = this.fold { it.maxOrNull()!! }
 
     override fun StructureND<Double>.max(dim: Int, keepDim: Boolean): DoubleTensor =
-        foldDim({ x -> x.maxOrNull()!! }, dim, keepDim)
+        foldDim(dim, keepDim) { x -> x.maxOrNull()!! }.toDoubleTensor()
 
-    override fun StructureND<Double>.argMax(dim: Int, keepDim: Boolean): DoubleTensor =
+
-        foldDim({ x ->
+    override fun StructureND<Double>.argMax(dim: Int, keepDim: Boolean): IntTensor =
-            x.withIndex().maxByOrNull { it.value }?.index!!.toDouble()
+        foldDim(dim, keepDim) { x ->
-        }, dim, keepDim)
+            x.withIndex().maxByOrNull { it.value }?.index!!
+        }.toIntTensor()
 
 
     override fun StructureND<Double>.mean(): Double = this.fold { it.sum() / tensor.numElements }
 
-    override fun StructureND<Double>.mean(dim: Int, keepDim: Boolean): DoubleTensor =
+    override fun StructureND<Double>.mean(dim: Int, keepDim: Boolean): DoubleTensor = foldDim(dim, keepDim) { arr ->
-        foldDim(
+        check(dim < dimension) { "Dimension $dim out of range $dimension" }
-            { arr ->
+        arr.sum() / shape[dim]
-                check(dim < dimension) { "Dimension $dim out of range $dimension" }
+    }.toDoubleTensor()
-                arr.sum() / shape[dim]
-            },
-            dim,
-            keepDim
-        )
 
-    override fun StructureND<Double>.std(): Double = this.fold { arr ->
+    override fun StructureND<Double>.std(): Double = fold { arr ->
         val mean = arr.sum() / tensor.numElements
         sqrt(arr.sumOf { (it - mean) * (it - mean) } / (tensor.numElements - 1))
     }
 
     override fun StructureND<Double>.std(dim: Int, keepDim: Boolean): DoubleTensor = foldDim(
-        { arr ->
-            check(dim < dimension) { "Dimension $dim out of range $dimension" }
-            val mean = arr.sum() / shape[dim]
-            sqrt(arr.sumOf { (it - mean) * (it - mean) } / (shape[dim] - 1))
-        },
         dim,
         keepDim
-    )
+    ) { arr ->
+        check(dim < dimension) { "Dimension $dim out of range $dimension" }
+        val mean = arr.sum() / shape[dim]
+        sqrt(arr.sumOf { (it - mean) * (it - mean) } / (shape[dim] - 1))
+    }.toDoubleTensor()
 
-    override fun StructureND<Double>.variance(): Double = this.fold { arr ->
+    override fun StructureND<Double>.variance(): Double = fold { arr ->
         val mean = arr.sum() / tensor.numElements
         arr.sumOf { (it - mean) * (it - mean) } / (tensor.numElements - 1)
     }
 
     override fun StructureND<Double>.variance(dim: Int, keepDim: Boolean): DoubleTensor = foldDim(
-        { arr ->
-            check(dim < dimension) { "Dimension $dim out of range $dimension" }
-            val mean = arr.sum() / shape[dim]
-            arr.sumOf { (it - mean) * (it - mean) } / (shape[dim] - 1)
-        },
         dim,
         keepDim
-    )
+    ) { arr ->
+        check(dim < dimension) { "Dimension $dim out of range $dimension" }
+        val mean = arr.sum() / shape[dim]
+        arr.sumOf { (it - mean) * (it - mean) } / (shape[dim] - 1)
+    }.toDoubleTensor()
 
     private fun cov(x: DoubleTensor, y: DoubleTensor): Double {
         val n = x.shape[0]

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/IntTensor.kt

@@ -6,6 +6,7 @@
 package space.kscience.kmath.tensors.core
 
 import space.kscience.kmath.structures.IntBuffer
+import space.kscience.kmath.tensors.core.internal.array
 
 /**
  * Default [BufferedTensor] implementation for [Int] values
@@ -14,4 +15,7 @@ public class IntTensor internal constructor(
     shape: IntArray,
     buffer: IntArray,
     offset: Int = 0
-) : BufferedTensor<Int>(shape, IntBuffer(buffer), offset)
+) : BufferedTensor<Int>(shape, IntBuffer(buffer), offset){
+    public fun asDouble() : DoubleTensor =
+        DoubleTensor(shape, mutableBuffer.array().map{ it.toDouble()}.toDoubleArray(), bufferStart)
+}

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/TensorLinearStructure.kt

@@ -0,0 +1,74 @@
+/*
+ * 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.
+ */
+
+package space.kscience.kmath.tensors.core
+
+import space.kscience.kmath.nd.Strides
+import kotlin.math.max
+
+/**
+ * This [Strides] implementation follows the last dimension first convention
+ * For more information: https://numpy.org/doc/stable/reference/generated/numpy.ndarray.strides.html
+ *
+ * @param shape the shape of the tensor.
+ */
+public class TensorLinearStructure(override val shape: IntArray) : Strides() {
+    override val strides: IntArray
+        get() = stridesFromShape(shape)
+
+    override fun index(offset: Int): IntArray =
+        indexFromOffset(offset, strides, shape.size)
+
+    override val linearSize: Int
+        get() = shape.reduce(Int::times)
+
+    override fun equals(other: Any?): Boolean {
+        if (this === other) return true
+        if (other == null || this::class != other::class) return false
+
+        other as TensorLinearStructure
+
+        if (!shape.contentEquals(other.shape)) return false
+
+        return true
+    }
+
+    override fun hashCode(): Int {
+        return shape.contentHashCode()
+    }
+
+    public companion object {
+
+        public fun stridesFromShape(shape: IntArray): IntArray {
+            val nDim = shape.size
+            val res = IntArray(nDim)
+            if (nDim == 0)
+                return res
+
+            var current = nDim - 1
+            res[current] = 1
+
+            while (current > 0) {
+                res[current - 1] = max(1, shape[current]) * res[current]
+                current--
+            }
+            return res
+        }
+
+        public fun indexFromOffset(offset: Int, strides: IntArray, nDim: Int): IntArray {
+            val res = IntArray(nDim)
+            var current = offset
+            var strideIndex = 0
+
+            while (strideIndex < nDim) {
+                res[strideIndex] = (current / strides[strideIndex])
+                current %= strides[strideIndex]
+                strideIndex++
+            }
+            return res
+        }
+    }
+
+}

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/TensorLinearStructure.kt

@@ -1,71 +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 file.
- */
-
-package space.kscience.kmath.tensors.core.internal
-
-import space.kscience.kmath.nd.Strides
-import kotlin.math.max
-
-
-internal fun stridesFromShape(shape: IntArray): IntArray {
-    val nDim = shape.size
-    val res = IntArray(nDim)
-    if (nDim == 0)
-        return res
-
-    var current = nDim - 1
-    res[current] = 1
-
-    while (current > 0) {
-        res[current - 1] = max(1, shape[current]) * res[current]
-        current--
-    }
-    return res
-}
-
-internal fun indexFromOffset(offset: Int, strides: IntArray, nDim: Int): IntArray {
-    val res = IntArray(nDim)
-    var current = offset
-    var strideIndex = 0
-
-    while (strideIndex < nDim) {
-        res[strideIndex] = (current / strides[strideIndex])
-        current %= strides[strideIndex]
-        strideIndex++
-    }
-    return res
-}
-
-/**
- * This [Strides] implementation follows the last dimension first convention
- * For more information: https://numpy.org/doc/stable/reference/generated/numpy.ndarray.strides.html
- *
- * @param shape the shape of the tensor.
- */
-internal class TensorLinearStructure(override val shape: IntArray) : Strides() {
-    override val strides: IntArray
-        get() = stridesFromShape(shape)
-
-    override fun index(offset: Int): IntArray =
-        indexFromOffset(offset, strides, shape.size)
-
-    override val linearSize: Int
-        get() = shape.reduce(Int::times)
-
-    override fun equals(other: Any?): Boolean {
-        if (this === other) return true
-        if (other == null || this::class != other::class) return false
-
-        other as TensorLinearStructure
-
-        if (!shape.contentEquals(other.shape)) return false
-
-        return true
-    }
-
-    override fun hashCode(): Int {
-        return shape.contentHashCode()
-    }
-}

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt

@@ -53,7 +53,7 @@ internal val <T> BufferedTensor<T>.matrices: VirtualBuffer<BufferedTensor<T>>
 
 internal fun <T> BufferedTensor<T>.matrixSequence(): Sequence<BufferedTensor<T>> = matrices.asSequence()
 
-internal fun dotHelper(
+internal fun dotTo(
     a: MutableStructure2D<Double>,
     b: MutableStructure2D<Double>,
     res: MutableStructure2D<Double>,

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/tensorCastsUtils.kt

@@ -12,6 +12,7 @@ import space.kscience.kmath.tensors.api.Tensor
 import space.kscience.kmath.tensors.core.BufferedTensor
 import space.kscience.kmath.tensors.core.DoubleTensor
 import space.kscience.kmath.tensors.core.IntTensor
+import space.kscience.kmath.tensors.core.TensorLinearStructure
 
 internal fun BufferedTensor<Int>.asTensor(): IntTensor =
     IntTensor(this.shape, this.mutableBuffer.array(), this.bufferStart)

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/tensorAlgebraExtensions.kt

@@ -3,8 +3,11 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
  */
 
+@file:OptIn(PerformancePitfall::class)
+
 package space.kscience.kmath.tensors.core
 
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.nd.Shape
 import kotlin.jvm.JvmName
 

kmath-viktor/src/main/kotlin/space/kscience/kmath/viktor/ViktorFieldOpsND.kt

@@ -6,17 +6,20 @@
 package space.kscience.kmath.viktor
 
 import org.jetbrains.bio.viktor.F64Array
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.ExtendedFieldOps
 import space.kscience.kmath.operations.NumbersAddOps
+import space.kscience.kmath.operations.PowerOperations
 
 @OptIn(UnstableKMathAPI::class)
 @Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
 public open class ViktorFieldOpsND :
     FieldOpsND<Double, DoubleField>,
-    ExtendedFieldOps<StructureND<Double>> {
+    ExtendedFieldOps<StructureND<Double>>,
+    PowerOperations<StructureND<Double>> {
 
     public val StructureND<Double>.f64Buffer: F64Array
         get() = when (this) {
@@ -35,6 +38,7 @@ public open class ViktorFieldOpsND :
 
     override fun StructureND<Double>.unaryMinus(): StructureND<Double> = -1 * this
 
+    @PerformancePitfall
     override fun StructureND<Double>.map(transform: DoubleField.(Double) -> Double): ViktorStructureND =
         F64Array(*shape).apply {
             DefaultStrides(shape).asSequence().forEach { index ->
@@ -42,6 +46,7 @@ public open class ViktorFieldOpsND :
             }
         }.asStructure()
 
+    @PerformancePitfall
     override fun StructureND<Double>.mapIndexed(
         transform: DoubleField.(index: IntArray, Double) -> Double,
     ): ViktorStructureND = F64Array(*shape).apply {
@@ -50,6 +55,7 @@ public open class ViktorFieldOpsND :
         }
     }.asStructure()
 
+    @PerformancePitfall
     override fun zip(
         left: StructureND<Double>,
         right: StructureND<Double>,
@@ -110,7 +116,7 @@ public open class ViktorFieldOpsND :
 public val DoubleField.viktorAlgebra: ViktorFieldOpsND get() = ViktorFieldOpsND
 
 public open class ViktorFieldND(
-    override val shape: Shape
+    override val shape: Shape,
 ) : ViktorFieldOpsND(), FieldND<Double, DoubleField>, NumbersAddOps<StructureND<Double>> {
     override val zero: ViktorStructureND by lazy { F64Array.full(init = 0.0, shape = shape).asStructure() }
     override val one: ViktorStructureND by lazy { F64Array.full(init = 1.0, shape = shape).asStructure() }