diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/AbstractPolynomial.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/AbstractPolynomial.kt new file mode 100644 index 000000000..f6a617656 --- /dev/null +++ b/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 + +/** + * 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> : Ring

{ + // 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 Group.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 Group.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") + } + } +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt index 16af7f555..cfd21d552 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt @@ -117,16 +117,16 @@ public fun > PiecewisePolynomial( * Return a value of polynomial function with given [ring] a given [arg] or null if argument is outside piecewise * definition. */ -public fun , C : Ring> PiecewisePolynomial.value(ring: C, arg: T): T? = - findPiece(arg)?.value(ring, arg) +public fun , C : Ring> PiecewisePolynomial.substitute(ring: C, arg: T): T? = + findPiece(arg)?.substitute(ring, arg) /** * Convert this polynomial to a function returning nullable value (null if argument is outside piecewise range). */ -public fun , C : Ring> PiecewisePolynomial.asFunction(ring: C): (T) -> T? = { value(ring, it) } +public fun , C : Ring> PiecewisePolynomial.asFunction(ring: C): (T) -> T? = { substitute(ring, it) } /** * Convert this polynomial to a function using [defaultValue] for arguments outside the piecewise range. */ public fun , C : Ring> PiecewisePolynomial.asFunction(ring: C, defaultValue: T): (T) -> T = - { value(ring, it) ?: defaultValue } + { substitute(ring, it) ?: defaultValue } 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 a36d36f52..30280a396 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 @@ -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.contracts.contract +import kotlin.jvm.JvmName 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(public val coefficients: List) { +public class Polynomial(public val coefficients: List) : AbstractPolynomial { 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( +// val quotient: Polynomial, +// val reminder: Polynomial +// ) +// } } /** - * 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 Polynomial(vararg coefficients: T): Polynomial = Polynomial(coefficients.toList()) +public fun Polynomial(coefficients: List, reverse: Boolean = false): Polynomial = + 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.value(arg: Double): Double = coefficients.reduceIndexed { index, acc, c -> - acc + c * arg.pow(index) -} +@Suppress("FunctionName") +public fun Polynomial(vararg coefficients: T, reverse: Boolean = false): Polynomial = + Polynomial(with(coefficients) { if (reverse) reversed() else toList() }) + +public fun T.asPolynomial() : Polynomial = Polynomial(listOf(this)) + +// endregion /** - * Evaluates the value of the given polynomial for given argument. - * https://en.wikipedia.org/wiki/Horner%27s_method - */ -public fun > Polynomial.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 > Polynomial.asFunction(ring: C): (T) -> T = { value(ring, it) } - -/** - * Create a polynomial witch represents differentiated version of this polynomial - */ -@UnstableKMathAPI -public fun Polynomial.differentiate( - algebra: A, -): Polynomial where A : Ring, A : NumericAlgebra = 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 Polynomial.integrate( - algebra: A, -): Polynomial where A : Field, A : NumericAlgebra = 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 > Polynomial.integrate( - algebra: Field, - range: ClosedRange, -): T = algebra { - val integral = integrate(algebra) - integral.value(algebra, range.endInclusive) - integral.value(algebra, range.start) -} - -/** - * Space of polynomials. + * Space of polynomials constructed over ring. * - * @param T the type of operated polynomials. - * @param C the intersection of [Ring] of [T] and [ScaleOperations] of [T]. - * @param ring the [C] instance. + * @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 PolynomialSpace( - private val ring: C, -) : Group>, ScaleOperations> where C : Ring, C : ScaleOperations { - override val zero: Polynomial = Polynomial(emptyList()) +@Suppress("INAPPLICABLE_JVM_NAME") // KT-31420 +public open class PolynomialSpace>( + public val ring: A, +) : AbstractPolynomialSpace>{ + // 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.unaryMinus(): Polynomial = 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.plus(other: Int): Polynomial = + 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.minus(other: Int): Polynomial = + 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.times(other: Int): Polynomial = + 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): Polynomial = + 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): Polynomial = + 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): Polynomial = + 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): Polynomial = + 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): Polynomial = + 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): Polynomial = + Polynomial( + other.coefficients +// .subList(0, other.degree + 1) + .map { it * this } + ) + // endregion + + // region Polynomial-constant relation + public override operator fun Polynomial.plus(other: C): Polynomial = + 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.minus(other: C): Polynomial = + 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.times(other: C): Polynomial = + Polynomial( + coefficients +// .subList(0, degree + 1) + .map { it * other } + ) + // endregion + + // region Polynomial-polynomial relation + public override operator fun Polynomial.unaryMinus(): Polynomial = Polynomial(coefficients.map { -it }) - } - - override fun add(left: Polynomial, right: Polynomial): Polynomial { - val dim = max(left.coefficients.size, right.coefficients.size) - - return ring { - Polynomial(List(dim) { index -> - left.coefficients.getOrElse(index) { zero } + right.coefficients.getOrElse(index) { zero } - }) + public override operator fun Polynomial.plus(other: Polynomial): Polynomial = + Polynomial( + (0..max(degree, other.degree)) + .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.minus(other: Polynomial): Polynomial = + Polynomial( + (0..max(degree, other.degree)) + .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.times(other: Polynomial): Polynomial { + 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, value: Double): Polynomial = - ring { Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * value }) } + public override fun Polynomial.isZero(): Boolean = coefficients.all { it.isZero() } + public override fun Polynomial.isNotZero(): Boolean = coefficients.any { it.isNotZero() } + public override fun Polynomial.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.isNotOne(): Boolean = !isOne() + public override fun Polynomial.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.isNotMinusOne(): Boolean = !isMinusOne() + + override val zero: Polynomial = Polynomial(emptyList()) + override val one: Polynomial = Polynomial(listOf(ring.one)) + + @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "CovariantEquals") + public override fun Polynomial.equals(other: Polynomial): 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.degree: Int get() = coefficients.indexOfLast { it != ring.zero } + + public override fun Polynomial.asConstantOrNull(): C? = + with(coefficients) { + when { + isEmpty() -> ring.zero + degree > 0 -> null + else -> first() + } + } + public override fun Polynomial.asConstant(): C = asConstantOrNull() ?: error("Can not represent non-constant polynomial as a constant") + + @Suppress("NOTHING_TO_INLINE") + public inline fun Polynomial.substitute(argument: C): C = this.substitute(ring, argument) + @Suppress("NOTHING_TO_INLINE") + public inline fun Polynomial.substitute(argument: Polynomial): Polynomial = this.substitute(ring, argument) + + @Suppress("NOTHING_TO_INLINE") + public inline fun Polynomial.asFunction(): (C) -> C = { this.substitute(ring, it) } + @Suppress("NOTHING_TO_INLINE") + public inline fun Polynomial.asFunctionOnConstants(): (C) -> C = { this.substitute(ring, it) } + @Suppress("NOTHING_TO_INLINE") + public inline fun Polynomial.asFunctionOnPolynomials(): (Polynomial) -> Polynomial = { this.substitute(ring, it) } /** * Evaluates the polynomial for the given value [arg]. */ - public operator fun Polynomial.invoke(arg: T): T = value(ring, arg) - - public fun Polynomial.asFunction(): (T) -> T = asFunction(ring) + @Suppress("NOTHING_TO_INLINE") + public inline operator fun Polynomial.invoke(argument: C): C = this.substitute(ring, argument) + @Suppress("NOTHING_TO_INLINE") + public inline operator fun Polynomial.invoke(argument: Polynomial): Polynomial = this.substitute(ring, argument) + // endregion } -public inline fun C.polynomial(block: PolynomialSpace.() -> R): R where C : Ring, C : ScaleOperations { - contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) } - return PolynomialSpace(this).block() +/** + * Space of polynomials constructed over ring. + * + * @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( + ring: A, +) : PolynomialSpace(ring), ScaleOperations> where A : Ring, A : ScaleOperations { + + override fun scale(a: Polynomial, value: Double): Polynomial = + ring { Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * value }) } + } 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 new file mode 100644 index 000000000..1a3eb7874 --- /dev/null +++ b/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) +//fun > Polynomial.removeZeros() : Polynomial = +// if (degree > -1) Polynomial(coefficients.subList(0, degree + 1)) else zero + +/** + * Crates a [PolynomialSpace] over received ring. + */ +public fun > A.polynomial(): PolynomialSpace = + PolynomialSpace(this) + +/** + * Crates a [PolynomialSpace]'s scope over received ring. + */ +public inline fun , R> A.polynomial(block: PolynomialSpace.() -> R): R { + contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) } + return PolynomialSpace(this).block() +} + +/** + * Crates a [ScalablePolynomialSpace] over received scalable ring. + */ +public fun A.scalablePolynomial(): ScalablePolynomialSpace where A : Ring, A : ScaleOperations = + ScalablePolynomialSpace(this) + +/** + * Crates a [ScalablePolynomialSpace]'s scope over received scalable ring. + */ +public inline fun A.scalablePolynomial(block: ScalablePolynomialSpace.() -> R): R where A : Ring, A : ScaleOperations { + 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.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 Polynomial.substitute(ring: Ring, 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 Polynomial.substitute(ring: Ring, arg: Polynomial) : Polynomial = ring.polynomial { + if (coefficients.isEmpty()) return zero + var result: Polynomial = 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 > Polynomial.asFunction(ring: A): (C) -> C = { substitute(ring, it) } + +/** + * Represent the polynomial as a regular context-less function. + */ +public fun > Polynomial.asPolynomialFunctionOver(ring: A): (Polynomial) -> Polynomial = { substitute(ring, it) } + +// endregion + +// region Algebraic derivative and antiderivative + +/** + * Returns algebraic derivative of received polynomial. + */ +@UnstableKMathAPI +public fun Polynomial.derivative( + algebra: A, +): Polynomial where A : Ring, A : NumericAlgebra = 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 Polynomial.antiderivative( + algebra: A, +): Polynomial where A : Field, A : NumericAlgebra = 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 > Polynomial.integrate( + algebra: Field, + range: ClosedRange, +): C = algebra { + val integral = antiderivative(algebra) + integral.substitute(algebra, range.endInclusive) - integral.substitute(algebra, range.start) +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt index eb88d9ae0..0fcd4c6e5 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt +++ b/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 > PiecewisePolynomial.integrate(algebra: Field): PiecewisePolynomial = - 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] diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt index 2266092a3..62819be0c 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt +++ b/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> : Interpolator): PiecewisePolynomial override fun interpolate(points: XYColumnarData): (T) -> T = { x -> - interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue() + interpolatePolynomials(points).substitute(algebra, x) ?: getDefaultValue() } } diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/PolynomialTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/PolynomialTest.kt index 05c16d17e..e0f0e32a4 100644 --- a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/PolynomialTest.kt +++ b/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) } } \ No newline at end of file