diff --git a/examples/src/main/kotlin/space/kscience/kmath/functions/polynomials.kt b/examples/src/main/kotlin/space/kscience/kmath/functions/polynomials.kt
new file mode 100644
index 000000000..273fe5cb9
--- /dev/null
+++ b/examples/src/main/kotlin/space/kscience/kmath/functions/polynomials.kt
@@ -0,0 +1,396 @@
+/*
+ * 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.functions
+
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.operations.algebra
+
+
+/**
+ * Shows [ListPolynomial]s' and [ListRationalFunction]s' capabilities.
+ */
+fun listPolynomialsExample() {
+    // [ListPolynomial] is a representation of a univariate polynomial as a list of coefficients from the least term to
+    // the greatest term. For example,
+    val polynomial1: ListPolynomial<Int> = ListPolynomial(listOf(2, -3, 1))
+    // represents polynomial 2 + (-3) x + x^2
+
+    // There are also shortcut fabrics:
+    val polynomial2: ListPolynomial<Int> = ListPolynomial(2, -3, 1)
+    println(polynomial1 == polynomial2) // true
+    // and even
+    val polynomial3: ListPolynomial<Int> = 57.asListPolynomial()
+    val polynomial4: ListPolynomial<Int> = ListPolynomial(listOf(57))
+    println(polynomial3 == polynomial4) // true
+
+    val polynomial5: ListPolynomial<Int> = ListPolynomial(3, -1)
+    // For every ring there can be provided a polynomial ring:
+    Int.algebra.listPolynomialSpace {
+        println(-polynomial5 == ListPolynomial(-3, 1)) // true
+        println(polynomial1 + polynomial5 == ListPolynomial(5, -4, 1)) // true
+        println(polynomial1 - polynomial5 == ListPolynomial(-1, -2, 1)) // true
+        println(polynomial1 * polynomial5 == ListPolynomial(6, -11, 6, -1)) // true
+    }
+    // You can even write
+    val x: ListPolynomial<Double> = ListPolynomial(0.0, 1.0)
+    val polynomial6: ListPolynomial<Double> = ListPolynomial(2.0, -3.0, 1.0)
+    Double.algebra.listPolynomialSpace {
+        println(2 - 3 * x + x * x == polynomial6)
+        println(2.0 - 3.0 * x + x * x == polynomial6)
+    }
+
+    // Also there are some utilities for polynomials:
+    println(polynomial1.substitute(Int.algebra, 1) == 0) // true, because 2 + (-3) * 1 + 1^2 = 0
+    println(polynomial1.substitute(Int.algebra, polynomial5) == polynomial1) // true, because 2 + (-3) * (3-x) + (3-x)^2 = 2 - 3x + x^2
+    println(polynomial1.derivative(Int.algebra) == ListPolynomial(-3, 2)) // true, (2 - 3x + x^2)' = -3 + 2x
+    println(polynomial1.nthDerivative(Int.algebra, 2) == 2.asListPolynomial()) // true, (2 - 3x + x^2)'' = 2
+
+    // Lastly, there are rational functions and some other utilities:
+    Double.algebra.listRationalFunctionSpace {
+        val rationalFunction1: ListRationalFunction<Double> = ListRationalFunction(listOf(2.0, -3.0, 1.0), listOf(3.0, -1.0))
+        // It's just (2 - 3x + x^2)/(3 - x)
+
+        val rationalFunction2 : ListRationalFunction<Double> = ListRationalFunction(listOf(5.0, -4.0, 1.0), listOf(3.0, -1.0))
+        // It's just (5 - 4x + x^2)/(3 - x)
+
+        println(rationalFunction1 + 1 == rationalFunction2)
+    }
+}
+
+/**
+ * Shows [NumberedPolynomial]s' and [NumberedRationalFunction]s' capabilities.
+ */
+fun numberedPolynomialsExample() {
+    // Consider polynomial
+    // 3 + 5 x_2 - 7 x_1^2 x_3
+    // Consider, for example, its term -7 x_1^2 x_3. -7 is a coefficient of the term, whereas (2, 0, 1, 0, 0, ...) is
+    // description of degrees of variables x_1, x_2, ... in the term. Such description with removed leading zeros
+    // [2, 0, 1] is called "signature" of the term -7 x_1^2 x_3.
+
+    val polynomial1: NumberedPolynomial<Int>
+    with(Int.algebra) {
+        // [NumberedPolynomial] is a representation of a multivariate polynomial, that stores terms in a map with terms'
+        // signatures as the map's keys and terms' coefficients as corresponding values. For example,
+        polynomial1 = NumberedPolynomial(
+            mapOf(
+                listOf<UInt>() to 3,
+                listOf(0u, 1u) to 5,
+                listOf(2u, 0u, 1u) to -7,
+            )
+        )
+        // represents polynomial 3 + 5 x_2 - 7 x_1^2 x_3
+
+        // This `NumberedPolynomial` function needs context of either ring of constant (as `Int.algebra` in this example)
+        // or space of NumberedPolynomials over it. To understand why it is like this see documentations of functions
+        // NumberedPolynomial and NumberedPolynomialWithoutCheck
+
+        // There are also shortcut fabrics:
+        val polynomial2: NumberedPolynomial<Int> = NumberedPolynomial(
+            listOf<UInt>() to 3,
+            listOf(0u, 1u) to 5,
+            listOf(2u, 0u, 1u) to -7,
+        )
+        println(polynomial1 == polynomial2) // true
+        // and even
+        val polynomial3: NumberedPolynomial<Int> = 57.asNumberedPolynomial() // This one actually does not algebraic context!
+        val polynomial4: NumberedPolynomial<Int> = NumberedPolynomial(listOf<UInt>() to 57)
+        println(polynomial3 == polynomial4) // true
+
+        numberedPolynomialSpace {
+            // Also there is DSL for constructing NumberedPolynomials:
+            val polynomial5: NumberedPolynomial<Int> = NumberedPolynomial {
+                3 {}
+                5 { 2 inPowerOf 1u }
+                -7 with { 1 pow 2u; 3 pow 1u }
+                // `pow` and `inPowerOf` are the same
+                // `with` is omittable
+            }
+            println(polynomial1 == polynomial5) // true
+
+            // Unfortunately the DSL does not work good in bare context of constants' ring, so for now it's disabled and
+            // works only in NumberedPolynomialSpace and NumberedRationalFunctionSpace
+        }
+    }
+
+    val polynomial6: NumberedPolynomial<Int> = with(Int.algebra) {
+        NumberedPolynomial(
+            listOf<UInt>() to 7,
+            listOf(0u, 1u) to -5,
+            listOf(2u, 0u, 1u) to 0,
+            listOf(0u, 0u, 0u, 4u) to 4,
+        )
+    }
+    // For every ring there can be provided a polynomial ring:
+    Int.algebra.numberedPolynomialSpace {
+        println(
+            -polynomial6 == NumberedPolynomial {
+                (-7) {}
+                5 { 2 pow 1u }
+                0 { 1 pow 2u; 3 pow 1u }
+                (-4) { 4 pow 4u }
+            }
+        ) // true
+        println(
+            polynomial1 + polynomial6 == NumberedPolynomial {
+                10 {}
+                0 { 2 pow 1u }
+                (-7) { 1 pow 2u; 3 pow 1u }
+                4 { 4 pow 4u }
+            }
+        ) // true
+        println(
+            polynomial1 - polynomial6 == NumberedPolynomial {
+                (-4) {}
+                10 { 2 pow 1u }
+                (-7) { 1 pow 2u; 3 pow 1u }
+                (-4) { 4 pow 4u }
+            }
+        ) // true
+
+        polynomial1 * polynomial6 // Multiplication works too
+    }
+
+    Double.algebra.numberedPolynomialSpace {
+        // You can even write
+        val x_1: NumberedPolynomial<Double> = NumberedPolynomial { 1.0 { 1 pow 1u } }
+        val x_2: NumberedPolynomial<Double> = NumberedPolynomial { 1.0 { 2 pow 1u } }
+        val x_3: NumberedPolynomial<Double> = NumberedPolynomial { 1.0 { 3 pow 1u } }
+        val polynomial7: NumberedPolynomial<Double> = NumberedPolynomial {
+            3.0 {}
+            5.0 { 2 pow 1u }
+            (-7.0) { 1 pow 2u; 3 pow 1u }
+        }
+        Double.algebra.listPolynomialSpace {
+            println(3 + 5 * x_2 - 7 * x_1 * x_1 * x_3 == polynomial7)
+            println(3.0 + 5.0 * x_2 - 7.0 * x_1 * x_1 * x_3 == polynomial7)
+        }
+    }
+
+    Int.algebra.numberedPolynomialSpace {
+        val x_4: NumberedPolynomial<Int> = NumberedPolynomial { 1 { 4 pow 1u } }
+        // Also there are some utilities for polynomials:
+        println(polynomial1.substitute(mapOf(0 to 1, 1 to -2, 2 to -1)) == 0.asNumberedPolynomial()) // true,
+            // because it's substitution x_1 -> 1, x_2 -> -2, x_3 -> -1,
+            // so 3 + 5 x_2 - 7 x_1^2 x_3 = 3 + 5 * (-2) - 7 * 1^2 * (-1) = 3 - 10 + 7 = 0
+        println(
+            polynomial1.substitute(mapOf(1 to x_4)) == NumberedPolynomial {
+                3 {}
+                5 { 4 pow 1u }
+                (-7) { 1 pow 2u; 3 pow 1u }
+            }
+        ) // true, because it's substitution x_2 -> x_4, so result is 3 + 5 x_4 - 7 x_1^2 x_3
+        println(
+            polynomial1.derivativeWithRespectTo(Int.algebra, 1) ==
+                    NumberedPolynomial { 5 {} }
+        ) // true, d/dx_2 (3 + 5 x_2 - 7 x_1^2 x_3) = 5
+    }
+
+    // Lastly, there are rational functions and some other utilities:
+    Double.algebra.numberedRationalFunctionSpace {
+        val rationalFunction1: NumberedRationalFunction<Double> = NumberedRationalFunction(
+            NumberedPolynomial {
+                2.0 {}
+                (-3.0) { 1 pow 1u }
+                1.0 { 1 pow 2u }
+            },
+            NumberedPolynomial {
+                3.0 {}
+                (-1.0) { 1 pow 1u }
+            }
+        )
+        // It's just (2 - 3x + x^2)/(3 - x) where x = x_1
+
+        val rationalFunction2: NumberedRationalFunction<Double> = NumberedRationalFunction(
+            NumberedPolynomial {
+                5.0 {}
+                (-4.0) { 1 pow 1u }
+                1.0 { 1 pow 2u }
+            },
+            NumberedPolynomial {
+                3.0 {}
+                (-1.0) { 1 pow 1u }
+            }
+        )
+        // It's just (5 - 4x + x^2)/(3 - x) where x = x_1
+
+        println(rationalFunction1 + 1 == rationalFunction2)
+    }
+}
+
+/**
+ * Shows [LabeledPolynomial]s' and [LabeledRationalFunction]s' capabilities.
+ */
+fun labeledPolynomialsExample() {
+    val x by symbol
+    val y by symbol
+    val z by symbol
+    val t by symbol
+
+    // Consider polynomial
+    // 3 + 5 y - 7 x^2 z
+    // Consider, for example, its term -7 x^2 z. -7 is a coefficient of the term, whereas matching (x -> 2, z -> 3) is
+    // description of degrees of variables x_1, x_2, ... in the term. Such description is called "signature" of the
+    // term -7 x_1^2 x_3.
+
+    val polynomial1: LabeledPolynomial<Int>
+    with(Int.algebra) {
+        // [LabeledPolynomial] is a representation of a multivariate polynomial, that stores terms in a map with terms'
+        // signatures as the map's keys and terms' coefficients as corresponding values. For example,
+        polynomial1 = LabeledPolynomial(
+            mapOf(
+                mapOf<Symbol, UInt>() to 3,
+                mapOf(y to 1u) to 5,
+                mapOf(x to 2u, z to 1u) to -7,
+            )
+        )
+        // represents polynomial 3 + 5 y - 7 x^2 z
+
+        // This `LabeledPolynomial` function needs context of either ring of constant (as `Int.algebra` in this example)
+        // or space of LabeledPolynomials over it. To understand why it is like this see documentations of functions
+        // LabeledPolynomial and LabeledPolynomialWithoutCheck
+
+        // There are also shortcut fabrics:
+        val polynomial2: LabeledPolynomial<Int> = LabeledPolynomial(
+            mapOf<Symbol, UInt>() to 3,
+            mapOf(y to 1u) to 5,
+            mapOf(x to 2u, z to 1u) to -7,
+        )
+        println(polynomial1 == polynomial2) // true
+        // and even
+        val polynomial3: LabeledPolynomial<Int> = 57.asLabeledPolynomial() // This one actually does not algebraic context!
+        val polynomial4: LabeledPolynomial<Int> = LabeledPolynomial(mapOf<Symbol, UInt>() to 57)
+        println(polynomial3 == polynomial4) // true
+
+        labeledPolynomialSpace {
+            // Also there is DSL for constructing NumberedPolynomials:
+            val polynomial5: LabeledPolynomial<Int> = LabeledPolynomial {
+                3 {}
+                5 { y inPowerOf 1u }
+                -7 with { x pow 2u; z pow 1u }
+                // `pow` and `inPowerOf` are the same
+                // `with` is omittable
+            }
+            println(polynomial1 == polynomial5) // true
+
+            // Unfortunately the DSL does not work good in bare context of constants' ring, so for now it's disabled and
+            // works only in NumberedPolynomialSpace and NumberedRationalFunctionSpace
+        }
+    }
+
+    val polynomial6: LabeledPolynomial<Int> = with(Int.algebra) {
+        LabeledPolynomial(
+            mapOf<Symbol, UInt>() to 7,
+            mapOf(y to 1u) to -5,
+            mapOf(x to 2u, z to 1u) to 0,
+            mapOf(t to 4u) to 4,
+        )
+    }
+    // For every ring there can be provided a polynomial ring:
+    Int.algebra.labeledPolynomialSpace {
+        println(
+            -polynomial6 == LabeledPolynomial {
+                (-7) {}
+                5 { y pow 1u }
+                0 { x pow 2u; z pow 1u }
+                (-4) { t pow 4u }
+            }
+        ) // true
+        println(
+            polynomial1 + polynomial6 == LabeledPolynomial {
+                10 {}
+                0 { y pow 1u }
+                (-7) { x pow 2u; z pow 1u }
+                4 { t pow 4u }
+            }
+        ) // true
+        println(
+            polynomial1 - polynomial6 == LabeledPolynomial {
+                (-4) {}
+                10 { y pow 1u }
+                (-7) { x pow 2u; z pow 1u }
+                (-4) { t pow 4u }
+            }
+        ) // true
+
+        polynomial1 * polynomial6 // Multiplication works too
+    }
+
+    Double.algebra.labeledPolynomialSpace {
+        // You can even write
+        val polynomial7: LabeledPolynomial<Double> = LabeledPolynomial {
+            3.0 {}
+            5.0 { y pow 1u }
+            (-7.0) { x pow 2u; z pow 1u }
+        }
+        Double.algebra.listPolynomialSpace {
+            println(3 + 5 * y - 7 * x * x * z == polynomial7)
+            println(3.0 + 5.0 * y - 7.0 * x * x * z == polynomial7)
+        }
+    }
+
+    Int.algebra.labeledPolynomialSpace {
+        // Also there are some utilities for polynomials:
+        println(polynomial1.substitute(mapOf(x to 1, y to -2, z to -1)) == 0.asLabeledPolynomial()) // true,
+        // because it's substitution x -> 1, y -> -2, z -> -1,
+        // so 3 + 5 y - 7 x^2 z = 3 + 5 * (-2) - 7 * 1^2 * (-1) = 3 - 10 + 7 = 0
+        println(
+            polynomial1.substitute(mapOf(y to t.asPolynomial())) == LabeledPolynomial {
+                3 {}
+                5 { t pow 1u }
+                (-7) { x pow 2u; z pow 1u }
+            }
+        ) // true, because it's substitution y -> t, so result is 3 + 5 t - 7 x^2 z
+        println(
+            polynomial1.derivativeWithRespectTo(Int.algebra, y) == LabeledPolynomial { 5 {} }
+        ) // true, d/dy (3 + 5 y - 7 x^2 z) = 5
+    }
+
+    // Lastly, there are rational functions and some other utilities:
+    Double.algebra.labeledRationalFunctionSpace {
+        val rationalFunction1: LabeledRationalFunction<Double> = LabeledRationalFunction(
+            LabeledPolynomial {
+                2.0 {}
+                (-3.0) { x pow 1u }
+                1.0 { x pow 2u }
+            },
+            LabeledPolynomial {
+                3.0 {}
+                (-1.0) { x pow 1u }
+            }
+        )
+        // It's just (2 - 3x + x^2)/(3 - x)
+
+        val rationalFunction2: LabeledRationalFunction<Double> = LabeledRationalFunction(
+            LabeledPolynomial {
+                5.0 {}
+                (-4.0) { x pow 1u }
+                1.0 { x pow 2u }
+            },
+            LabeledPolynomial {
+                3.0 {}
+                (-1.0) { x pow 1u }
+            }
+        )
+        // It's just (5 - 4x + x^2)/(3 - x)
+
+        println(rationalFunction1 + 1 == rationalFunction2)
+    }
+}
+
+fun main() {
+    println("ListPolynomials:")
+    listPolynomialsExample()
+    println()
+
+    println("NumberedPolynomials:")
+    numberedPolynomialsExample()
+    println()
+
+    println("ListPolynomials:")
+    labeledPolynomialsExample()
+    println()
+}
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt
index b0c54502d..061325b60 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt
@@ -522,7 +522,7 @@ public class LabeledPolynomialSpace<C, A : Ring<C>>(
     override operator fun LabeledPolynomial<C>.plus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
         LabeledPolynomialAsIs(
             buildMap(coefficients.size + other.coefficients.size) {
-                other.coefficients.mapValuesTo(this) { it.value }
+                coefficients.mapValuesTo(this) { it.value }
                 other.coefficients.mapValuesTo(this) { (key, value) -> if (key in this) this[key]!! + value else value }
             }
         )
@@ -532,7 +532,7 @@ public class LabeledPolynomialSpace<C, A : Ring<C>>(
     override operator fun LabeledPolynomial<C>.minus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
         LabeledPolynomialAsIs(
             buildMap(coefficients.size + other.coefficients.size) {
-                other.coefficients.mapValuesTo(this) { it.value }
+                coefficients.mapValuesTo(this) { it.value }
                 other.coefficients.mapValuesTo(this) { (key, value) -> if (key in this) this[key]!! - value else -value }
             }
         )
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledRationalFunction.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledRationalFunction.kt
index 03f323813..1b0bbbec8 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledRationalFunction.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledRationalFunction.kt
@@ -13,7 +13,7 @@ import kotlin.jvm.JvmName
 /**
  * Represents multivariate rational function that stores its numerator and denominator as [LabeledPolynomial]s.
  */
-public class LabeledRationalFunction<C>(
+public data class LabeledRationalFunction<C>(
     public override val numerator: LabeledPolynomial<C>,
     public override val denominator: LabeledPolynomial<C>
 ) : RationalFunction<C, LabeledPolynomial<C>> {
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/NumberedRationalFunction.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/NumberedRationalFunction.kt
index 92f507735..497755126 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/NumberedRationalFunction.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/NumberedRationalFunction.kt
@@ -15,7 +15,7 @@ import kotlin.math.max
 /**
  * Represents multivariate rational function that stores its numerator and denominator as [NumberedPolynomial]s.
  */
-public class NumberedRationalFunction<C> internal constructor(
+public data class NumberedRationalFunction<C>(
     public override val numerator: NumberedPolynomial<C>,
     public override val denominator: NumberedPolynomial<C>
 ) : RationalFunction<C, NumberedPolynomial<C>> {
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledConstructors.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledConstructors.kt
index 47325c4bb..7586b7012 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledConstructors.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledConstructors.kt
@@ -73,6 +73,8 @@ public inline fun <C> LabeledPolynomialWithoutCheck(vararg pairs: Pair<Map<Symbo
  * 1. Zeros at the ends of terms' signatures (e.g. [coefs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public fun <C> LabeledPolynomial(coefs: Map<Map<Symbol, UInt>, C>, add: (C, C) -> C) : LabeledPolynomial<C> {
     val fixedCoefs = mutableMapOf<Map<Symbol, UInt>, C>()
@@ -93,6 +95,8 @@ public fun <C> LabeledPolynomial(coefs: Map<Map<Symbol, UInt>, C>, add: (C, C) -
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public fun <C> LabeledPolynomial(pairs: Collection<Pair<Map<Symbol, UInt>, C>>, add: (C, C) -> C) : LabeledPolynomial<C> {
     val fixedCoefs = mutableMapOf<Map<Symbol, UInt>, C>()
@@ -113,6 +117,8 @@ public fun <C> LabeledPolynomial(pairs: Collection<Pair<Map<Symbol, UInt>, C>>,
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public fun <C> LabeledPolynomial(vararg pairs: Pair<Map<Symbol, UInt>, C>, add: (C, C) -> C) : LabeledPolynomial<C> {
     val fixedCoefs = mutableMapOf<Map<Symbol, UInt>, C>()
@@ -135,6 +141,8 @@ public fun <C> LabeledPolynomial(vararg pairs: Pair<Map<Symbol, UInt>, C>, add:
  * 1. Zeros at the ends of terms' signatures (e.g. [coefs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> A.LabeledPolynomial(coefs: Map<Map<Symbol, UInt>, C>) : LabeledPolynomial<C> = LabeledPolynomial(coefs, ::add)
 /**
@@ -144,6 +152,8 @@ public inline fun <C, A: Ring<C>> A.LabeledPolynomial(coefs: Map<Map<Symbol, UIn
  * 1. Zeros at the ends of terms' signatures (e.g. [coefs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> LabeledPolynomialSpace<C, A>.LabeledPolynomial(coefs: Map<Map<Symbol, UInt>, C>) : LabeledPolynomial<C> = LabeledPolynomial(coefs) { left: C, right: C -> left + right }
 
@@ -154,6 +164,8 @@ public inline fun <C, A: Ring<C>> LabeledPolynomialSpace<C, A>.LabeledPolynomial
  * 1. Zeros at the ends of terms' signatures (e.g. [coefs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> LabeledRationalFunctionSpace<C, A>.LabeledPolynomial(coefs: Map<Map<Symbol, UInt>, C>) : LabeledPolynomial<C> = LabeledPolynomial(coefs) { left: C, right: C -> left + right }
 
@@ -164,6 +176,8 @@ public inline fun <C, A: Ring<C>> LabeledRationalFunctionSpace<C, A>.LabeledPoly
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> A.LabeledPolynomial(pairs: Collection<Pair<Map<Symbol, UInt>, C>>) : LabeledPolynomial<C> = LabeledPolynomial(pairs, ::add)
 
@@ -174,6 +188,8 @@ public inline fun <C, A: Ring<C>> A.LabeledPolynomial(pairs: Collection<Pair<Map
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> LabeledPolynomialSpace<C, A>.LabeledPolynomial(pairs: Collection<Pair<Map<Symbol, UInt>, C>>) : LabeledPolynomial<C> = LabeledPolynomial(pairs) { left: C, right: C -> left + right }
 /**
@@ -183,6 +199,8 @@ public inline fun <C, A: Ring<C>> LabeledPolynomialSpace<C, A>.LabeledPolynomial
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> LabeledRationalFunctionSpace<C, A>.LabeledPolynomial(pairs: Collection<Pair<Map<Symbol, UInt>, C>>) : LabeledPolynomial<C> = LabeledPolynomial(pairs) { left: C, right: C -> left + right }
 
@@ -193,6 +211,8 @@ public inline fun <C, A: Ring<C>> LabeledRationalFunctionSpace<C, A>.LabeledPoly
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> A.LabeledPolynomial(vararg pairs: Pair<Map<Symbol, UInt>, C>) : LabeledPolynomial<C> = LabeledPolynomial(*pairs) { left: C, right: C -> left + right }
 /**
@@ -202,6 +222,8 @@ public inline fun <C, A: Ring<C>> A.LabeledPolynomial(vararg pairs: Pair<Map<Sym
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> LabeledPolynomialSpace<C, A>.LabeledPolynomial(vararg pairs: Pair<Map<Symbol, UInt>, C>) : LabeledPolynomial<C> = LabeledPolynomial(*pairs) { left: C, right: C -> left + right }
 /**
@@ -211,6 +233,8 @@ public inline fun <C, A: Ring<C>> LabeledPolynomialSpace<C, A>.LabeledPolynomial
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see LabeledPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> LabeledRationalFunctionSpace<C, A>.LabeledPolynomial(vararg pairs: Pair<Map<Symbol, UInt>, C>) : LabeledPolynomial<C> = LabeledPolynomial(*pairs) { left: C, right: C -> left + right }
 
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedConstructors.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedConstructors.kt
index 4850e6cec..3cd398623 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedConstructors.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedConstructors.kt
@@ -72,6 +72,8 @@ public inline fun <C> NumberedPolynomialWithoutCheck(vararg pairs: Pair<List<UIn
  * 1. Zeros at the ends of terms' signatures (e.g. [coefs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public fun <C> NumberedPolynomial(coefs: Map<List<UInt>, C>, add: (C, C) -> C) : NumberedPolynomial<C> {
     val fixedCoefs = mutableMapOf<List<UInt>, C>()
@@ -92,6 +94,8 @@ public fun <C> NumberedPolynomial(coefs: Map<List<UInt>, C>, add: (C, C) -> C) :
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public fun <C> NumberedPolynomial(pairs: Collection<Pair<List<UInt>, C>>, add: (C, C) -> C) : NumberedPolynomial<C> {
     val fixedCoefs = mutableMapOf<List<UInt>, C>()
@@ -112,6 +116,8 @@ public fun <C> NumberedPolynomial(pairs: Collection<Pair<List<UInt>, C>>, add: (
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public fun <C> NumberedPolynomial(vararg pairs: Pair<List<UInt>, C>, add: (C, C) -> C) : NumberedPolynomial<C> {
     val fixedCoefs = mutableMapOf<List<UInt>, C>()
@@ -134,6 +140,8 @@ public fun <C> NumberedPolynomial(vararg pairs: Pair<List<UInt>, C>, add: (C, C)
  * 1. Zeros at the ends of terms' signatures (e.g. [coefs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> A.NumberedPolynomial(coefs: Map<List<UInt>, C>) : NumberedPolynomial<C> = NumberedPolynomial(coefs, ::add)
 /**
@@ -143,6 +151,8 @@ public inline fun <C, A: Ring<C>> A.NumberedPolynomial(coefs: Map<List<UInt>, C>
  * 1. Zeros at the ends of terms' signatures (e.g. [coefs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> NumberedPolynomialSpace<C, A>.NumberedPolynomial(coefs: Map<List<UInt>, C>) : NumberedPolynomial<C> = NumberedPolynomial(coefs) { left: C, right: C -> left + right }
 
@@ -153,6 +163,8 @@ public inline fun <C, A: Ring<C>> NumberedPolynomialSpace<C, A>.NumberedPolynomi
  * 1. Zeros at the ends of terms' signatures (e.g. [coefs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> NumberedRationalFunctionSpace<C, A>.NumberedPolynomial(coefs: Map<List<UInt>, C>) : NumberedPolynomial<C> = NumberedPolynomial(coefs) { left: C, right: C -> left + right }
 
@@ -163,6 +175,8 @@ public inline fun <C, A: Ring<C>> NumberedRationalFunctionSpace<C, A>.NumberedPo
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> A.NumberedPolynomial(pairs: Collection<Pair<List<UInt>, C>>) : NumberedPolynomial<C> = NumberedPolynomial(pairs, ::add)
 
@@ -173,6 +187,8 @@ public inline fun <C, A: Ring<C>> A.NumberedPolynomial(pairs: Collection<Pair<Li
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> NumberedPolynomialSpace<C, A>.NumberedPolynomial(pairs: Collection<Pair<List<UInt>, C>>) : NumberedPolynomial<C> = NumberedPolynomial(pairs) { left: C, right: C -> left + right }
 /**
@@ -192,6 +208,8 @@ public inline fun <C, A: Ring<C>> NumberedRationalFunctionSpace<C, A>.NumberedPo
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> A.NumberedPolynomial(vararg pairs: Pair<List<UInt>, C>) : NumberedPolynomial<C> = NumberedPolynomial(*pairs) { left: C, right: C -> left + right }
 /**
@@ -201,6 +219,8 @@ public inline fun <C, A: Ring<C>> A.NumberedPolynomial(vararg pairs: Pair<List<U
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> NumberedPolynomialSpace<C, A>.NumberedPolynomial(vararg pairs: Pair<List<UInt>, C>) : NumberedPolynomial<C> = NumberedPolynomial(*pairs) { left: C, right: C -> left + right }
 /**
@@ -210,6 +230,8 @@ public inline fun <C, A: Ring<C>> NumberedPolynomialSpace<C, A>.NumberedPolynomi
  * 1. Zeros at the ends of terms' signatures (e.g. [pairs] keys) will be removed. (See [cleanUp].)
  * 1. Terms that happen to have the same signature will be summed up.
  * 1. New map will be formed of resulting terms.
+ *
+ * @see NumberedPolynomialWithoutCheck
  */
 public inline fun <C, A: Ring<C>> NumberedRationalFunctionSpace<C, A>.NumberedPolynomial(vararg pairs: Pair<List<UInt>, C>) : NumberedPolynomial<C> = NumberedPolynomial(*pairs) { left: C, right: C -> left + right }
 
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedUtil.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedUtil.kt
index 06911feca..96f5f480b 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedUtil.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedUtil.kt
@@ -353,7 +353,7 @@ public fun <C, A : Ring<C>> NumberedPolynomial<C>.derivativeWithRespectTo(
         buildMap(coefficients.size) {
             coefficients
                 .forEach { (degs, c) ->
-                    if (degs.size > variable) return@forEach
+                    if (degs.lastIndex < variable) return@forEach
                     put(
                         degs.mapIndexed { index, deg ->
                             when {
@@ -383,7 +383,7 @@ public fun <C, A : Ring<C>> NumberedPolynomial<C>.nthDerivativeWithRespectTo(
         buildMap(coefficients.size) {
             coefficients
                 .forEach { (degs, c) ->
-                    if (degs.size > variable) return@forEach
+                    if (degs.lastIndex < variable) return@forEach
                     put(
                         degs.mapIndexed { index, deg ->
                             when {
@@ -417,7 +417,7 @@ public fun <C, A : Ring<C>> NumberedPolynomial<C>.nthDerivativeWithRespectTo(
         buildMap(coefficients.size) {
             coefficients
                 .forEach { (degs, c) ->
-                    if (degs.size > maxRespectedVariable) return@forEach
+                    if (degs.lastIndex < maxRespectedVariable) return@forEach
                     put(
                         degs.mapIndexed { index, deg ->
                             if (index !in filteredVariablesAndOrders) return@mapIndexed deg