kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/AbstractPolynomial.kt

@@ -0,0 +1,354 @@
+/*
+ * 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.operations.*
+import kotlin.js.JsName
+import kotlin.jvm.JvmName
+
+
+/**
+ * Abstraction of polynomials.
+ */
+public interface AbstractPolynomial<C>
+
+/**
+ * Abstraction of ring of polynomials of type [P] over ring of constants of type [C].
+ */
+@Suppress("INAPPLICABLE_JVM_NAME", "PARAMETER_NAME_CHANGED_ON_OVERRIDE")
+public interface AbstractPolynomialSpace<C, P: AbstractPolynomial<C>> : Ring<P> {
+    // region Constant-integer relation
+    /**
+     * Returns sum of the constant and the integer represented as constant (member of underlying ring).
+     *
+     * The operation is equivalent to adding [other] copies of unit of underlying ring to [this].
+     */
+    @JvmName("constantIntPlus")
+    public operator fun C.plus(other: Int): C
+    /**
+     * Returns difference between the constant and the integer represented as constant (member of underlying ring).
+     *
+     * The operation is equivalent to subtraction [other] copies of unit of underlying ring from [this].
+     */
+    @JvmName("constantIntMinus")
+    public operator fun C.minus(other: Int): C
+    /**
+     * Returns product of the constant and the integer represented as constant (member of underlying ring).
+     *
+     * The operation is equivalent to sum of [other] copies of [this].
+     */
+    @JvmName("constantIntTimes")
+    public operator fun C.times(other: Int): C
+    // endregion
+
+    // region Integer-constant relation
+    /**
+     * Returns sum of the integer represented as constant (member of underlying ring) and the constant.
+     *
+     * The operation is equivalent to adding [this] copies of unit of underlying ring to [other].
+     */
+    @JvmName("intConstantPlus")
+    public operator fun Int.plus(other: C): C
+    /**
+     * Returns difference between the integer represented as constant (member of underlying ring) and the constant.
+     *
+     * The operation is equivalent to subtraction [this] copies of unit of underlying ring from [other].
+     */
+    @JvmName("intConstantMinus")
+    public operator fun Int.minus(other: C): C
+    /**
+     * Returns product of the integer represented as constant (member of underlying ring) and the constant.
+     *
+     * The operation is equivalent to sum of [this] copies of [other].
+     */
+    @JvmName("intConstantTimes")
+    public operator fun Int.times(other: C): C
+    // endregion
+
+    // region Polynomial-integer relation
+    /**
+     * Returns sum of the constant and the integer represented as polynomial.
+     *
+     * The operation is equivalent to adding [other] copies of unit polynomial to [this].
+     */
+    public operator fun P.plus(other: Int): P = optimizedAddMultiplied(this, one, other)
+    /**
+     * Returns difference between the constant and the integer represented as polynomial.
+     *
+     * The operation is equivalent to subtraction [other] copies of unit polynomial from [this].
+     */
+    public operator fun P.minus(other: Int): P = optimizedAddMultiplied(this, one, -other)
+    /**
+     * Returns product of the constant and the integer represented as polynomial.
+     *
+     * The operation is equivalent to sum of [other] copies of [this].
+     */
+    public operator fun P.times(other: Int): P = optimizedMultiply(this, other)
+    // endregion
+
+    // region Integer-polynomial relation
+    /**
+     * Returns sum of the integer represented as polynomial and the constant.
+     *
+     * The operation is equivalent to adding [this] copies of unit polynomial to [other].
+     */
+    public operator fun Int.plus(other: P): P = optimizedAddMultiplied(other, one, this)
+    /**
+     * Returns difference between the integer represented as polynomial and the constant.
+     *
+     * The operation is equivalent to subtraction [this] copies of unit polynomial from [other].
+     */
+    public operator fun Int.minus(other: P): P = optimizedAddMultiplied(-other, one, this)
+    /**
+     * Returns product of the integer represented as polynomial and the constant.
+     *
+     * The operation is equivalent to sum of [this] copies of [other].
+     */
+    public operator fun Int.times(other: P): P = optimizedMultiply(other, this)
+    // endregion
+
+    // region Constant-constant relation
+    /**
+     * Returns the same constant.
+     */
+    @JvmName("constantUnaryPlus")
+    @JsName("constantUnaryPlus")
+    public operator fun C.unaryPlus(): C = this
+    /**
+     * Returns negation of the constant.
+     */
+    @JvmName("constantUnaryMinus")
+    @JsName("constantUnaryMinus")
+    public operator fun C.unaryMinus(): C
+    /**
+     * Returns sum of the constants.
+     */
+    @JvmName("constantPlus")
+    @JsName("constantPlus")
+    public operator fun C.plus(other: C): C
+    /**
+     * Returns difference of the constants.
+     */
+    @JvmName("constantMinus")
+    @JsName("constantMinus")
+    public operator fun C.minus(other: C): C
+    /**
+     * Returns product of the constants.
+     */
+    @JvmName("constantTimes")
+    @JsName("constantTimes")
+    public operator fun C.times(other: C): C
+
+    /**
+     * Check if the instant is zero constant.
+     */
+    @JvmName("constantIsZero")
+    public fun C.isZero(): Boolean
+    /**
+     * Check if the instant is NOT zero constant.
+     */
+    @JvmName("constantIsNotZero")
+    public fun C.isNotZero(): Boolean
+    /**
+     * Check if the instant is unit constant.
+     */
+    @JvmName("constantIsOne")
+    public fun C.isOne(): Boolean
+    /**
+     * Check if the instant is NOT unit constant.
+     */
+    @JvmName("constantIsNotOne")
+    public fun C.isNotOne(): Boolean
+    /**
+     * Check if the instant is minus unit constant.
+     */
+    @JvmName("constantIsMinusOne")
+    public fun C.isMinusOne(): Boolean
+    /**
+     * Check if the instant is NOT minus unit constant.
+     */
+    @JvmName("constantIsNotMinusOne")
+    public fun C.isNotMinusOne(): Boolean
+    // endregion
+
+    // region Constant-polynomial relation
+    /**
+     * Returns sum of the constant represented as polynomial and the polynomial.
+     */
+    public operator fun C.plus(other: P): P
+    /**
+     * Returns difference between the constant represented as polynomial and the polynomial.
+     */
+    public operator fun C.minus(other: P): P
+    /**
+     * Returns product of the constant represented as polynomial and the polynomial.
+     */
+    public operator fun C.times(other: P): P
+    // endregion
+
+    // region Polynomial-constant relation
+    /**
+     * Returns sum of the constant represented as polynomial and the polynomial.
+     */
+    public operator fun P.plus(other: C): P
+    /**
+     * Returns difference between the constant represented as polynomial and the polynomial.
+     */
+    public operator fun P.minus(other: C): P
+    /**
+     * Returns product of the constant represented as polynomial and the polynomial.
+     */
+    public operator fun P.times(other: C): P
+    // endregion
+
+    // region Polynomial-polynomial relation
+    /**
+     * Returns the same polynomial.
+     */
+    public override operator fun P.unaryPlus(): P = this
+    /**
+     * Returns negation of the polynomial.
+     */
+    public override operator fun P.unaryMinus(): P
+    /**
+     * Returns sum of the polynomials.
+     */
+    public override operator fun P.plus(other: P): P
+    /**
+     * Returns difference of the polynomials.
+     */
+    public override operator fun P.minus(other: P): P
+    /**
+     * Returns product of the polynomials.
+     */
+    public override operator fun P.times(other: P): P
+
+    /**
+     * Check if the instant is zero polynomial.
+     */
+    public fun P.isZero(): Boolean = this == zero
+    /**
+     * Check if the instant is NOT zero polynomial.
+     */
+    public fun P.isNotZero(): Boolean = this != zero
+    /**
+     * Check if the instant is unit polynomial.
+     */
+    public fun P.isOne(): Boolean = this == one
+    /**
+     * Check if the instant is NOT unit polynomial.
+     */
+    public fun P.isNotOne(): Boolean = this != one
+    /**
+     * Check if the instant is minus unit polynomial.
+     */
+    public fun P.isMinusOne(): Boolean = this == -one
+    /**
+     * Check if the instant is NOT minus unit polynomial.
+     */
+    public fun P.isNotMinusOne(): Boolean = this != -one
+
+    /**
+     * Instance of zero polynomial (zero of the polynomial ring).
+     */
+    override val zero: P
+    /**
+     * Instance of unit polynomial (unit of the polynomial ring).
+     */
+    override val one: P
+
+    /**
+     * Checks equality of the polynomials.
+     */
+    @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "CovariantEquals")
+    public fun P.equals(other: P): Boolean
+    // endregion
+
+    // Not sure is it necessary...
+    // region Polynomial properties
+    /**
+     * Degree of the polynomial, [see also](https://en.wikipedia.org/wiki/Degree_of_a_polynomial). If the polynomial is
+     * zero, degree is -1.
+     */
+    public val P.degree: Int
+
+    /**
+     * Checks if the instant is constant polynomial (of degree no more than 0) over considered ring.
+     */
+    public fun P.isConstant(): Boolean = degree <= 0
+    /**
+     * Checks if the instant is **not** constant polynomial (of degree no more than 0) over considered ring.
+     */
+    public fun P.isNotConstant(): Boolean = !isConstant()
+    /**
+     * Checks if the instant is constant non-zero polynomial (of degree no more than 0) over considered ring.
+     */
+    public fun P.isNonZeroConstant(): Boolean = degree == 0
+    /**
+     * Checks if the instant is **not** constant non-zero polynomial (of degree no more than 0) over considered ring.
+     */
+    public fun P.isNotNonZeroConstant(): Boolean = !isNonZeroConstant()
+    /**
+     * If polynomial is a constant polynomial represents and returns it as constant.
+     * Otherwise, (when the polynomial is not constant polynomial) returns `null`.
+     */
+    public fun P.asConstantOrNull(): C?
+    /**
+     * If polynomial is a constant polynomial represents and returns it as constant.
+     * Otherwise, (when the polynomial is not constant polynomial) raises corresponding exception.
+     */
+    public fun P.asConstant(): C = asConstantOrNull() ?: error("Can not represent non-constant polynomial as a constant")
+    // endregion
+
+    // region Legacy of Ring interface
+    override fun add(left: P, right: P): P = left + right
+    override fun multiply(left: P, right: P): P = left * right
+    // endregion
+
+    public companion object {
+        // TODO: All of this should be moved to algebraic structures' place for utilities
+        // TODO: Move receiver to context receiver
+        /**
+         * Multiplication of element and integer.
+         *
+         * @receiver the multiplicand.
+         * @param other the multiplier.
+         * @return the difference.
+         */
+        internal tailrec fun <C> Group<C>.optimizedMultiply(arg: C, other: Int): C =
+            when {
+                other == 0 -> zero
+                other == 1 -> arg
+                other == -1 -> -arg
+                other % 2 == 0 -> optimizedMultiply(arg + arg, other / 2)
+                other % 2 == 1 -> optimizedAddMultiplied(arg, arg + arg, other / 2)
+                other % 2 == -1 -> optimizedAddMultiplied(-arg, arg + arg, other / 2)
+                else -> error("Error in multiplication group instant by integer: got reminder by division by 2 different from 0, 1 and -1")
+            }
+
+        // TODO: Move receiver to context receiver
+        /**
+         * Adds product of [arg] and [multiplier] to [base].
+         *
+         * @receiver the algebra to provide multiplication.
+         * @param base the augend.
+         * @param arg the multiplicand.
+         * @param multiplier the multiplier.
+         * @return sum of the augend [base] and product of the multiplicand [arg] and the multiplier [multiplier].
+         * @author Gleb Minaev
+         */
+        internal tailrec fun <C> Group<C>.optimizedAddMultiplied(base: C, arg: C, multiplier: Int): C =
+            when {
+                multiplier == 0 -> base
+                multiplier == 1 -> base + arg
+                multiplier == -1 -> base - arg
+                multiplier % 2 == 0 -> optimizedAddMultiplied(base, arg + arg, multiplier / 2)
+                multiplier % 2 == 1 -> optimizedAddMultiplied(base + arg, arg + arg, multiplier / 2)
+                multiplier % 2 == -1 -> optimizedAddMultiplied(base + arg, arg + arg, multiplier / 2)
+                else -> error("Error in multiplication group instant by integer: got reminder by division by 2 different from 0, 1 and -1")
+            }
+    }
+}

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt

