diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt
index bdc725e63..1ac413818 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt
@@ -14,7 +14,7 @@ import kotlin.math.min
  *
  * @param coefficients constant is the leftmost coefficient.
  */
-public data class Polynomial<T>(public val coefficients: List<T>) : AbstractPolynomial<T> {
+public data class Polynomial<C>(public val coefficients: List<C>) : AbstractPolynomial<C> {
     override fun toString(): String = "Polynomial$coefficients"
 }
 
@@ -44,7 +44,7 @@ internal fun polynomialError(message: Any): Nothing = throw PolynomialError(mess
  * [reverse] parameter is true.
  */
 @Suppress("FunctionName")
-public fun <T> Polynomial(coefficients: List<T>, reverse: Boolean = false): Polynomial<T> =
+public fun <C> Polynomial(coefficients: List<C>, reverse: Boolean = false): Polynomial<C> =
     Polynomial(with(coefficients) { if (reverse) reversed() else this })
 
 /**
@@ -52,10 +52,10 @@ public fun <T> Polynomial(coefficients: List<T>, reverse: Boolean = false): Poly
  * [reverse] parameter is true.
  */
 @Suppress("FunctionName")
-public fun <T> Polynomial(vararg coefficients: T, reverse: Boolean = false): Polynomial<T> =
+public fun <C> Polynomial(vararg coefficients: C, reverse: Boolean = false): Polynomial<C> =
     Polynomial(with(coefficients) { if (reverse) reversed() else toList() })
 
-public fun <T> T.asPolynomial() : Polynomial<T> = Polynomial(listOf(this))
+public fun <C> C.asPolynomial() : Polynomial<C> = Polynomial(listOf(this))
 
 // endregion
 
@@ -352,7 +352,7 @@ public open class PolynomialSpace<C, A : Ring<C>>(
                 Polynomial(
                     (0..(thisDegree + otherDegree))
                         .map { d ->
-                            (max(0, d - otherDegree)..(min(thisDegree, d)))
+                            (max(0, d - otherDegree)..min(thisDegree, d))
                                 .map { coefficients[it] * other.coefficients[d - it] }
                                 .reduce { acc, rational -> acc + rational }
                         }
@@ -370,16 +370,14 @@ public open class PolynomialSpace<C, A : Ring<C>>(
      */
     public override fun Polynomial<C>.isOne(): Boolean =
         with(coefficients) {
-            isNotEmpty() &&
-                    asSequence().withIndex().any { (index, c) -> if (index == 0) c.isOne() else c.isZero() } // TODO: It's better to write new methods like `anyIndexed`. But what's better way to do it?
+            isNotEmpty() && withIndex().any { (index, c) -> if (index == 0) c.isOne() else c.isZero() }
         }
     /**
      * Check if the instant is minus unit polynomial.
      */
     public override fun Polynomial<C>.isMinusOne(): Boolean =
         with(coefficients) {
-            isNotEmpty() &&
-                    asSequence().withIndex().any { (index, c) -> if (index == 0) c.isMinusOne() else c.isZero() } // TODO: It's better to write new methods like `anyIndexed`. But what's better way to do it?
+            isNotEmpty() && withIndex().any { (index, c) -> if (index == 0) c.isMinusOne() else c.isZero() }
         }
 
     /**
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt
index 7e40b6246..a41a6cbd6 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt
@@ -9,7 +9,8 @@ import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.*
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract
-import kotlin.jvm.JvmName
+import kotlin.math.max
+import kotlin.math.min
 import kotlin.math.pow
 
 
@@ -77,15 +78,32 @@ public fun <C> Polynomial<C>.substitute(ring: Ring<C>, arg: C): C = ring {
     return result
 }
 
-// TODO: (Waiting for hero) Replace with optimisation: the [result] may be unboxed, and all operations may be performed
-//  as soon as possible on it
 public fun <C> Polynomial<C>.substitute(ring: Ring<C>, arg: Polynomial<C>) : Polynomial<C> = ring.polynomial {
     if (coefficients.isEmpty()) return zero
-    var result: Polynomial<C> = coefficients.last().asPolynomial()
-    for (j in coefficients.size - 2 downTo 0) {
-        result = (arg * result) + coefficients[j]
+
+    val thisDegree = degree
+    if (thisDegree == -1) return zero
+    val argDegree = arg.degree
+    if (argDegree == -1) return coefficients[0].asPolynomial()
+    val constantZero = constantZero
+    val resultCoefs: MutableList<C> = MutableList(thisDegree + argDegree + 1) { constantZero }
+    val resultCoefsUpdate: MutableList<C> = MutableList(thisDegree + argDegree + 1) { constantZero }
+    var resultDegree = 0
+    for (deg in thisDegree downTo 0) {
+        resultCoefsUpdate[0] = coefficients[deg]
+        for (updateDeg in 0 .. resultDegree + argDegree) {
+            var newC = resultCoefsUpdate[updateDeg]
+            for (deg1 in max(0, updateDeg - argDegree)..min(resultDegree, updateDeg))
+                newC += resultCoefs[deg1] * arg.coefficients[updateDeg - deg1]
+            resultCoefsUpdate[updateDeg] = newC
+        }
+        resultDegree += argDegree
+        for (updateDeg in 0 .. resultDegree + argDegree) {
+            resultCoefs[updateDeg] = resultCoefsUpdate[updateDeg]
+            resultCoefsUpdate[updateDeg] = constantZero
+        }
     }
-    return result
+    return Polynomial<C>(resultCoefs)
 }
 
 /**
@@ -109,7 +127,7 @@ public fun <C, A : Ring<C>> Polynomial<C>.asPolynomialFunctionOver(ring: A): (Po
 public fun <C, A> Polynomial<C>.derivative(
     algebra: A,
 ): Polynomial<C> where  A : Ring<C>, A : NumericAlgebra<C> = algebra {
-    Polynomial(coefficients.drop(1).mapIndexed { index, c -> number(index) * c })
+    Polynomial(coefficients.drop(1).mapIndexed { index, c -> number(index + 1) * c })
 }
 
 /**
@@ -120,8 +138,7 @@ public fun <C, A> Polynomial<C>.nthDerivative(
     algebra: A,
     order: UInt,
 ): Polynomial<C> where  A : Ring<C>, A : NumericAlgebra<C> = algebra {
-    TODO()
-    Polynomial(coefficients.drop(order.toInt()).mapIndexed { index, c ->  number(index) * c })
+    Polynomial(coefficients.drop(order.toInt()).mapIndexed { index, c -> (1..order.toInt()).fold(c) { acc, i -> acc * number(index + i) } })
 }
 
 /**
@@ -133,7 +150,7 @@ public fun <C, A> Polynomial<C>.antiderivative(
 ): Polynomial<C> where  A : Field<C>, A : NumericAlgebra<C> = algebra {
     val integratedCoefficients = buildList(coefficients.size + 1) {
         add(zero)
-        coefficients.forEachIndexed{ index, t -> add(t / number(index + 1)) }
+        coefficients.mapIndexedTo(this) { index, t -> t / number(index + 1) }
     }
     Polynomial(integratedCoefficients)
 }
@@ -146,12 +163,11 @@ public fun <C, A> Polynomial<C>.nthAntiderivative(
     algebra: A,
     order: UInt,
 ): Polynomial<C> where  A : Field<C>, A : NumericAlgebra<C> = algebra {
-    TODO()
-    val integratedCoefficients = buildList(coefficients.size + 1) {
-        add(zero)
-        coefficients.forEachIndexed{ index, t -> add(t / number(index + 1)) }
+    val newCoefficients = buildList(coefficients.size + order.toInt()) {
+        repeat(order.toInt()) { add(zero) }
+        coefficients.mapIndexedTo(this) { index, c -> (1..order.toInt()).fold(c) { acc, i -> acc / number(index + i) } }
     }
-    Polynomial(integratedCoefficients)
+    return Polynomial(newCoefficients)
 }
 
 /**