diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt
index 85150ad32..d1101d503 100644
--- a/benchmarks/src/main/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt
+++ b/benchmarks/src/main/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt
@@ -25,7 +25,7 @@ fun main(args: Array<String>) {
         bufferedField.run {
             var res: NDBuffer<Double> = one
             repeat(n) {
-                res += 1.0
+                res += one
             }
         }
     }
@@ -34,7 +34,7 @@ fun main(args: Array<String>) {
 
 
     val elementTime = measureTimeMillis {
-        var res = bufferedField.produce { one }
+        var res = bufferedField.run{one.toElement()}
         repeat(n) {
             res += 1.0
         }
@@ -77,7 +77,7 @@ fun main(args: Array<String>) {
         genericField.run {
             var res: NDBuffer<Double> = one
             repeat(n) {
-                res += 1.0
+                res += one
             }
         }
     }
diff --git a/benchmarks/src/main/kotlin/scientifik/kmath/structures/StructureWriteBenchmark.kt b/benchmarks/src/main/kotlin/scientifik/kmath/structures/StructureWriteBenchmark.kt
index 7bb4c8267..a2be2aa9c 100644
--- a/benchmarks/src/main/kotlin/scientifik/kmath/structures/StructureWriteBenchmark.kt
+++ b/benchmarks/src/main/kotlin/scientifik/kmath/structures/StructureWriteBenchmark.kt
@@ -1,5 +1,6 @@
 package scientifik.kmath.structures
 
+import scientifik.kmath.structures.Buffer.Companion.DoubleBufferFactory
 import kotlin.system.measureTimeMillis
 
 
diff --git a/doc/algebra.md b/doc/algebra.md
new file mode 100644
index 000000000..2888f9484
--- /dev/null
+++ b/doc/algebra.md
@@ -0,0 +1,70 @@
+# Algebra and algebra elements
+
+The mathematical operations in `kmath` are generally separated from mathematical objects.
+This means that in order to perform an operation, say `+`, one needs two objects of a type `T` and
+and algebra context which defines appropriate operation, say `Space<T>`. Next one needs to run actual operation
+in the context:
+
+```kotlin
+val a: T
+val b: T
+val space: Space<T>
+
+val c = space.run{a + b}
+```
+
+From the first glance, this distinction seems to be a needless complication, but in fact one needs
+to remember that in mathematics, one could define different operations on the same objects. For example,
+one could use different types of geometry for vectors.
+
+## Algebra hierarchy
+
+Mathematical contexts have the following hierarchy:
+
+**Space** <- **Ring** <- **Field**
+
+All classes follow abstract mathematical constructs.
+[Space](http://mathworld.wolfram.com/Space.html) defines `zero` element, addition operation and multiplication by constant,
+[Ring](http://mathworld.wolfram.com/Ring.html) adds multiplication and unit `one` element,
+[Field](http://mathworld.wolfram.com/Field.html) adds division operation.
+
+Typical case of `Field` is the `RealField` which works on doubles. And typical case of `Space` is a `VectorSpace`.
+
+In some cases algebra context could hold additional operation like `exp` or `sin`, in this case it inherits appropriate
+interface. Also a context could have an operation which produces an element outside of its context. For example
+`Matrix` `dot` operation produces a matrix with new dimensions which could not be compatible with initial matrix in
+terms of linear operations.
+
+## Algebra element
+
+In order to achieve more familiar behavior (where you apply operations directly to mathematica objects), without involving contexts
+`kmath` introduces special type objects called `MathElement`. A `MathElement` is basically some object coupled to
+a mathematical context. For example `Complex` is the pair of real numbers representing real and imaginary parts,
+but it also holds reference to the `ComplexField` singleton which allows to perform direct operations on `Complex`
+numbers without explicit involving the context like:
+
+```kotlin
+    val c1 = Complex(1.0, 1.0)
+    val c2 = Complex(1.0, -1.0)
+    val c3 = c1 + c2 + 3.0.toComplex()
+    //or with field notation:
+    val c4 = ComplexField.run{c1 + i - 2.0}
+```
+
+Both notations have their pros and cons.
+
+The hierarchy for algebra elements follows the hierarchy for the corresponding algebra.
+
+**MathElement** <- **SpaceElement** <- **RingElement** <- **FieldElement**
+
+**MathElement** is the generic common ancestor of the class with context.
+
+One important distinction between algebra elements and algebra contexts is that algebra element has three type parameters:
+
+1. The type of elements, field operates on.
+2. The self-type of the element returned from operation (must be algebra element).
+3. The type of the algebra over first type-parameter.
+
+The middle type is needed in case algebra members do not store context. For example, it is not possible to add
+a context to regular `Double`. The element performs automatic conversions from context types and back.
+One should used context operations in all important places. The performance of element operations is not guaranteed.
diff --git a/doc/nd-performance.md b/doc/nd-performance.md
new file mode 100644
index 000000000..8ddfc355b
--- /dev/null
+++ b/doc/nd-performance.md
@@ -0,0 +1,13 @@
+# Performance for n-dimensional structures operations
+
+One of the most sought after features of mathematical libraries is the high-performance operations on n-dimensional
+structures. In `kmath` performance depends on which particular context was used for operation.
+
+Let us consider following contexts:
+```kotlin
+    // automatically build context
+    val bufferedField = NDField.auto(intArrayOf(dim, dim), RealField)
+    val specializedField = NDField.real(intArrayOf(dim, dim))
+    val genericField = NDField.buffered(intArrayOf(dim, dim), RealField)
+    val lazyNDField = NDField.lazy(intArrayOf(dim, dim), RealField)
+```
\ No newline at end of file
diff --git a/doc/operations.md b/doc/operations.md
index e00b38ecc..ff4a407ee 100644
--- a/doc/operations.md
+++ b/doc/operations.md
@@ -20,7 +20,7 @@ val c3 = ComplexField.run{ c1 - i*2.0}
 ```
 
 **Note**: In theory it is possible to add behaviors directly to the context, but currently kotlin syntax does not support
-that. Watch [KT-10468](https://youtrack.jetbrains.com/issue/KT-10468) for updates.
+that. Watch [KT-10468](https://youtrack.jetbrains.com/issue/KT-10468) and [KEEP-176](https://github.com/Kotlin/KEEP/pull/176) for updates.
 
 ## Nested fields
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
index 5618b3100..402947288 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
@@ -2,7 +2,11 @@ package scientifik.kmath.linear
 
 import scientifik.kmath.operations.Field
 import scientifik.kmath.operations.RealField
-import scientifik.kmath.structures.*
+import scientifik.kmath.structures.MutableBuffer.Companion.boxing
+import scientifik.kmath.structures.MutableNDStructure
+import scientifik.kmath.structures.NDStructure
+import scientifik.kmath.structures.get
+import scientifik.kmath.structures.mutableNdStructure
 import kotlin.math.absoluteValue
 
 /**
@@ -106,7 +110,7 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
         //TODO fix performance
         val lu: MutableNDStructure<T> = mutableNdStructure(
             intArrayOf(matrix.numRows, matrix.numCols),
-            ::boxingMutableBuffer
+            ::boxing
         ) { index: IntArray -> matrix[index[0], index[1]] }
 
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
index 11f1f7dd8..72a859111 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
@@ -4,8 +4,9 @@ import scientifik.kmath.operations.Field
 import scientifik.kmath.operations.Norm
 import scientifik.kmath.operations.RealField
 import scientifik.kmath.operations.Ring
+import scientifik.kmath.structures.Buffer.Companion.boxing
 import scientifik.kmath.structures.asSequence
-import scientifik.kmath.structures.boxingBuffer
+
 
 
 /**
@@ -52,7 +53,7 @@ fun <T : Any, R : Ring<T>> Vector<T, R>.toMatrix(): Matrix<T, R> {
 //        matrix(size, 1, context.field) { i, j -> get(i) }
 //    }
     //return Matrix.of(size, 1, context.space) { i, _ -> get(i) }
-    return StructureMatrixSpace(size, 1, context.space, ::boxingBuffer).produce { i, _ -> get(i) }
+    return StructureMatrixSpace(size, 1, context.space, ::boxing).produce { i, _ -> get(i) }
 }
 
 object VectorL2Norm : Norm<Vector<out Number, *>, Double> {
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
index 1fc48264a..2a5765e30 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
@@ -5,6 +5,9 @@ import scientifik.kmath.operations.Ring
 import scientifik.kmath.operations.Space
 import scientifik.kmath.operations.SpaceElement
 import scientifik.kmath.structures.*
+import scientifik.kmath.structures.Buffer.Companion.DoubleBufferFactory
+import scientifik.kmath.structures.Buffer.Companion.auto
+import scientifik.kmath.structures.Buffer.Companion.boxing
 
 
 interface MatrixSpace<T : Any, R : Ring<T>> : Space<Matrix<T, R>> {
@@ -52,14 +55,14 @@ interface MatrixSpace<T : Any, R : Ring<T>> : Space<Matrix<T, R>> {
             rows: Int,
             columns: Int,
             ring: R,
-            bufferFactory: BufferFactory<T> = ::boxingBuffer
+            bufferFactory: BufferFactory<T> = ::boxing
         ): MatrixSpace<T, R> = StructureMatrixSpace(rows, columns, ring, bufferFactory)
 
         /**
          * Automatic buffered matrix, unboxed if it is possible
          */
         inline fun <reified T : Any, R : Ring<T>> smart(rows: Int, columns: Int, ring: R): MatrixSpace<T, R> =
-            buffered(rows, columns, ring, ::autoBuffer)
+            buffered(rows, columns, ring, ::auto)
     }
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
index 4cbed0b42..a673c0197 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
@@ -3,7 +3,9 @@ package scientifik.kmath.linear
 import scientifik.kmath.operations.RealField
 import scientifik.kmath.operations.Space
 import scientifik.kmath.operations.SpaceElement
-import scientifik.kmath.structures.*
+import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.BufferFactory
+import scientifik.kmath.structures.asSequence
 
 typealias Point<T> = Buffer<T>
 
@@ -40,7 +42,7 @@ interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
          * Non-boxing double vector space
          */
         fun real(size: Int): BufferVectorSpace<Double, RealField> {
-            return realSpaceCache.getOrPut(size) { BufferVectorSpace(size, RealField, DoubleBufferFactory) }
+            return realSpaceCache.getOrPut(size) { BufferVectorSpace(size, RealField, Buffer.DoubleBufferFactory) }
         }
 
         /**
@@ -49,14 +51,14 @@ interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
         fun <T : Any, S : Space<T>> buffered(
             size: Int,
             space: S,
-            bufferFactory: BufferFactory<T> = ::boxingBuffer
+            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
         ): VectorSpace<T, S> = BufferVectorSpace(size, space, bufferFactory)
 
         /**
          * Automatic buffered vector, unboxed if it is possible
          */
         inline fun <reified T : Any, S : Space<T>> smart(size: Int, space: S): VectorSpace<T, S> =
-            buffered(size, space, ::autoBuffer)
+            buffered(size, space, Buffer.Companion::auto)
     }
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt
index 51d0815ec..219af01fd 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt
@@ -68,14 +68,14 @@ interface Ring<T> : Space<T> {
 
     operator fun T.times(b: T): T = multiply(this, b)
 
-//    operator fun T.plus(b: Number) = this.plus(b * one)
-//    operator fun Number.plus(b: T) = b + this
-//
-//    operator fun T.minus(b: Number) = this.minus(b * one)
-//    operator fun Number.minus(b: T) = -b + this
+    operator fun T.plus(b: Number) = this.plus(b * one)
+    operator fun Number.plus(b: T) = b + this
+
+    operator fun T.minus(b: Number) = this.minus(b * one)
+    operator fun Number.minus(b: T) = -b + this
 }
 
-abstract class AbstractRing<T: Any> : AbstractSpace<T>(), Ring<T> {
+abstract class AbstractRing<T : Any> : AbstractSpace<T>(), Ring<T> {
     final override operator fun T.times(b: T): T = multiply(this, b)
 }
 
@@ -89,7 +89,7 @@ interface Field<T> : Ring<T> {
     operator fun Number.div(b: T) = this * divide(one, b)
 }
 
-abstract class AbstractField<T: Any> : AbstractRing<T>(), Field<T> {
+abstract class AbstractField<T : Any> : AbstractRing<T>(), Field<T> {
     final override operator fun T.div(b: T): T = divide(this, b)
     final override operator fun Number.div(b: T) = this * divide(one, b)
 }
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt
index 60b0e9869..018675be3 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt
@@ -22,6 +22,9 @@ abstract class StridedNDField<T, F : Field<T>>(shape: IntArray, elementField: F)
             produce { index -> get(index) }
         }
     }
+
+    fun NDBuffer<T>.toElement(): StridedNDElement<T, F> =
+        StridedNDElement(this@StridedNDField, buffer)
 }
 
 
@@ -40,14 +43,14 @@ class BufferNDField<T, F : Field<T>>(
     override val zero by lazy { produce { zero } }
     override val one by lazy { produce { one } }
 
-    override fun produce(initializer: F.(IntArray) -> T): BufferNDElement<T, F> =
-        BufferNDElement(
+    override fun produce(initializer: F.(IntArray) -> T): StridedNDElement<T, F> =
+        StridedNDElement(
             this,
             buildBuffer(strides.linearSize) { offset -> elementField.initializer(strides.index(offset)) })
 
-    override fun map(arg: NDBuffer<T>, transform: F.(T) -> T): BufferNDElement<T, F> {
+    override fun map(arg: NDBuffer<T>, transform: F.(T) -> T): StridedNDElement<T, F> {
         check(arg)
-        return BufferNDElement(
+        return StridedNDElement(
             this,
             buildBuffer(arg.strides.linearSize) { offset -> elementField.transform(arg.buffer[offset]) })
     }
@@ -55,9 +58,9 @@ class BufferNDField<T, F : Field<T>>(
     override fun mapIndexed(
         arg: NDBuffer<T>,
         transform: F.(index: IntArray, T) -> T
-    ): BufferNDElement<T, F> {
+    ): StridedNDElement<T, F> {
         check(arg)
-        return BufferNDElement(
+        return StridedNDElement(
             this,
             buildBuffer(arg.strides.linearSize) { offset ->
                 elementField.transform(
@@ -71,17 +74,17 @@ class BufferNDField<T, F : Field<T>>(
         a: NDBuffer<T>,
         b: NDBuffer<T>,
         transform: F.(T, T) -> T
-    ): BufferNDElement<T, F> {
+    ): StridedNDElement<T, F> {
         check(a, b)
-        return BufferNDElement(
+        return StridedNDElement(
             this,
             buildBuffer(strides.linearSize) { offset -> elementField.transform(a.buffer[offset], b.buffer[offset]) })
     }
 }
 
-class BufferNDElement<T, F : Field<T>>(override val context: StridedNDField<T, F>, override val buffer: Buffer<T>) :
+class StridedNDElement<T, F : Field<T>>(override val context: StridedNDField<T, F>, override val buffer: Buffer<T>) :
     NDBuffer<T>,
-    FieldElement<NDBuffer<T>, BufferNDElement<T, F>, StridedNDField<T, F>>,
+    FieldElement<NDBuffer<T>, StridedNDElement<T, F>, StridedNDField<T, F>>,
     NDElement<T, F> {
 
     override val elementField: F
@@ -90,8 +93,8 @@ class BufferNDElement<T, F : Field<T>>(override val context: StridedNDField<T, F
     override fun unwrap(): NDBuffer<T> =
         this
 
-    override fun NDBuffer<T>.wrap(): BufferNDElement<T, F> =
-        BufferNDElement(context, this.buffer)
+    override fun NDBuffer<T>.wrap(): StridedNDElement<T, F> =
+        StridedNDElement(context, this.buffer)
 
     override val strides
         get() = context.strides
@@ -106,42 +109,42 @@ class BufferNDElement<T, F : Field<T>>(override val context: StridedNDField<T, F
         strides.indices().map { it to get(it) }
 
     override fun map(action: F.(T) -> T) =
-        context.run { map(this@BufferNDElement, action) }.wrap()
+        context.run { map(this@StridedNDElement, action) }.wrap()
 
     override fun mapIndexed(transform: F.(index: IntArray, T) -> T) =
-        context.run { mapIndexed(this@BufferNDElement, transform) }.wrap()
+        context.run { mapIndexed(this@StridedNDElement, transform) }.wrap()
 }
 
 /**
  * Element by element application of any operation on elements to the whole array. Just like in numpy
  */
-operator fun <T : Any, F : Field<T>> Function1<T, T>.invoke(ndElement: BufferNDElement<T, F>) =
+operator fun <T : Any, F : Field<T>> Function1<T, T>.invoke(ndElement: StridedNDElement<T, F>) =
     ndElement.context.run { ndElement.map { invoke(it) } }
 
 /* plus and minus */
 
 /**
- * Summation operation for [BufferNDElement] and single element
+ * Summation operation for [StridedNDElement] and single element
  */
-operator fun <T : Any, F : Field<T>> BufferNDElement<T, F>.plus(arg: T) =
+operator fun <T : Any, F : Field<T>> StridedNDElement<T, F>.plus(arg: T) =
     context.run { map { it + arg } }
 
 /**
- * Subtraction operation between [BufferNDElement] and single element
+ * Subtraction operation between [StridedNDElement] and single element
  */
-operator fun <T : Any, F : Field<T>> BufferNDElement<T, F>.minus(arg: T) =
+operator fun <T : Any, F : Field<T>> StridedNDElement<T, F>.minus(arg: T) =
     context.run { map { it - arg } }
 
 /* prod and div */
 
 /**
- * Product operation for [BufferNDElement] and single element
+ * Product operation for [StridedNDElement] and single element
  */
-operator fun <T : Any, F : Field<T>> BufferNDElement<T, F>.times(arg: T) =
+operator fun <T : Any, F : Field<T>> StridedNDElement<T, F>.times(arg: T) =
     context.run { map { it * arg } }
 
 /**
- * Division operation between [BufferNDElement] and single element
+ * Division operation between [StridedNDElement] and single element
  */
-operator fun <T : Any, F : Field<T>> BufferNDElement<T, F>.div(arg: T) =
+operator fun <T : Any, F : Field<T>> StridedNDElement<T, F>.div(arg: T) =
     context.run { map { it / arg } }
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
index 2c08427b3..d93a0d46e 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
@@ -1,6 +1,10 @@
 package scientifik.kmath.structures
 
 
+typealias BufferFactory<T> = (Int, (Int) -> T) -> Buffer<T>
+typealias MutableBufferFactory<T> = (Int, (Int) -> T) -> MutableBuffer<T>
+
+
 /**
  * A generic random access structure for both primitives and objects
  */
@@ -14,6 +18,31 @@ interface Buffer<T> {
 
     fun contentEquals(other: Buffer<*>): Boolean =
         asSequence().mapIndexed { index, value -> value == other[index] }.all { it }
+
+    companion object {
+
+        /**
+         * Create a boxing buffer of given type
+         */
+        inline fun <T> boxing(size: Int, initializer: (Int) -> T): Buffer<T> = ListBuffer(List(size, initializer))
+
+        /**
+         * Create most appropriate immutable buffer for given type avoiding boxing wherever possible
+         */
+        @Suppress("UNCHECKED_CAST")
+        inline fun <reified T : Any> auto(size: Int, initializer: (Int) -> T): Buffer<T> {
+            return when (T::class) {
+                Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as Buffer<T>
+                Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as Buffer<T>
+                Long::class -> LongBuffer(LongArray(size) { initializer(it) as Long }) as Buffer<T>
+                else -> boxing(size, initializer)
+            }
+        }
+
+        val DoubleBufferFactory: BufferFactory<Double> = { size, initializer -> DoubleBuffer(DoubleArray(size, initializer)) }
+        val IntBufferFactory: BufferFactory<Int> = { size, initializer -> IntBuffer(IntArray(size, initializer)) }
+        val LongBufferFactory: BufferFactory<Long> = { size, initializer -> LongBuffer(LongArray(size, initializer)) }
+    }
 }
 
 fun <T> Buffer<T>.asSequence(): Sequence<T> = iterator().asSequence()
@@ -27,6 +56,27 @@ interface MutableBuffer<T> : Buffer<T> {
      * A shallow copy of the buffer
      */
     fun copy(): MutableBuffer<T>
+
+    companion object {
+        /**
+         * Create a boxing mutable buffer of given type
+         */
+        inline fun <T : Any> boxing(size: Int, initializer: (Int) -> T): MutableBuffer<T> =
+            MutableListBuffer(MutableList(size, initializer))
+
+        /**
+         * Create most appropriate mutable buffer for given type avoiding boxing wherever possible
+         */
+        @Suppress("UNCHECKED_CAST")
+        inline fun <reified T : Any> auto(size: Int, initializer: (Int) -> T): MutableBuffer<T> {
+            return when (T::class) {
+                Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as MutableBuffer<T>
+                Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as MutableBuffer<T>
+                Long::class -> LongBuffer(LongArray(size) { initializer(it) as Long }) as MutableBuffer<T>
+                else -> boxing(size, initializer)
+            }
+        }
+    }
 }
 
 
@@ -130,48 +180,3 @@ fun <T> Buffer<T>.asReadOnly(): Buffer<T> = if (this is MutableBuffer) {
 } else {
     this
 }
-
-/**
- * Create a boxing buffer of given type
- */
-inline fun <T> boxingBuffer(size: Int, initializer: (Int) -> T): Buffer<T> = ListBuffer(List(size, initializer))
-
-/**
- * Create most appropriate immutable buffer for given type avoiding boxing wherever possible
- */
-@Suppress("UNCHECKED_CAST")
-inline fun <reified T : Any> autoBuffer(size: Int, initializer: (Int) -> T): Buffer<T> {
-    return when (T::class) {
-        Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as Buffer<T>
-        Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as Buffer<T>
-        Long::class -> LongBuffer(LongArray(size) { initializer(it) as Long }) as Buffer<T>
-        else -> boxingBuffer(size, initializer)
-    }
-}
-
-/**
- * Create a boxing mutable buffer of given type
- */
-inline fun <T : Any> boxingMutableBuffer(size: Int, initializer: (Int) -> T): MutableBuffer<T> =
-    MutableListBuffer(MutableList(size, initializer))
-
-/**
- * Create most appropriate mutable buffer for given type avoiding boxing wherever possible
- */
-@Suppress("UNCHECKED_CAST")
-inline fun <reified T : Any> autoMutableBuffer(size: Int, initializer: (Int) -> T): MutableBuffer<T> {
-    return when (T::class) {
-        Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as MutableBuffer<T>
-        Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as MutableBuffer<T>
-        Long::class -> LongBuffer(LongArray(size) { initializer(it) as Long }) as MutableBuffer<T>
-        else -> boxingMutableBuffer(size, initializer)
-    }
-}
-
-typealias BufferFactory<T> = (Int, (Int) -> T) -> Buffer<T>
-typealias MutableBufferFactory<T> = (Int, (Int) -> T) -> MutableBuffer<T>
-
-val DoubleBufferFactory: BufferFactory<Double> = { size, initializer -> DoubleBuffer(DoubleArray(size, initializer)) }
-val IntBufferFactory: BufferFactory<Int> = { size, initializer -> IntBuffer(IntArray(size, initializer)) }
-val LongBufferFactory: BufferFactory<Long> = { size, initializer -> LongBuffer(LongArray(size, initializer)) }
-
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt
index 022a4e61f..18ada3449 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDElement.kt
@@ -10,47 +10,46 @@ interface NDElement<T, F : Field<T>> : NDStructure<T> {
 
     fun mapIndexed(transform: F.(index: IntArray, T) -> T): NDElement<T, F>
     fun map(action: F.(T) -> T) = mapIndexed { _, value -> action(value) }
-}
+
+    companion object {
+        /**
+         * Create a optimized NDArray of doubles
+         */
+        fun real(shape: IntArray, initializer: RealField.(IntArray) -> Double = { 0.0 }) =
+            NDField.real(shape).produce(initializer)
 
 
-object NDElements {
-    /**
-     * Create a optimized NDArray of doubles
-     */
-    fun real(shape: IntArray, initializer: RealField.(IntArray) -> Double = { 0.0 }) =
-        NDField.real(shape).produce(initializer)
+        fun real1D(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }) =
+            real(intArrayOf(dim)) { initializer(it[0]) }
 
 
-    fun real1D(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }) =
-        real(intArrayOf(dim)) { initializer(it[0]) }
+        fun real2D(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }) =
+            real(intArrayOf(dim1, dim2)) { initializer(it[0], it[1]) }
+
+        fun real3D(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }) =
+            real(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
 
 
-    fun real2D(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }) =
-        real(intArrayOf(dim1, dim2)) { initializer(it[0], it[1]) }
+        /**
+         * Simple boxing NDArray
+         */
+        fun <T : Any, F : Field<T>> buffered(
+            shape: IntArray,
+            field: F,
+            initializer: F.(IntArray) -> T
+        ): StridedNDElement<T, F> {
+            val ndField = BufferNDField(shape, field, Buffer.Companion::boxing)
+            return ndField.produce(initializer)
+        }
 
-    fun real3D(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }) =
-        real(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
-
-
-    /**
-     * Simple boxing NDArray
-     */
-    fun <T : Any, F : Field<T>> buffered(
-        shape: IntArray,
-        field: F,
-        initializer: F.(IntArray) -> T
-    ): BufferNDElement<T, F> {
-        val ndField = BufferNDField(shape, field, ::boxingBuffer)
-        return ndField.produce(initializer)
-    }
-
-    inline fun <reified T : Any, F : Field<T>> auto(
-        shape: IntArray,
-        field: F,
-        noinline initializer: F.(IntArray) -> T
-    ): BufferNDElement<T, F> {
-        val ndField = NDField.auto(shape, field)
-        return ndField.produce(initializer)
+        inline fun <reified T : Any, F : Field<T>> auto(
+            shape: IntArray,
+            field: F,
+            noinline initializer: F.(IntArray) -> T
+        ): StridedNDElement<T, F> {
+            val ndField = NDField.auto(shape, field)
+            return StridedNDElement(ndField, ndField.produce(initializer).buffer)
+        }
     }
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
index c1abf8acd..b859b6dd3 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
@@ -2,6 +2,7 @@ package scientifik.kmath.structures
 
 import scientifik.kmath.operations.AbstractField
 import scientifik.kmath.operations.Field
+import scientifik.kmath.structures.Buffer.Companion.boxing
 
 /**
  * An exception is thrown when the expected ans actual shape of NDArray differs
@@ -36,13 +37,18 @@ interface NDField<T, F : Field<T>, N : NDStructure<T>> : Field<N> {
          * Create a nd-field with boxing generic buffer
          */
         fun <T : Any, F : Field<T>> buffered(shape: IntArray, field: F) =
-            BufferNDField(shape, field, ::boxingBuffer)
+            BufferNDField(shape, field, Buffer.Companion::boxing)
 
         /**
          * Create a most suitable implementation for nd-field using reified class.
          */
-        inline fun <reified T : Any, F : Field<T>> auto(shape: IntArray, field: F) =
-            BufferNDField(shape, field, ::autoBuffer)
+        inline fun <reified T : Any, F : Field<T>> auto(shape: IntArray, field: F): StridedNDField<T, F> =
+            if (T::class == Double::class) {
+                @Suppress("UNCHECKED_CAST")
+                real(shape) as StridedNDField<T, F>
+            } else {
+                BufferNDField(shape, field, Buffer.Companion::auto)
+            }
     }
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
index 0852271e4..74e3a59a2 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
@@ -163,7 +163,7 @@ data class BufferNDStructure<T>(
  * Transform structure to a new structure using provided [BufferFactory] and optimizing if argument is [BufferNDStructure]
  */
 inline fun <T, reified R : Any> NDStructure<T>.mapToBuffer(
-    factory: BufferFactory<R> = ::autoBuffer,
+    factory: BufferFactory<R> = Buffer.Companion::auto,
     crossinline transform: (T) -> R
 ): BufferNDStructure<R> {
     return if (this is BufferNDStructure<T>) {
@@ -179,16 +179,16 @@ inline fun <T, reified R : Any> NDStructure<T>.mapToBuffer(
  *
  * Strides should be reused if possible
  */
-fun <T> ndStructure(strides: Strides, bufferFactory: BufferFactory<T> = ::boxingBuffer, initializer: (IntArray) -> T) =
+fun <T> ndStructure(strides: Strides, bufferFactory: BufferFactory<T> = Buffer.Companion::boxing, initializer: (IntArray) -> T) =
     BufferNDStructure(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
 
 /**
  * Inline create NDStructure with non-boxing buffer implementation if it is possible
  */
 inline fun <reified T : Any> inlineNDStructure(strides: Strides, crossinline initializer: (IntArray) -> T) =
-    BufferNDStructure(strides, autoBuffer(strides.linearSize) { i -> initializer(strides.index(i)) })
+    BufferNDStructure(strides, Buffer.Companion.auto(strides.linearSize) { i -> initializer(strides.index(i)) })
 
-fun <T> ndStructure(shape: IntArray, bufferFactory: BufferFactory<T> = ::boxingBuffer, initializer: (IntArray) -> T) =
+fun <T> ndStructure(shape: IntArray, bufferFactory: BufferFactory<T> = Buffer.Companion::boxing, initializer: (IntArray) -> T) =
     ndStructure(DefaultStrides(shape), bufferFactory, initializer)
 
 inline fun <reified T : Any> inlineNdStructure(shape: IntArray, crossinline initializer: (IntArray) -> T) =
@@ -216,17 +216,17 @@ class MutableBufferNDStructure<T>(
  */
 fun <T : Any> mutableNdStructure(
     strides: Strides,
-    bufferFactory: MutableBufferFactory<T> = ::boxingMutableBuffer,
+    bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::boxing,
     initializer: (IntArray) -> T
 ) =
     MutableBufferNDStructure(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
 
 inline fun <reified T : Any> inlineMutableNdStructure(strides: Strides, crossinline initializer: (IntArray) -> T) =
-    MutableBufferNDStructure(strides, autoMutableBuffer(strides.linearSize) { i -> initializer(strides.index(i)) })
+    MutableBufferNDStructure(strides, MutableBuffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) })
 
 fun <T : Any> mutableNdStructure(
     shape: IntArray,
-    bufferFactory: MutableBufferFactory<T> = ::boxingMutableBuffer,
+    bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::boxing,
     initializer: (IntArray) -> T
 ) =
     mutableNdStructure(DefaultStrides(shape), bufferFactory, initializer)
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
index 82d399436..b8bab7c07 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
@@ -2,7 +2,7 @@ package scientifik.kmath.structures
 
 import scientifik.kmath.operations.RealField
 
-typealias RealNDElement = BufferNDElement<Double, RealField>
+typealias RealNDElement = StridedNDElement<Double, RealField>
 
 class RealNDField(shape: IntArray) :
     StridedNDField<Double, RealField>(shape, RealField),
@@ -21,20 +21,20 @@ class RealNDField(shape: IntArray) :
     ): RealNDElement {
         check(arg)
         val array = buildBuffer(arg.strides.linearSize) { offset -> RealField.transform(arg.buffer[offset]) }
-        return BufferNDElement(this, array)
+        return StridedNDElement(this, array)
     }
 
     override fun produce(initializer: RealField.(IntArray) -> Double): RealNDElement {
         val array = buildBuffer(strides.linearSize) { offset -> elementField.initializer(strides.index(offset)) }
-        return BufferNDElement(this, array)
+        return StridedNDElement(this, array)
     }
 
     override fun mapIndexed(
         arg: NDBuffer<Double>,
         transform: RealField.(index: IntArray, Double) -> Double
-    ): BufferNDElement<Double, RealField> {
+    ): StridedNDElement<Double, RealField> {
         check(arg)
-        return BufferNDElement(
+        return StridedNDElement(
             this,
             buildBuffer(arg.strides.linearSize) { offset ->
                 elementField.transform(
@@ -48,9 +48,9 @@ class RealNDField(shape: IntArray) :
         a: NDBuffer<Double>,
         b: NDBuffer<Double>,
         transform: RealField.(Double, Double) -> Double
-    ): BufferNDElement<Double, RealField> {
+    ): StridedNDElement<Double, RealField> {
         check(a, b)
-        return BufferNDElement(
+        return StridedNDElement(
             this,
             buildBuffer(strides.linearSize) { offset -> elementField.transform(a.buffer[offset], b.buffer[offset]) })
     }
@@ -80,7 +80,7 @@ class RealNDField(shape: IntArray) :
  */
 inline fun StridedNDField<Double, RealField>.produceInline(crossinline initializer: RealField.(Int) -> Double): RealNDElement {
     val array = DoubleArray(strides.linearSize) { offset -> elementField.initializer(offset) }
-    return BufferNDElement(this, DoubleBuffer(array))
+    return StridedNDElement(this, DoubleBuffer(array))
 }
 
 /**
@@ -93,13 +93,13 @@ operator fun Function1<Double, Double>.invoke(ndElement: RealNDElement) =
 /* plus and minus */
 
 /**
- * Summation operation for [BufferNDElement] and single element
+ * Summation operation for [StridedNDElement] and single element
  */
 operator fun RealNDElement.plus(arg: Double) =
     context.produceInline { i -> buffer[i] + arg }
 
 /**
- * Subtraction operation between [BufferNDElement] and single element
+ * Subtraction operation between [StridedNDElement] and single element
  */
 operator fun RealNDElement.minus(arg: Double) =
     context.produceInline { i -> buffer[i] - arg }
\ No newline at end of file
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/ExpressionFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/ExpressionFieldTest.kt
index 5c5002e06..033b2792f 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/ExpressionFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/ExpressionFieldTest.kt
@@ -12,7 +12,7 @@ class ExpressionFieldTest {
         val context = ExpressionField(RealField)
         val expression = with(context) {
             val x = variable("x", 2.0)
-            x * x + 2 * x + 1.0
+            x * x + 2 * x + one
         }
         assertEquals(expression("x" to 1.0), 4.0)
         assertEquals(expression(), 9.0)
@@ -44,7 +44,7 @@ class ExpressionFieldTest {
     fun valueExpression() {
         val expressionBuilder: ExpressionField<Double>.() -> Expression<Double> = {
             val x = variable("x")
-            x * x + 2 * x + 1.0
+            x * x + 2 * x + one
         }
 
         val expression = ExpressionField(RealField).expressionBuilder()
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NDFieldTest.kt
index 637147629..39cce5c67 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NDFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NDFieldTest.kt
@@ -7,7 +7,7 @@ import kotlin.test.assertEquals
 class NDFieldTest {
     @Test
     fun testStrides() {
-        val ndArray = NDElements.real(intArrayOf(10, 10)) { (it[0] + it[1]).toDouble() }
+        val ndArray = NDElement.real(intArrayOf(10, 10)) { (it[0] + it[1]).toDouble() }
         assertEquals(ndArray[5, 5], 10.0)
     }
 }
\ No newline at end of file
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
index 3bbd620fa..3c4160329 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
@@ -1,7 +1,7 @@
 package scientifik.kmath.structures
 
 import scientifik.kmath.operations.Norm
-import scientifik.kmath.structures.NDElements.real2D
+import scientifik.kmath.structures.NDElement.Companion.real2D
 import kotlin.math.abs
 import kotlin.math.pow
 import kotlin.test.Test
diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt
index 68e989e81..35ee5e6b7 100644
--- a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt
+++ b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt
@@ -23,3 +23,4 @@ object ComplexBufferSpec : FixedSizeBufferSpec<Complex> {
  */
 fun Complex.Companion.createBuffer(size: Int) = ObjectBuffer.create(ComplexBufferSpec, size)
 
+