@@ -117,16 +117,16 @@ public fun <T : Comparable<T>> PiecewisePolynomial(
  * Return a value of polynomial function with given [ring] a given [arg] or null if argument is outside piecewise
  * definition.
  */
-public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.value(ring: C, arg: T): T? =
+public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.substitute(ring: C, arg: T): T? =
-    findPiece(arg)?.value(ring, arg)
+    findPiece(arg)?.substitute(ring, arg)
 
 /**
  * Convert this polynomial to a function returning nullable value (null if argument is outside piecewise range).
  */
-public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { value(ring, it) }
+public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { substitute(ring, it) }
 
 /**
  * Convert this polynomial to a function using [defaultValue] for arguments outside the piecewise range.
  */
 public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C, defaultValue: T): (T) -> T =
-    { value(ring, it) ?: defaultValue }
+    { substitute(ring, it) ?: defaultValue }

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt

@@ -5,128 +5,371 @@
 
 package space.kscience.kmath.functions
 
-import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.functions.AbstractPolynomialSpace.Companion.optimizedAddMultiplied
+import space.kscience.kmath.functions.AbstractPolynomialSpace.Companion.optimizedMultiply
 import space.kscience.kmath.operations.*
-import kotlin.contracts.InvocationKind
+import kotlin.jvm.JvmName
-import kotlin.contracts.contract
 import kotlin.math.max
