diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledPolynomialUtil.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledPolynomialUtil.kt
index 5689f6e1b..544cca410 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledPolynomialUtil.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledPolynomialUtil.kt
@@ -5,8 +5,10 @@
 
 package space.kscience.kmath.functions
 
+import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.*
 import kotlin.contracts.*
+import kotlin.math.max
 
 
 // TODO: Docs
@@ -321,6 +323,251 @@ public inline fun <C, A : Ring<C>, R> A.labeledPolynomial(block: LabeledPolynomi
 //
 //// endregion
 
-//// region Algebraic derivative and antiderivative
-//// TODO
-//// endregion
\ No newline at end of file
+// region Algebraic derivative and antiderivative
+
+/**
+ * Returns algebraic derivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Ring<C>> LabeledPolynomial<C>.derivativeWithRespectTo(
+    algebra: A,
+    variable: Variable,
+): LabeledPolynomial<C> = algebra {
+    LabeledPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    if (variable !in degs) return@forEach
+                    put(
+                        buildMap {
+                            degs.forEach { (vari, deg) ->
+                                when {
+                                    vari != variable -> put(vari, deg)
+                                    deg > 1u -> put(vari, deg - 1u)
+                                }
+                            }
+                        },
+                        optimizedMultiply(c, degs[variable]!!)
+                    )
+                }
+        }
+    )
+}
+
+/**
+ * Returns algebraic derivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Ring<C>> LabeledPolynomial<C>.derivativeWithRespectTo(
+    algebra: A,
+    variables: Collection<Variable>,
+): LabeledPolynomial<C> = algebra {
+    val cleanedVariables = variables.toSet()
+    if (cleanedVariables.isEmpty()) return this@derivativeWithRespectTo
+    LabeledPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    if (!degs.keys.containsAll(cleanedVariables)) return@forEach
+                    put(
+                        buildMap {
+                            degs.forEach { (vari, deg) ->
+                                when {
+                                    vari !in cleanedVariables -> put(vari, deg)
+                                    deg > 1u -> put(vari, deg - 1u)
+                                }
+                            }
+                        },
+                        cleanedVariables.fold(c) { acc, variable -> optimizedMultiply(acc, degs[variable]!!) }
+                    )
+                }
+        }
+    )
+}
+
+/**
+ * Returns algebraic derivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Ring<C>> LabeledPolynomial<C>.nthDerivativeWithRespectTo(
+    algebra: A,
+    variable: Variable,
+    order: UInt
+): LabeledPolynomial<C> = algebra {
+    if (order == 0u) return this@nthDerivativeWithRespectTo
+    LabeledPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    if (degs.getOrElse(variable) { 0u } < order) return@forEach
+                    put(
+                        buildMap {
+                            degs.forEach { (vari, deg) ->
+                                when {
+                                    vari != variable -> put(vari, deg)
+                                    deg > order -> put(vari, deg - order)
+                                }
+                            }
+                        },
+                        degs[variable]!!.let { deg ->
+                            (deg downTo deg - order + 1u)
+                                .fold(c) { acc, ord -> optimizedMultiply(acc, ord) }
+                        }
+                    )
+                }
+        }
+    )
+}
+
+/**
+ * Returns algebraic derivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Ring<C>> LabeledPolynomial<C>.nthDerivativeWithRespectTo(
+    algebra: A,
+    variablesAndOrders: Map<Variable, UInt>,
+): LabeledPolynomial<C> = algebra {
+    val filteredVariablesAndOrders = variablesAndOrders.filterValues { it != 0u }
+    if (filteredVariablesAndOrders.isEmpty()) return this@nthDerivativeWithRespectTo
+    LabeledPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    if (filteredVariablesAndOrders.any { (variable, order) -> degs.getOrElse(variable) { 0u } < order }) return@forEach
+                    put(
+                        buildMap {
+                            degs.forEach { (vari, deg) ->
+                                if (vari !in filteredVariablesAndOrders) put(vari, deg)
+                                else {
+                                    val order = filteredVariablesAndOrders[vari]!!
+                                    if (deg > order) put(vari, deg - order)
+                                }
+                            }
+                        },
+                        filteredVariablesAndOrders.entries.fold(c) { acc1, (index, order) ->
+                            degs[index]!!.let { deg ->
+                                (deg downTo deg - order + 1u)
+                                    .fold(acc1) { acc2, ord -> optimizedMultiply(acc2, ord) }
+                            }
+                        }
+                    )
+                }
+        }
+    )
+}
+
+/**
+ * Returns algebraic antiderivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Field<C>> LabeledPolynomial<C>.antiderivativeWithRespectTo(
+    algebra: A,
+    variable: Variable,
+): LabeledPolynomial<C> = algebra {
+    LabeledPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    val newDegs = buildMap<Variable, UInt>(degs.size + 1) {
+                        put(variable, 1u)
+                        for ((vari, deg) in degs) put(vari, deg + getOrElse(vari) { 0u })
+                    }
+                    put(
+                        newDegs,
+                        c / optimizedMultiply(one, newDegs[variable]!!)
+                    )
+                }
+        }
+    )
+}
+
+/**
+ * Returns algebraic antiderivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Field<C>> LabeledPolynomial<C>.antiderivativeWithRespectTo(
+    algebra: A,
+    variables: Collection<Variable>,
+): LabeledPolynomial<C> = algebra {
+    val cleanedVariables = variables.toSet()
+    if (cleanedVariables.isEmpty()) return this@antiderivativeWithRespectTo
+    LabeledPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    val newDegs = buildMap<Variable, UInt>(degs.size + 1) {
+                        for (variable in cleanedVariables) put(variable, 1u)
+                        for ((vari, deg) in degs) put(vari, deg + getOrElse(vari) { 0u })
+                    }
+                    put(
+                        newDegs,
+                        cleanedVariables.fold(c) { acc, variable -> acc / optimizedMultiply(one, newDegs[variable]!!) }
+                    )
+                }
+        }
+    )
+}
+
+/**
+ * Returns algebraic derivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Field<C>> LabeledPolynomial<C>.nthAntiderivativeWithRespectTo(
+    algebra: A,
+    variable: Variable,
+    order: UInt
+): LabeledPolynomial<C> = algebra {
+    if (order == 0u) return this@nthAntiderivativeWithRespectTo
+    LabeledPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    val newDegs = buildMap<Variable, UInt>(degs.size + 1) {
+                        put(variable, order)
+                        for ((vari, deg) in degs) put(vari, deg + getOrElse(vari) { 0u })
+                    }
+                    put(
+                        newDegs,
+                        newDegs[variable]!!.let { deg ->
+                            (deg downTo  deg - order + 1u)
+                                .fold(c) { acc, ord -> acc / optimizedMultiply(one, ord) }
+                        }
+                    )
+                }
+        }
+    )
+}
+
+/**
+ * Returns algebraic derivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Field<C>> LabeledPolynomial<C>.nthAntiderivativeWithRespectTo(
+    algebra: A,
+    variablesAndOrders: Map<Variable, UInt>,
+): LabeledPolynomial<C> = algebra {
+    val filteredVariablesAndOrders = variablesAndOrders.filterValues { it != 0u }
+    if (filteredVariablesAndOrders.isEmpty()) return this@nthAntiderivativeWithRespectTo
+    LabeledPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    val newDegs = buildMap<Variable, UInt>(degs.size + 1) {
+                        for ((variable, order) in filteredVariablesAndOrders) put(variable, order)
+                        for ((vari, deg) in degs) put(vari, deg + getOrElse(vari) { 0u })
+                    }
+                    put(
+                        newDegs,
+                        filteredVariablesAndOrders.entries.fold(c) { acc1, (index, order) ->
+                            newDegs[index]!!.let { deg ->
+                                (deg downTo deg - order + 1u)
+                                    .fold(acc1) { acc2, ord -> acc2 / optimizedMultiply(one, ord) }
+                            }
+                        }
+                    )
+                }
+        }
+    )
+}
+
+// endregion
\ No newline at end of file
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedPolynomialUtil.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedPolynomialUtil.kt
index de95d0151..47644b26a 100644
--- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedPolynomialUtil.kt
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedPolynomialUtil.kt
@@ -4,7 +4,7 @@ import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.*
 import kotlin.contracts.*
 import kotlin.jvm.JvmName
-
+import kotlin.math.max
 
 
 // TODO: Docs
@@ -382,70 +382,51 @@ public fun <C, A : Ring<C>> NumberedPolynomial<C>.derivativeWithRespectTo(
 ): NumberedPolynomial<C> = algebra {
     NumberedPolynomial<C>(
         buildMap(coefficients.size) {
-            coefficients.forEach { (degs, c) ->
-                put(
-                    degs.mapIndexed { index, deg ->
-                        when {
-                            index != variable -> deg
-                            deg > 0u -> deg - 1u
-                            else -> return@forEach
-                        }
-                    },
-                    optimizedMultiply(c, degs.getOrElse(variable) { 1u }.toInt())
-                )
-            }
+            coefficients
+                .forEach { (degs, c) ->
+                    if (degs.size > variable) return@forEach
+                    put(
+                        degs.mapIndexed { index, deg ->
+                            when {
+                                index != variable -> deg
+                                deg > 0u -> deg - 1u
+                                else -> return@forEach
+                            }
+                        }.cleanUp(),
+                        optimizedMultiply(c, degs[variable])
+                    )
+                }
         }
     )
 }
 
 /**
  * Returns algebraic derivative of received polynomial.
- */ // TODO: This one does not work!!!
-@UnstableKMathAPI
-public fun <C, A : Ring<C>> NumberedPolynomial<C>.derivativeWithRespectTo(
-    algebra: A,
-    variables: IntArray,
-): NumberedPolynomial<C> = algebra {
-    NumberedPolynomial<C>(
-        buildMap(coefficients.size) {
-            coefficients.forEach { (degs, c) ->
-                put(
-                    degs.mapIndexed { index, deg ->
-                        when {
-                            index !in variables -> deg
-                            deg > 0u -> deg - 1u
-                            else -> return@forEach
-                        }
-                    },
-                    optimizedMultiply(c, variables.fold(1u) { acc, variable -> acc * degs.getOrElse(variable) { 1u } }.toInt())
-                )
-            }
-        }
-    )
-}
-
-/**
- * Returns algebraic derivative of received polynomial.
- */ // TODO: This one does not work!!!
+ */
 @UnstableKMathAPI
 public fun <C, A : Ring<C>> NumberedPolynomial<C>.derivativeWithRespectTo(
     algebra: A,
     variables: Collection<Int>,
 ): NumberedPolynomial<C> = algebra {
+    val cleanedVariables = variables.toSet()
+    if (cleanedVariables.isEmpty()) return this@derivativeWithRespectTo
+    val maxRespectedVariable = cleanedVariables.maxOrNull()!!
     NumberedPolynomial<C>(
         buildMap(coefficients.size) {
-            coefficients.forEach { (degs, c) ->
-                put(
-                    degs.mapIndexed { index, deg ->
-                        when {
-                            index !in variables -> deg
-                            deg > 0u -> deg - 1u
-                            else -> return@forEach
-                        }
-                    },
-                    optimizedMultiply(c, variables.fold(1u) { acc, variable -> acc * degs.getOrElse(variable) { 1u } }.toInt())
-                )
-            }
+            coefficients
+                .forEach { (degs, c) ->
+                    if (degs.size > maxRespectedVariable) return@forEach
+                    put(
+                        degs.mapIndexed { index, deg ->
+                            when {
+                                index !in cleanedVariables -> deg
+                                deg > 0u -> deg - 1u
+                                else -> return@forEach
+                            }
+                        }.cleanUp(),
+                        cleanedVariables.fold(c) { acc, variable -> optimizedMultiply(acc, degs[variable]) }
+                    )
+                }
         }
     )
 }
@@ -459,24 +440,60 @@ public fun <C, A : Ring<C>> NumberedPolynomial<C>.nthDerivativeWithRespectTo(
     variable: Int,
     order: UInt
 ): NumberedPolynomial<C> = algebra {
+    if (order == 0u) return this@nthDerivativeWithRespectTo
     NumberedPolynomial<C>(
         buildMap(coefficients.size) {
-            coefficients.forEach { (degs, c) ->
-                put(
-                    degs.mapIndexed { index, deg ->
-                        when {
-                            index != variable -> deg
-                            deg >= order -> deg - order
-                            else -> return@forEach
+            coefficients
+                .forEach { (degs, c) ->
+                    if (degs.size > variable) return@forEach
+                    put(
+                        degs.mapIndexed { index, deg ->
+                            when {
+                                index != variable -> deg
+                                deg >= order -> deg - order
+                                else -> return@forEach
+                            }
+                        }.cleanUp(),
+                        degs[variable].let { deg ->
+                            (deg downTo deg - order + 1u)
+                                .fold(c) { acc, ord -> optimizedMultiply(acc, ord) }
                         }
-                    },
-                    degs.getOrElse(variable) { 1u }.toInt().let {
-                        (0u until order).fold(c) { acc, ord ->
-                            optimizedMultiply(acc, ord.toInt())
+                    )
+                }
+        }
+    )
+}
+
+/**
+ * Returns algebraic derivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Ring<C>> NumberedPolynomial<C>.nthDerivativeWithRespectTo(
+    algebra: A,
+    variablesAndOrders: Map<Int, UInt>,
+): NumberedPolynomial<C> = algebra {
+    val filteredVariablesAndOrders = variablesAndOrders.filterValues { it != 0u }
+    if (filteredVariablesAndOrders.isEmpty()) return this@nthDerivativeWithRespectTo
+    val maxRespectedVariable = filteredVariablesAndOrders.keys.maxOrNull()!!
+    NumberedPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    if (degs.size > maxRespectedVariable) return@forEach
+                    put(
+                        degs.mapIndexed { index, deg ->
+                            if (index !in filteredVariablesAndOrders) return@mapIndexed deg
+                            val order = filteredVariablesAndOrders[index]!!
+                            if (deg >= order) deg - order else return@forEach
+                        }.cleanUp(),
+                        filteredVariablesAndOrders.entries.fold(c) { acc1, (index, order) ->
+                            degs[index].let { deg ->
+                                (deg downTo deg - order + 1u)
+                                    .fold(acc1) { acc2, ord -> optimizedMultiply(acc2, ord) }
+                            }
                         }
-                    }
-                )
-            }
+                    )
+                }
         }
     )
 }
@@ -491,52 +508,92 @@ public fun <C, A : Field<C>> NumberedPolynomial<C>.antiderivativeWithRespectTo(
 ): NumberedPolynomial<C> = algebra {
     NumberedPolynomial<C>(
         buildMap(coefficients.size) {
-            coefficients.forEach { (degs, c) ->
-                put(
-                    degs.mapIndexed { index, deg -> if(index != variable) deg else deg + 1u },
-                    c / optimizedMultiply(one, degs.getOrElse(variable) { 1u }.toInt())
-                )
-            }
+            coefficients
+                .forEach { (degs, c) ->
+                    put(
+                        List(max(variable + 1, degs.size)) { if (it != variable) degs[it] else degs[it] + 1u },
+                        c / optimizedMultiply(one, degs[variable])
+                    )
+                }
         }
     )
 }
 
 /**
  * Returns algebraic antiderivative of received polynomial.
- */ // TODO: This one does not work!!!
-@UnstableKMathAPI
-public fun <C, A : Field<C>> NumberedPolynomial<C>.antiderivativeWithRespectTo(
-    algebra: A,
-    variables: IntArray,
-): NumberedPolynomial<C> = algebra {
-    NumberedPolynomial<C>(
-        buildMap(coefficients.size) {
-            coefficients.forEach { (degs, c) ->
-                put(
-                    degs.mapIndexed { index, deg -> if(index !in variables) deg else deg + 1u },
-                    c / optimizedMultiply(one, variables.fold(1u) { acc, variable -> acc * degs.getOrElse(variable) { 1u } }.toInt())
-                )
-            }
-        }
-    )
-}
-
-/**
- * Returns algebraic antiderivative of received polynomial.
- */ // TODO: This one does not work!!!
+ */
 @UnstableKMathAPI
 public fun <C, A : Field<C>> NumberedPolynomial<C>.antiderivativeWithRespectTo(
     algebra: A,
     variables: Collection<Int>,
 ): NumberedPolynomial<C> = algebra {
+    val cleanedVariables = variables.toSet()
+    if (cleanedVariables.isEmpty()) return this@antiderivativeWithRespectTo
+    val maxRespectedVariable = cleanedVariables.maxOrNull()!!
     NumberedPolynomial<C>(
         buildMap(coefficients.size) {
-            coefficients.forEach { (degs, c) ->
-                put(
-                    degs.mapIndexed { index, deg -> if(index !in variables) deg else deg + 1u },
-                    c / optimizedMultiply(one, variables.fold(1u) { acc, variable -> acc * degs.getOrElse(variable) { 1u } }.toInt())
-                )
-            }
+            coefficients
+                .forEach { (degs, c) ->
+                    put(
+                        List(max(maxRespectedVariable + 1, degs.size)) { if (it !in variables) degs[it] else degs[it] + 1u },
+                        cleanedVariables.fold(c) { acc, variable -> acc / optimizedMultiply(one, degs[variable]) }
+                    )
+                }
+        }
+    )
+}
+
+/**
+ * Returns algebraic derivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Field<C>> NumberedPolynomial<C>.nthAntiderivativeWithRespectTo(
+    algebra: A,
+    variable: Int,
+    order: UInt
+): NumberedPolynomial<C> = algebra {
+    if (order == 0u) return this@nthAntiderivativeWithRespectTo
+    NumberedPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    put(
+                        List(max(variable + 1, degs.size)) { if (it != variable) degs[it] else degs[it] + order },
+                        degs[variable].let { deg ->
+                            (deg downTo deg - order + 1u)
+                                .fold(c) { acc, ord -> acc / optimizedMultiply(one, ord) }
+                        }
+                    )
+                }
+        }
+    )
+}
+
+/**
+ * Returns algebraic derivative of received polynomial.
+ */
+@UnstableKMathAPI
+public fun <C, A : Field<C>> NumberedPolynomial<C>.nthAntiderivativeWithRespectTo(
+    algebra: A,
+    variablesAndOrders: Map<Int, UInt>,
+): NumberedPolynomial<C> = algebra {
+    val filteredVariablesAndOrders = variablesAndOrders.filterValues { it != 0u }
+    if (filteredVariablesAndOrders.isEmpty()) return this@nthAntiderivativeWithRespectTo
+    val maxRespectedVariable = filteredVariablesAndOrders.keys.maxOrNull()!!
+    NumberedPolynomial<C>(
+        buildMap(coefficients.size) {
+            coefficients
+                .forEach { (degs, c) ->
+                    put(
+                        List(max(maxRespectedVariable + 1, degs.size)) { degs[it] + filteredVariablesAndOrders.getOrElse(it) { 0u } },
+                        filteredVariablesAndOrders.entries.fold(c) { acc1, (index, order) ->
+                            degs[index].let { deg ->
+                                (deg downTo deg - order + 1u)
+                                    .fold(acc1) { acc2, ord -> acc2 / optimizedMultiply(one, ord) }
+                            }
+                        }
+                    )
+                }
         }
     )
 }
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 a41a6cbd6..9c8cc0090 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
@@ -138,7 +138,7 @@ public fun <C, A> Polynomial<C>.nthDerivative(
     algebra: A,
     order: UInt,
 ): Polynomial<C> where  A : Ring<C>, A : NumericAlgebra<C> = algebra {
-    Polynomial(coefficients.drop(order.toInt()).mapIndexed { index, c -> (1..order.toInt()).fold(c) { acc, i -> acc * number(index + i) } })
+    Polynomial(coefficients.drop(order.toInt()).mapIndexed { index, c -> (index + 1..index + order.toInt()).fold(c) { acc, i -> acc * number(i) } })
 }
 
 /**