From 2483c56f1c95180e168bd564be6131122a6c8cfe Mon Sep 17 00:00:00 2001 From: Gleb Minaev <43728100+lounres@users.noreply.github.com> Date: Thu, 3 Mar 2022 20:45:35 +0300 Subject: [PATCH] Restructured Polynomial --- .../kmath/functions/AbstractPolynomial.kt | 354 +++++++++++++++ .../kscience/kmath/functions/Piecewise.kt | 8 +- .../kscience/kmath/functions/Polynomial.kt | 425 ++++++++++++++---- .../kmath/functions/polynomialUtil.kt | 139 ++++++ .../kmath/integration/SplineIntegrator.kt | 3 +- .../kmath/interpolation/Interpolator.kt | 4 +- .../kmath/functions/PolynomialTest.kt | 11 +- kotlin-js-store/yarn.lock | 56 +-- 8 files changed, 846 insertions(+), 154 deletions(-) create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/AbstractPolynomial.kt create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt 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 15d548641..2e6873ece 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 b13adefa5..bbd76c87e 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 @@ -9,7 +9,7 @@ package space.kscience.kmath.interpolation import space.kscience.kmath.data.XYColumnarData import space.kscience.kmath.functions.PiecewisePolynomial -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 @@ -33,7 +33,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 diff --git a/kotlin-js-store/yarn.lock b/kotlin-js-store/yarn.lock index e21abe604..9fc75720e 100644 --- a/kotlin-js-store/yarn.lock +++ b/kotlin-js-store/yarn.lock @@ -274,11 +274,6 @@ argparse@^2.0.1: resolved "https://registry.yarnpkg.com/argparse/-/argparse-2.0.1.tgz#246f50f3ca78a3240f6c997e8a9bd1eac49e4b38" integrity sha512-8+9WqebbFzpX9OR+Wa6O29asIogeRMzcGtAINdpMHHyAg10f05aSFVBbcEqGf/PXw1EjAZ+q2/bEBg3DvurK3Q== -astring@1.7.5: - version "1.7.5" - resolved "https://registry.yarnpkg.com/astring/-/astring-1.7.5.tgz#a7d47fceaf32b052d33a3d07c511efeec67447ca" - integrity sha512-lobf6RWXb8c4uZ7Mdq0U12efYmpD1UFnyOWVJPTa3ukqZrMopav+2hdNu0hgBF0JIBFK9QgrBDfwYvh3DFJDAA== - balanced-match@^1.0.0: version "1.0.2" resolved "https://registry.yarnpkg.com/balanced-match/-/balanced-match-1.0.2.tgz#e83e3a7e3f300b34cb9d87f615fa0cbf357690ee" @@ -294,24 +289,11 @@ base64id@2.0.0, base64id@~2.0.0: resolved "https://registry.yarnpkg.com/base64id/-/base64id-2.0.0.tgz#2770ac6bc47d312af97a8bf9a634342e0cd25cb6" integrity sha512-lGe34o6EHj9y3Kts9R4ZYs/Gr+6N7MCaMlIFA3F1R2O5/m7K06AxfSeO5530PEERE6/WyEg3lsuyw4GHlPZHog== -benchmark@*: - version "2.1.4" - resolved "https://registry.yarnpkg.com/benchmark/-/benchmark-2.1.4.tgz#09f3de31c916425d498cc2ee565a0ebf3c2a5629" - integrity sha1-CfPeMckWQl1JjMLuVloOvzwqVik= - dependencies: - lodash "^4.17.4" - platform "^1.3.3" - binary-extensions@^2.0.0: version "2.2.0" resolved "https://registry.yarnpkg.com/binary-extensions/-/binary-extensions-2.2.0.tgz#75f502eeaf9ffde42fc98829645be4ea76bd9e2d" integrity sha512-jDctJ/IVQbZoJykoeHbhXpOlNBqGNcwXJKJog42E5HDPUwQTSdjCHdihjj0DlnheQ7blbT6dHOafNAiS8ooQKA== -binaryen@101.0.0: - version "101.0.0" - resolved "https://registry.yarnpkg.com/binaryen/-/binaryen-101.0.0.tgz#42a9e4cc7a22e2c1d75a31d28005a9b518b2c555" - integrity sha512-FRmVxvrR8jtcf0qcukNAPZDM3dZ2sc9GmA/hKxBI7k3fFzREKh1cAs+ruQi+ITTKz7u/AuFMuVnbJwTh0ef/HQ== - body-parser@^1.19.0: version "1.19.1" resolved "https://registry.yarnpkg.com/body-parser/-/body-parser-1.19.1.tgz#1499abbaa9274af3ecc9f6f10396c995943e31d4" @@ -599,14 +581,6 @@ dom-serialize@^2.2.1: extend "^3.0.0" void-elements "^2.0.0" -dukat@0.5.8-rc.4: - version "0.5.8-rc.4" - resolved "https://registry.yarnpkg.com/dukat/-/dukat-0.5.8-rc.4.tgz#90384dcb50b14c26f0e99dae92b2dea44f5fce21" - integrity sha512-ZnMt6DGBjlVgK2uQamXfd7uP/AxH7RqI0BL9GLrrJb2gKdDxvJChWy+M9AQEaL+7/6TmxzJxFOsRiInY9oGWTA== - dependencies: - google-protobuf "3.12.2" - typescript "3.9.5" - ee-first@1.1.1: version "1.1.1" resolved "https://registry.yarnpkg.com/ee-first/-/ee-first-1.1.1.tgz#590c61156b0ae2f4f0255732a158b266bc56b21d" @@ -881,11 +855,6 @@ glob@^7.1.3, glob@^7.1.7: once "^1.3.0" path-is-absolute "^1.0.0" -google-protobuf@3.12.2: - version "3.12.2" - resolved "https://registry.yarnpkg.com/google-protobuf/-/google-protobuf-3.12.2.tgz#50ce9f9b6281235724eb243d6a83e969a2176e53" - integrity sha512-4CZhpuRr1d6HjlyrxoXoocoGFnRYgKULgMtikMddA9ztRyYR59Aondv2FioyxWVamRo0rF2XpYawkTCBEQOSkA== - graceful-fs@^4.1.2, graceful-fs@^4.1.6, graceful-fs@^4.2.0, graceful-fs@^4.2.4, graceful-fs@^4.2.6: version "4.2.9" resolved "https://registry.yarnpkg.com/graceful-fs/-/graceful-fs-4.2.9.tgz#041b05df45755e587a24942279b9d113146e1c96" @@ -1065,11 +1034,6 @@ jest-worker@^27.4.1: merge-stream "^2.0.0" supports-color "^8.0.0" -js-base64@3.6.1: - version "3.6.1" - resolved "https://registry.yarnpkg.com/js-base64/-/js-base64-3.6.1.tgz#555aae398b74694b4037af1f8a5a6209d170efbe" - integrity sha512-Frdq2+tRRGLQUIQOgsIGSCd1VePCS2fsddTG5dTCqR0JHgltXWfsxnY0gIXPoMeRmdom6Oyq+UMOFg5suduOjQ== - js-yaml@4.1.0: version "4.1.0" resolved "https://registry.yarnpkg.com/js-yaml/-/js-yaml-4.1.0.tgz#c1fb65f8f5017901cdd2c951864ba18458a10602" @@ -1179,7 +1143,7 @@ locate-path@^6.0.0: dependencies: p-locate "^5.0.0" -lodash@^4.17.15, lodash@^4.17.21, lodash@^4.17.4: +lodash@^4.17.15, lodash@^4.17.21: version "4.17.21" resolved "https://registry.yarnpkg.com/lodash/-/lodash-4.17.21.tgz#679591c564c3bffaae8454cf0b3df370c3d6911c" integrity sha512-v2kDEe57lecTulaDIuNTPy3Ry4gLGJ6Z1O3vE1krgXZNrsQ+LFTGHVxVjcXPs17LhbZVGedAJv8XZ1tvj5FvSg== @@ -1437,11 +1401,6 @@ pkg-dir@^4.2.0: dependencies: find-up "^4.0.0" -platform@^1.3.3: - version "1.3.6" - resolved "https://registry.yarnpkg.com/platform/-/platform-1.3.6.tgz#48b4ce983164b209c2d45a107adb31f473a6e7a7" - integrity sha512-fnWVljUchTro6RiCFvCXBbNhJc2NijN7oIQxbwsyL0buWJPG85v81ehlHI9fXrJsMNgTofEoWIQeClKpgxFLrg== - postcss-modules-extract-imports@^3.0.0: version "3.0.0" resolved "https://registry.yarnpkg.com/postcss-modules-extract-imports/-/postcss-modules-extract-imports-3.0.0.tgz#cda1f047c0ae80c97dbe28c3e76a43b88025741d" @@ -1696,14 +1655,6 @@ source-map-loader@3.0.0: iconv-lite "^0.6.2" source-map-js "^0.6.2" -source-map-support@0.5.20: - version "0.5.20" - resolved "https://registry.yarnpkg.com/source-map-support/-/source-map-support-0.5.20.tgz#12166089f8f5e5e8c56926b377633392dd2cb6c9" - integrity sha512-n1lZZ8Ve4ksRqizaBQgxXDgKwttHDhyfQjA6YZZn8+AroHbsIz+JjwxQDxbp+7y5OYCI8t1Yk7etjD9CRd2hIw== - dependencies: - buffer-from "^1.0.0" - source-map "^0.6.0" - source-map-support@~0.5.20: version "0.5.21" resolved "https://registry.yarnpkg.com/source-map-support/-/source-map-support-0.5.21.tgz#04fe7c7f9e1ed2d662233c28cb2b35b9f63f6e4f" @@ -1838,11 +1789,6 @@ type-is@~1.6.18: media-typer "0.3.0" mime-types "~2.1.24" -typescript@3.9.5: - version "3.9.5" - resolved "https://registry.yarnpkg.com/typescript/-/typescript-3.9.5.tgz#586f0dba300cde8be52dd1ac4f7e1009c1b13f36" - integrity sha512-hSAifV3k+i6lEoCJ2k6R2Z/rp/H3+8sdmcn5NrS3/3kE7+RyZXm9aqvxWqjEXHAd8b0pShatpcdMTvEdvAJltQ== - ua-parser-js@^0.7.28: version "0.7.31" resolved "https://registry.yarnpkg.com/ua-parser-js/-/ua-parser-js-0.7.31.tgz#649a656b191dffab4f21d5e053e27ca17cbff5c6"