-import kotlin.math.pow
+import kotlin.math.min
 
 /**
  * Polynomial coefficients model without fixation on specific context they are applied to.
  *
  * @param coefficients constant is the leftmost coefficient.
  */
-public class Polynomial<out T>(public val coefficients: List<T>) {
+public class Polynomial<T>(public val coefficients: List<T>) : AbstractPolynomial<T> {
     override fun toString(): String = "Polynomial$coefficients"
+
+//    public companion object {
+//        /**
+//         * Default name of variables used in string representations.
+//         *
+//         * @see Polynomial.toString
+//         */
+//        public var defaultVariableName: String = "x"
+//
+//        /**
+//         * Represents result of division with remainder.
+//         */
+//        public data class DividingResult<C>(
+//            val quotient: Polynomial<C>,
+//            val reminder: Polynomial<C>
+//        )
+//    }
 }
 
 /**
- * Returns a [Polynomial] instance with given [coefficients].
+ * Represents internal [Polynomial] errors.
+ */
+internal class PolynomialError(message: String): Error(message)
+
+// region Constructors and converters
+
+/**
+ * Returns a [Polynomial] instance with given [coefficients]. The collection of coefficients will be reversed if
+ * [reverse] parameter is true.
  */
 @Suppress("FunctionName")
-public fun <T> Polynomial(vararg coefficients: T): Polynomial<T> = Polynomial(coefficients.toList())
+public fun <T> Polynomial(coefficients: List<T>, reverse: Boolean = false): Polynomial<T> =
+    Polynomial(with(coefficients) { if (reverse) reversed() else this })
 
 /**
- * Evaluates the value of the given double polynomial for given double argument.
+ * Returns a [Polynomial] instance with given [coefficients]. The collection of coefficients will be reversed if
+ * [reverse] parameter is true.
  */
-public fun Polynomial<Double>.value(arg: Double): Double = coefficients.reduceIndexed { index, acc, c ->
+@Suppress("FunctionName")
-    acc + c * arg.pow(index)
+public fun <T> Polynomial(vararg coefficients: T, reverse: Boolean = false): Polynomial<T> =
-}
+    Polynomial(with(coefficients) { if (reverse) reversed() else toList() })
+
+public fun <T> T.asPolynomial() : Polynomial<T> = Polynomial(listOf(this))
+
+// endregion
 
 /**
- * Evaluates the value of the given polynomial for given argument.
+ * Space of polynomials constructed over ring.
- * https://en.wikipedia.org/wiki/Horner%27s_method
- */
-public fun <T, C : Ring<T>> Polynomial<T>.value(ring: C, arg: T): T = ring {
-    if (coefficients.isEmpty()) return@ring zero
-    var result: T = coefficients.last()
-    for (j in coefficients.size - 2 downTo 0) {
-        result = (arg * result) + coefficients[j]
-    }
-    return result
-}
-
-/**
- * Represent the polynomial as a regular context-less function.
- */
-public fun <T, C : Ring<T>> Polynomial<T>.asFunction(ring: C): (T) -> T = { value(ring, it) }
-
-/**
- * Create a polynomial witch represents differentiated version of this polynomial
- */
-@UnstableKMathAPI
-public fun <T, A> Polynomial<T>.differentiate(
-    algebra: A,
-): Polynomial<T> where  A : Ring<T>, A : NumericAlgebra<T> = algebra {
-    Polynomial(coefficients.drop(1).mapIndexed { index, t -> number(index) * t })
-}
-
-/**
- * Create a polynomial witch represents indefinite integral version of this polynomial
- */
-@UnstableKMathAPI
-public fun <T, A> Polynomial<T>.integrate(
-    algebra: A,
-): Polynomial<T> where  A : Field<T>, A : NumericAlgebra<T> = algebra {
-    val integratedCoefficients = buildList(coefficients.size + 1) {
-        add(zero)
-        coefficients.forEachIndexed{ index, t -> add(t / (number(index) + one)) }
-    }
-    Polynomial(integratedCoefficients)
-}
-
-/**
- * Compute a definite integral of a given polynomial in a [range]
- */
-@UnstableKMathAPI
-public fun <T : Comparable<T>> Polynomial<T>.integrate(
-    algebra: Field<T>,
-    range: ClosedRange<T>,
-): T = algebra {
-    val integral = integrate(algebra)
-    integral.value(algebra, range.endInclusive) - integral.value(algebra, range.start)
-}
-
-/**
- * Space of polynomials.
  *
- * @param T the type of operated polynomials.
+ * @param C the type of constants. Polynomials have them as a coefficients in their terms.
- * @param C the intersection of [Ring] of [T] and [ScaleOperations] of [T].
+ * @param A type of underlying ring of constants. It's [Ring] of [C].
- * @param ring the [C] instance.
+ * @param ring underlying ring of constants of type [A].
  */
-public class PolynomialSpace<T, C>(
+@Suppress("INAPPLICABLE_JVM_NAME") // KT-31420
-    private val ring: C,
+public open class PolynomialSpace<C, A : Ring<C>>(
-) : Group<Polynomial<T>>, ScaleOperations<Polynomial<T>> where C : Ring<T>, C : ScaleOperations<T> {
+    public val ring: A,
-    override val zero: Polynomial<T> = Polynomial(emptyList())
+) : AbstractPolynomialSpace<C, Polynomial<C>>{
+    // region Constant-integer relation
+    @JvmName("constantIntPlus")
+    public override operator fun C.plus(other: Int): C = ring { optimizedAddMultiplied(this@plus, one, other) }
+    @JvmName("constantIntMinus")
+    public override operator fun C.minus(other: Int): C = ring { optimizedAddMultiplied(this@minus, one, -other) }
+    @JvmName("constantIntTimes")
+    public override operator fun C.times(other: Int): C = ring { optimizedMultiply(this@times, other) }
+    // endregion
 
-    override fun Polynomial<T>.unaryMinus(): Polynomial<T> = ring {
+    // region Integer-constant relation
+    @JvmName("intConstantPlus")
+    public override operator fun Int.plus(other: C): C = ring { optimizedAddMultiplied(other, one, this@plus) }
+    @JvmName("intConstantMinus")
+    public override operator fun Int.minus(other: C): C = ring { optimizedAddMultiplied(-other, one, this@minus) }
+    @JvmName("intConstantTimes")
+    public override operator fun Int.times(other: C): C = ring { optimizedMultiply(other, this@times) }
+    // endregion
+
+    // region Polynomial-integer relation
+    public override operator fun Polynomial<C>.plus(other: Int): Polynomial<C> =
+        if (other == 0) this
+        else
+            Polynomial(
+                coefficients
+                    .toMutableList()
+                    .apply {
+                        if (isEmpty()) this[0] = ring.zero + other
+                        else this[0] = this[0]!! + other
+                    }
+            )
+    public override operator fun Polynomial<C>.minus(other: Int): Polynomial<C> =
+        if (other == 0) this
+        else
+            Polynomial(
+                coefficients
+                    .toMutableList()
+                    .apply {
+                        if (isEmpty()) this[0] = ring.zero - other
+                        else this[0] = this[0]!! - other
+                    }
+            )
+    public override operator fun Polynomial<C>.times(other: Int): Polynomial<C> =
+        if (other == 0) zero
+        else Polynomial(
+            coefficients
+                .subList(0, degree + 1)
+                .map { it * other }
+        )
+    // endregion
+
+    // region Integer-polynomial relation
+    public override operator fun Int.plus(other: Polynomial<C>): Polynomial<C> =
+        if (this == 0) other
+        else
+            Polynomial(
+                other.coefficients
+                    .toMutableList()
+                    .apply {
+                        if (isEmpty()) this[0] = ring.zero + this@plus
+                        else this[0] = this[0]!! + this@plus
+                    }
+            )
+    public override operator fun Int.minus(other: Polynomial<C>): Polynomial<C> =
+        if (this == 0) other
+        else
+            Polynomial(
+                other.coefficients
+                    .toMutableList()
+                    .apply {
+                        if (isEmpty()) this[0] = ring.zero - this@minus
+                        else this[0] = this[0]!! - this@minus
+                    }
+            )
+    public override operator fun Int.times(other: Polynomial<C>): Polynomial<C> =
+        if (this == 0) zero
+        else Polynomial(
+            other.coefficients
+                .subList(0, other.degree + 1)
+                .map { it * this }
+        )
+    // endregion
+
+    // region Constant-constant relation
+    @JvmName("constantUnaryMinus")
+    public override operator fun C.unaryMinus(): C = ring { -this@unaryMinus }
+    @JvmName("constantPlus")
+    public override operator fun C.plus(other: C): C = ring { this@plus + other }
+    @JvmName("constantMinus")
+    public override operator fun C.minus(other: C): C = ring { this@minus - other }
+    @JvmName("constantTimes")
+    public override operator fun C.times(other: C): C = ring { this@times * other }
+    @JvmName("constantIsZero")
+    public override fun C.isZero(): Boolean = ring { this == zero }
+    @JvmName("constantIsNotZero")
+    public override fun C.isNotZero(): Boolean = ring { this != zero }
+    @JvmName("constantIsOne")
+    public override fun C.isOne(): Boolean = ring { this == one }
+    @JvmName("constantIsNotOne")
+    public override fun C.isNotOne(): Boolean = ring { this != one }
+    @JvmName("constantIsMinusOne")
+    public override fun C.isMinusOne(): Boolean = ring { this == -one }
+    @JvmName("constantIsNotMinusOne")
+    public override fun C.isNotMinusOne(): Boolean = ring { this != -one }
+    // endregion
+
+    // region Constant-polynomial relation
+    public override operator fun C.plus(other: Polynomial<C>): Polynomial<C> =
+        with(other.coefficients) {
+            if (isEmpty()) Polynomial(listOf(this@plus))
+            else Polynomial(
+                toMutableList()
+                    .apply {
+                        this[0] += this@plus
+                    }
+            )
+        }
+//        if (degree == -1) UnivariatePolynomial(other) else UnivariatePolynomial(
+//            listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1)
+//        )
+    public override operator fun C.minus(other: Polynomial<C>): Polynomial<C> =
+        with(other.coefficients) {
+            if (isEmpty()) Polynomial(listOf(-this@minus))
+            else Polynomial(
+                toMutableList()
+                    .apply {
+                        this[0] -= this@minus
+                    }
+            )
+        }
+//        if (degree == -1) UnivariatePolynomial(other) else UnivariatePolynomial(
+//            listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1)
+//        )
+    public override operator fun C.times(other: Polynomial<C>): Polynomial<C> =
+        Polynomial(
+            other.coefficients
+//                .subList(0, other.degree + 1)
+                .map { it * this }
+        )
+    // endregion
+
+    // region Polynomial-constant relation
+    public override operator fun Polynomial<C>.plus(other: C): Polynomial<C> =
+        if (other.isZero()) this
+        else with(coefficients) {
+            if (isEmpty()) Polynomial(listOf(other))
+            else Polynomial(
+                toMutableList()
+                    .apply {
+                        this[0] += other
+                    }
+            )
+        }
+//        if (degree == -1) UnivariatePolynomial(other) else UnivariatePolynomial(
+//            listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1)
+//        )
+    public override operator fun Polynomial<C>.minus(other: C): Polynomial<C> =
+        with(coefficients) {
+            if (isEmpty()) Polynomial(listOf(other))
+            else Polynomial(
+                toMutableList()
+                    .apply {
+                        this[0] -= other
+                    }
+            )
+        }
+//        if (degree == -1) UnivariatePolynomial(-other) else UnivariatePolynomial(
+//            listOf(coefficients[0] - other) + coefficients.subList(1, degree + 1)
+//        )
+    public override operator fun Polynomial<C>.times(other: C): Polynomial<C> =
+        Polynomial(
+            coefficients
+//                .subList(0, degree + 1)
+                .map { it * other }
+        )
+    // endregion
+
+    // region Polynomial-polynomial relation
+    public override operator fun Polynomial<C>.unaryMinus(): Polynomial<C> =
         Polynomial(coefficients.map { -it })
+    public override operator fun Polynomial<C>.plus(other: Polynomial<C>): Polynomial<C> =
+        Polynomial(
+            (0..max(degree, other.degree))
+                .map {
+                    when {
+                        it > degree -> other.coefficients[it]
+                        it > other.degree -> coefficients[it]
+                        else -> coefficients[it] + other.coefficients[it]
                     }
-
+                }
-    override fun add(left: Polynomial<T>, right: Polynomial<T>): Polynomial<T> {
+                .ifEmpty { listOf(ring.zero) }
-        val dim = max(left.coefficients.size, right.coefficients.size)
+        )
-
+    public override operator fun Polynomial<C>.minus(other: Polynomial<C>): Polynomial<C> =
-        return ring {
+        Polynomial(
-            Polynomial(List(dim) { index ->
+            (0..max(degree, other.degree))
-                left.coefficients.getOrElse(index) { zero } + right.coefficients.getOrElse(index) { zero }
+                .map {
-            })
+                    when {
+                        it > degree -> -other.coefficients[it]
+                        it > other.degree -> coefficients[it]
+                        else -> coefficients[it] - other.coefficients[it]
+                    }
+                }
+                .ifEmpty { listOf(ring.zero) }
+        )
+    public override operator fun Polynomial<C>.times(other: Polynomial<C>): Polynomial<C> {
+        val thisDegree = degree
+        val otherDegree = other.degree
+        return when {
+            thisDegree == -1 -> this
+            otherDegree == -1 -> other
+            else ->
+                Polynomial(
+                    (0..(thisDegree + otherDegree))
+                        .map { d ->
+                            (max(0, d - otherDegree)..(min(thisDegree, d)))
+                                .map { coefficients[it] * other.coefficients[d - it] }
+                                .reduce { acc, rational -> acc + rational }
+                        }
+                )
         }
     }
 
-    override fun scale(a: Polynomial<T>, value: Double): Polynomial<T> =
+    public override fun Polynomial<C>.isZero(): Boolean = coefficients.all { it.isZero() }
-        ring { Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * value }) }
+    public override fun Polynomial<C>.isNotZero(): Boolean = coefficients.any { it.isNotZero() }
+    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?
+    public override fun Polynomial<C>.isNotOne(): Boolean = !isOne()
+    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?
+    public override fun Polynomial<C>.isNotMinusOne(): Boolean = !isMinusOne()
+
+    override val zero: Polynomial<C> = Polynomial(emptyList())
+    override val one: Polynomial<C> = Polynomial(listOf(ring.one))
+
+    @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "CovariantEquals")
+    public override fun Polynomial<C>.equals(other: Polynomial<C>): Boolean =
+        when {
+            this === other -> true
+            else -> {
+                if (this.degree == other.degree)
+                    (0..degree).all { coefficients[it] == other.coefficients[it] }
+                else false
+            }
+        }
+    // endregion
+
+    // Not sure is it necessary...
+    // region Polynomial properties
+    public override val Polynomial<C>.degree: Int get() = coefficients.indexOfLast { it != ring.zero }
+
+    public override fun Polynomial<C>.asConstantOrNull(): C? =
+        with(coefficients) {
+            when {
+                isEmpty() -> ring.zero
+                degree > 0 -> null
+                else -> first()
+            }
+        }
+    public override fun Polynomial<C>.asConstant(): C = asConstantOrNull() ?: error("Can not represent non-constant polynomial as a constant")
+
+    @Suppress("NOTHING_TO_INLINE")
+    public inline fun Polynomial<C>.substitute(argument: C): C = this.substitute(ring, argument)
+    @Suppress("NOTHING_TO_INLINE")
+    public inline fun Polynomial<C>.substitute(argument: Polynomial<C>): Polynomial<C> = this.substitute(ring, argument)
+
+    @Suppress("NOTHING_TO_INLINE")
+    public inline fun Polynomial<C>.asFunction(): (C) -> C = { this.substitute(ring, it) }
+    @Suppress("NOTHING_TO_INLINE")
+    public inline fun Polynomial<C>.asFunctionOnConstants(): (C) -> C = { this.substitute(ring, it) }
+    @Suppress("NOTHING_TO_INLINE")
+    public inline fun Polynomial<C>.asFunctionOnPolynomials(): (Polynomial<C>) -> Polynomial<C> = { this.substitute(ring, it) }
 
     /**
      * Evaluates the polynomial for the given value [arg].
      */
-    public operator fun Polynomial<T>.invoke(arg: T): T = value(ring, arg)
+    @Suppress("NOTHING_TO_INLINE")
-
+    public inline operator fun Polynomial<C>.invoke(argument: C): C = this.substitute(ring, argument)
-    public fun Polynomial<T>.asFunction(): (T) -> T = asFunction(ring)
+    @Suppress("NOTHING_TO_INLINE")
+    public inline operator fun Polynomial<C>.invoke(argument: Polynomial<C>): Polynomial<C> = this.substitute(ring, argument)
+    // endregion
 
 }
 
-public inline fun <T, C, R> C.polynomial(block: PolynomialSpace<T, C>.() -> R): R where C : Ring<T>, C : ScaleOperations<T> {
+/**
-    contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
+ * Space of polynomials constructed over ring.
-    return PolynomialSpace(this).block()
+ *
+ * @param C the type of constants. Polynomials have them as a coefficients in their terms.
+ * @param A type of underlying ring of constants. It's [Ring] of [C].
+ * @param ring underlying ring of constants of type [A].
+ */
+public class ScalablePolynomialSpace<C, A>(
+    ring: A,
+) : PolynomialSpace<C, A>(ring), ScaleOperations<Polynomial<C>> where A : Ring<C>, A : ScaleOperations<C> {
+
+    override fun scale(a: Polynomial<C>, value: Double): Polynomial<C> =
+        ring { Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * value }) }
+
 }

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt

@@ -0,0 +1,139 @@
+/*
+ * 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.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.*
+import kotlin.contracts.InvocationKind
+import kotlin.contracts.contract
+import kotlin.jvm.JvmName
+import kotlin.math.pow
+
+
+// region Utilities
+
+/**
+ * Removes zeros on the end of the coefficient list of polynomial.
+ */
+//context(PolynomialSpace<C, A>)
+//fun <C, A: Ring<C>> Polynomial<C>.removeZeros() : Polynomial<C> =
+//    if (degree > -1) Polynomial(coefficients.subList(0, degree + 1)) else zero
+
+/**
+ * Crates a [PolynomialSpace] over received ring.
+ */
+public fun <C, A : Ring<C>> A.polynomial(): PolynomialSpace<C, A> =
+    PolynomialSpace(this)
+
+/**
+ * Crates a [PolynomialSpace]'s scope over received ring.
+ */
+public inline fun <C, A : Ring<C>, R> A.polynomial(block: PolynomialSpace<C, A>.() -> R): R {
+    contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
+    return PolynomialSpace(this).block()
+}
+
+/**
+ * Crates a [ScalablePolynomialSpace] over received scalable ring.
+ */
+public fun <C, A> A.scalablePolynomial(): ScalablePolynomialSpace<C, A> where A : Ring<C>, A : ScaleOperations<C> =
+    ScalablePolynomialSpace(this)
+
+/**
+ * Crates a [ScalablePolynomialSpace]'s scope over received scalable ring.
+ */
+public inline fun <C, A, R> A.scalablePolynomial(block: ScalablePolynomialSpace<C, A>.() -> R): R where A : Ring<C>, A : ScaleOperations<C> {
+    contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
+    return ScalablePolynomialSpace(this).block()
+}
+
+// endregion
+
+// region Polynomial substitution and functional representation
+
+// TODO: May be apply Horner's method too?
+/**
+ * Evaluates the value of the given double polynomial for given double argument.
+ */
+public fun Polynomial<Double>.substitute(arg: Double): Double =
+    coefficients.reduceIndexedOrNull { index, acc, c ->
+        acc + c * arg.pow(index)
+    } ?: .0
+
+/**
+ * Evaluates the value of the given polynomial for given argument.
+ *
+ * It is an implementation of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
+ */
+public fun <C> Polynomial<C>.substitute(ring: Ring<C>, arg: C): C = ring {
+    if (coefficients.isEmpty()) return@ring zero
+    var result: C = coefficients.last()
+    for (j in coefficients.size - 2 downTo 0) {
+        result = (arg * result) + coefficients[j]
+    }
+    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]
+    }
+    return result
+}
+
+/**
+ * Represent the polynomial as a regular context-less function.
+ */
+public fun <C, A : Ring<C>> Polynomial<C>.asFunction(ring: A): (C) -> C = { substitute(ring, it) }
+
+/**
+ * Represent the polynomial as a regular context-less function.
+ */
+public fun <C, A : Ring<C>> Polynomial<C>.asPolynomialFunctionOver(ring: A): (Polynomial<C>) -> Polynomial<C> = { substitute(ring, it) }
+
+// endregion
+
+// region Algebraic derivative and antiderivative
+
+/**
+ * Returns algebraic derivative of received polynomial.
+ */
+@UnstableKMathAPI
+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, t -> number(index) * t })
+}
+
+/**
+ * Create a polynomial witch represents indefinite integral version of this polynomial
+ */
+@UnstableKMathAPI
+public fun <C, A> Polynomial<C>.antiderivative(
+    algebra: A,
+): 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) + one)) }
+    }
+    Polynomial(integratedCoefficients)
+}
+
+/**
+ * Compute a definite integral of a given polynomial in a [range]
+ */
+@UnstableKMathAPI
+public fun <C : Comparable<C>> Polynomial<C>.integrate(
+    algebra: Field<C>,
+    range: ClosedRange<C>,
+): C = algebra {
+    val integral = antiderivative(algebra)
+    integral.substitute(algebra, range.endInclusive) - integral.substitute(algebra, range.start)
+}

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt

@@ -7,6 +7,7 @@ package space.kscience.kmath.integration
 
 import space.kscience.kmath.functions.PiecewisePolynomial
 import space.kscience.kmath.functions.integrate
+import space.kscience.kmath.functions.antiderivative
 import space.kscience.kmath.interpolation.PolynomialInterpolator
 import space.kscience.kmath.interpolation.SplineInterpolator
 import space.kscience.kmath.interpolation.interpolatePolynomials
@@ -23,7 +24,7 @@ import space.kscience.kmath.structures.MutableBufferFactory
 @OptIn(PerformancePitfall::class)
 @UnstableKMathAPI
 public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(algebra: Field<T>): PiecewisePolynomial<T> =
-    PiecewisePolynomial(pieces.map { it.first to it.second.integrate(algebra) })
+    PiecewisePolynomial(pieces.map { it.first to it.second.antiderivative(algebra) })
 
 /**
  * Compute definite integral of given [PiecewisePolynomial] piece by piece in a given [range]

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt

@@ -10,7 +10,7 @@ package space.kscience.kmath.interpolation
 import space.kscience.kmath.data.XYColumnarData
 import space.kscience.kmath.functions.PiecewisePolynomial
 import space.kscience.kmath.functions.asFunction
-import space.kscience.kmath.functions.value
+import space.kscience.kmath.functions.substitute
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Ring
 import space.kscience.kmath.structures.Buffer
@@ -34,7 +34,7 @@ public interface PolynomialInterpolator<T : Comparable<T>> : Interpolator<T, T,
     public fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T>
 
     override fun interpolate(points: XYColumnarData<T, T, T>): (T) -> T = { x ->
-        interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue()
+        interpolatePolynomials(points).substitute(algebra, x) ?: getDefaultValue()
     }
 }
 

kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/PolynomialTest.kt

@@ -5,13 +5,22 @@
 
 package space.kscience.kmath.functions
 
+import space.kscience.kmath.operations.algebra
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
 class PolynomialTest {
+    @Test
+    fun simple_polynomial_test() {
+        Double.algebra.scalablePolynomial {
+            val x = Polynomial(listOf(0.0, 1.0))
+            val polynomial = x * x - 2 * x + 1
+            assertEquals(0.0, polynomial.substitute(1.0), 0.001)
+        }
+    }
     @Test
     fun testIntegration() {
         val polynomial = Polynomial(1.0, -2.0, 1.0)
-        assertEquals(0.0, polynomial.value(1.0), 0.001)
+        assertEquals(0.0, polynomial.substitute(1.0), 0.001)
     }
 }