From 843d63c76a46a55c683f4a8276663f5fecc6d6f5 Mon Sep 17 00:00:00 2001 From: Gleb Minaev <43728100+lounres@users.noreply.github.com> Date: Sun, 13 Mar 2022 03:27:00 +0300 Subject: [PATCH] Added support for all polynomials. But standard utilities still are not fully implemented. --- .../kmath/functions/AbstractPolynomial.kt | 47 +- .../functions/AbstractRationalFunction.kt | 2 - .../kmath/functions/LabeledPolynomial.kt | 927 ++++++++++++++++++ .../kmath/functions/NumberedPolynomial.kt | 689 +++++++++++++ .../kscience/kmath/functions/Polynomial.kt | 151 ++- .../kscience/kmath/functions/Variable.kt | 38 + .../kscience/kmath/functions/algebraicStub.kt | 51 + .../kmath/functions/labeledPolynomialUtil.kt | 490 +++++++++ .../kmath/functions/numberedPolynomialUtil.kt | 605 ++++++++++++ .../kmath/functions/polynomialUtil.kt | 6 +- 10 files changed, 2911 insertions(+), 95 deletions(-) create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/NumberedPolynomial.kt create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Variable.kt create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/algebraicStub.kt create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledPolynomialUtil.kt create mode 100644 kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedPolynomialUtil.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 index 237a72bcc..b7b7116f0 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/AbstractPolynomial.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/AbstractPolynomial.kt @@ -17,6 +17,9 @@ public interface AbstractPolynomial /** * Abstraction of ring of polynomials of type [P] over ring of constants of type [C]. + * + * @param C the type of constants. Polynomials have them as a coefficients in their terms. + * @param P the type of polynomials. */ @Suppress("INAPPLICABLE_JVM_NAME", "PARAMETER_NAME_CHANGED_ON_OVERRIDE") public interface AbstractPolynomialSpace> : Ring

{ @@ -307,48 +310,4 @@ public interface AbstractPolynomialSpace> : Ring

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/AbstractRationalFunction.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/AbstractRationalFunction.kt index 5cb570c9f..34050aa0f 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/AbstractRationalFunction.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/AbstractRationalFunction.kt @@ -5,8 +5,6 @@ package space.kscience.kmath.functions -import space.kscience.kmath.functions.AbstractPolynomialSpace.Companion.optimizedAddMultiplied -import space.kscience.kmath.functions.AbstractPolynomialSpace.Companion.optimizedMultiply import space.kscience.kmath.operations.* import kotlin.js.JsName import kotlin.jvm.JvmName diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt new file mode 100644 index 000000000..48f6f57fa --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt @@ -0,0 +1,927 @@ +package space.kscience.kmath.functions + +import space.kscience.kmath.operations.* +import kotlin.contracts.InvocationKind +import kotlin.contracts.contract +import kotlin.experimental.ExperimentalTypeInference +import kotlin.jvm.JvmName +import kotlin.math.max + + +/** + * Represents multivariate polynomials with labeled variables. + * + * @param C Ring in which the polynomial is considered. + */ +public class LabeledPolynomial +internal constructor( + /** + * Map that collects coefficients of the polynomial. Every non-zero monomial + * `a x_1^{d_1} ... x_n^{d_n}` is represented as pair "key-value" in the map, where value is coefficients `a` and + * key is map that associates variables in the monomial with multiplicity of them occurring in the monomial. + * For example polynomial + * ``` + * 5 a^2 c^3 - 6 b + 0 b c + * ``` + * has coefficients represented as + * ``` + * mapOf( + * mapOf( + * a to 2, + * c to 3 + * ) to 5, + * mapOf( + * b to 1 + * ) to (-6) + * ) + * ``` + * where `a`, `b` and `c` are corresponding [Variable] objects. + */ + public val coefficients: Map, C> +) : AbstractPolynomial { + override fun toString(): String = "LabeledPolynomial$coefficients" +} + +// region Internal utilities + +/** + * Represents internal [LabeledPolynomial] errors. + */ +internal class LabeledPolynomialError: Error { + constructor(): super() + constructor(message: String): super(message) + constructor(message: String?, cause: Throwable?): super(message, cause) + constructor(cause: Throwable?): super(cause) +} + +/** + * Throws an [LabeledPolynomialError] with the given [message]. + */ +internal fun labeledPolynomialError(message: Any): Nothing = throw LabeledPolynomialError(message.toString()) + +/** + * Returns the same degrees description of the monomial, but without zero degrees. + */ +internal fun Map.cleanUp() = filterValues { it > 0U } + +// endregion + +// region Constructors and converters + +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * @param toCheckInput If it's `true` cleaning of [coefficients] is executed otherwise it is not. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//internal fun > LabeledPolynomial(coefs: Map, C>, toCheckInput: Boolean = true): LabeledPolynomial { +// if (!toCheckInput) return LabeledPolynomial(coefs) +// +// // Map for cleaned coefficients. +// val fixedCoefs = mutableMapOf, C>() +// +// // Cleaning the degrees, summing monomials of the same degrees. +// for (entry in coefs) { +// val key = entry.key.cleanUp() +// val value = entry.value +// fixedCoefs[key] = if (key in fixedCoefs) fixedCoefs[key]!! + value else value +// } +// +// // Removing zero monomials. +// return LabeledPolynomial( +// fixedCoefs +// .filter { it.value.isNotZero() } +// ) +//} +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * @param toCheckInput If it's `true` cleaning of [coefficients] is executed otherwise it is not. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//internal fun > LabeledPolynomial(pairs: Collection, C>>, toCheckInput: Boolean): LabeledPolynomial { +// if (!toCheckInput) return LabeledPolynomial(pairs.toMap()) +// +// // Map for cleaned coefficients. +// val fixedCoefs = mutableMapOf, C>() +// +// // Cleaning the degrees, summing monomials of the same degrees. +// for (entry in pairs) { +// val key = entry.first.cleanUp() +// val value = entry.second +// fixedCoefs[key] = if (key in fixedCoefs) fixedCoefs[key]!! + value else value +// } +// +// // Removing zero monomials. +// return LabeledPolynomial( +// fixedCoefs.filterValues { it.isNotZero() } +// ) +//} +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param pairs Collection of pairs that represents monomials. +// * @param toCheckInput If it's `true` cleaning of [coefficients] is executed otherwise it is not. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//internal fun > LabeledPolynomial(vararg pairs: Pair, C>, toCheckInput: Boolean): LabeledPolynomial { +// if (!toCheckInput) return LabeledPolynomial(pairs.toMap()) +// +// // Map for cleaned coefficients. +// val fixedCoefs = mutableMapOf, C>() +// +// // Cleaning the degrees, summing monomials of the same degrees. +// for (entry in pairs) { +// val key = entry.first.cleanUp() +// val value = entry.second +// fixedCoefs[key] = if (key in fixedCoefs) fixedCoefs[key]!! + value else value +// } +// +// // Removing zero monomials. +// return LabeledPolynomial( +// fixedCoefs.filterValues { it.isNotZero() } +// ) +//} +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param coefs Coefficients of the instants. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//fun > LabeledPolynomial(coefs: Map, C>): LabeledPolynomial = LabeledPolynomial(coefs, toCheckInput = true) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param pairs Collection of pairs that represents monomials. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//fun > LabeledPolynomial(pairs: Collection, C>>): LabeledPolynomial = LabeledPolynomial(pairs, toCheckInput = true) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param pairs Collection of pairs that represents monomials. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//fun > LabeledPolynomial(vararg pairs: Pair, C>): LabeledPolynomial = LabeledPolynomial(*pairs, toCheckInput = true) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * @param toCheckInput If it's `true` cleaning of [coefficients] is executed otherwise it is not. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(LabeledPolynomialSpace) +//internal fun > LabeledPolynomial(coefs: Map, C>, toCheckInput: Boolean = true): LabeledPolynomial { +// if (!toCheckInput) return LabeledPolynomial(coefs) +// +// // Map for cleaned coefficients. +// val fixedCoefs = mutableMapOf, C>() +// +// // Cleaning the degrees, summing monomials of the same degrees. +// for (entry in coefs) { +// val key = entry.key.cleanUp() +// val value = entry.value +// fixedCoefs[key] = if (key in fixedCoefs) fixedCoefs[key]!! + value else value +// } +// +// // Removing zero monomials. +// return LabeledPolynomial( +// fixedCoefs +// .filter { it.value.isNotZero() } +// ) +//} +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * @param toCheckInput If it's `true` cleaning of [coefficients] is executed otherwise it is not. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(LabeledPolynomialSpace) +//internal fun > LabeledPolynomial(pairs: Collection, C>>, toCheckInput: Boolean): LabeledPolynomial { +// if (!toCheckInput) return LabeledPolynomial(pairs.toMap()) +// +// // Map for cleaned coefficients. +// val fixedCoefs = mutableMapOf, C>() +// +// // Cleaning the degrees, summing monomials of the same degrees. +// for (entry in pairs) { +// val key = entry.first.cleanUp() +// val value = entry.second +// fixedCoefs[key] = if (key in fixedCoefs) fixedCoefs[key]!! + value else value +// } +// +// // Removing zero monomials. +// return LabeledPolynomial( +// fixedCoefs.filterValues { it.isNotZero() } +// ) +//} +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param pairs Collection of pairs that represents monomials. +// * @param toCheckInput If it's `true` cleaning of [coefficients] is executed otherwise it is not. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(LabeledPolynomialSpace) +//internal fun > LabeledPolynomial(vararg pairs: Pair, C>, toCheckInput: Boolean): LabeledPolynomial { +// if (!toCheckInput) return LabeledPolynomial(pairs.toMap()) +// +// // Map for cleaned coefficients. +// val fixedCoefs = mutableMapOf, C>() +// +// // Cleaning the degrees, summing monomials of the same degrees. +// for (entry in pairs) { +// val key = entry.first.cleanUp() +// val value = entry.second +// fixedCoefs[key] = if (key in fixedCoefs) fixedCoefs[key]!! + value else value +// } +// +// // Removing zero monomials. +// return LabeledPolynomial( +// fixedCoefs.filterValues { it.isNotZero() } +// ) +//} +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param coefs Coefficients of the instants. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial(coefs: Map, C>): LabeledPolynomial = LabeledPolynomial(coefs, toCheckInput = true) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param pairs Collection of pairs that represents monomials. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial(pairs: Collection, C>>): LabeledPolynomial = LabeledPolynomial(pairs, toCheckInput = true) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from keys of received +// * map, sums up proportional monomials, removes aero monomials, and if result is zero map adds only element in it. +// * +// * @param pairs Collection of pairs that represents monomials. +// * +// * @throws LabeledPolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial(vararg pairs: Pair, C>): LabeledPolynomial = LabeledPolynomial(*pairs, toCheckInput = true) +// +//fun C.asLabeledPolynomial() : LabeledPolynomial = LabeledPolynomial(mapOf(emptyMap() to this)) +// +//context(A) +//fun > Variable.asLabeledPolynomial() : LabeledPolynomial = LabeledPolynomial(mapOf(mapOf(this to 1U) to one)) +// +//context(LabeledPolynomialSpace) +//fun > Variable.asLabeledPolynomial() : LabeledPolynomial = LabeledPolynomial(mapOf(mapOf(this to 1U) to ring.one)) +// +//context(A) +//fun > Variable.asLabeledPolynomial(c: C) : LabeledPolynomial = +// if(c.isZero()) LabeledPolynomial(emptyMap()) +// else LabeledPolynomial(mapOf(mapOf(this to 1U) to c)) +// +//context(LabeledPolynomialSpace) +//fun > Variable.asLabeledPolynomial(c: C) : LabeledPolynomial = +// if(c.isZero()) zero +// else LabeledPolynomial(mapOf(mapOf(this to 1U) to c)) + +// endregion + +/** + * Space of polynomials. + * + * @param C the type of operated polynomials. + * @param A the intersection of [Ring] of [C] and [ScaleOperations] of [C]. + * @param ring the [A] instance. + */ +@Suppress("PARAMETER_NAME_CHANGED_ON_OVERRIDE", "INAPPLICABLE_JVM_NAME") +public class LabeledPolynomialSpace>( + 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 + + // 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 Variable-integer relation + public operator fun Variable.plus(other: Int): LabeledPolynomial = + if (other == 0) LabeledPolynomial(mapOf( + mapOf(this@plus to 1U) to ring.one, + )) + else LabeledPolynomial(mapOf( + mapOf(this@plus to 1U) to ring.one, + emptyMap() to ring.one * other, + )) + public operator fun Variable.minus(other: Int): LabeledPolynomial = + if (other == 0) LabeledPolynomial(mapOf( + mapOf(this@minus to 1U) to -ring.one, + )) + else LabeledPolynomial(mapOf( + mapOf(this@minus to 1U) to -ring.one, + emptyMap() to ring.one * other, + )) + public operator fun Variable.times(other: Int): LabeledPolynomial = + if (other == 0) zero + else LabeledPolynomial(mapOf( + mapOf(this to 1U) to ring.one * other, + )) + // endregion + + // region Integer-variable relation + public operator fun Int.plus(other: Variable): LabeledPolynomial = + if (this == 0) LabeledPolynomial(mapOf( + mapOf(other to 1U) to ring.one, + )) + else LabeledPolynomial(mapOf( + mapOf(other to 1U) to ring.one, + emptyMap() to ring.one * this@plus, + )) + public operator fun Int.minus(other: Variable): LabeledPolynomial = + if (this == 0) LabeledPolynomial(mapOf( + mapOf(other to 1U) to -ring.one, + )) + else LabeledPolynomial(mapOf( + mapOf(other to 1U) to -ring.one, + emptyMap() to ring.one * this@minus, + )) + public operator fun Int.times(other: Variable): LabeledPolynomial = + if (this == 0) zero + else LabeledPolynomial(mapOf( + mapOf(other to 1U) to ring.one * this@times, + )) + // endregion + + // region Polynomial-integer relation + public override operator fun LabeledPolynomial.plus(other: Int): LabeledPolynomial = + if (other == 0) this + else + LabeledPolynomial( + coefficients + .toMutableMap() + .apply { + val degs = emptyMap() + + val result = getOrElse(degs) { ring.zero } + other + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + public override operator fun LabeledPolynomial.minus(other: Int): LabeledPolynomial = + if (other == 0) this + else + LabeledPolynomial( + coefficients + .toMutableMap() + .apply { + val degs = emptyMap() + + val result = getOrElse(degs) { ring.zero } - other + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + public override operator fun LabeledPolynomial.times(other: Int): LabeledPolynomial = + if (other == 0) zero + else LabeledPolynomial( + coefficients + .applyAndRemoveZeros { + mapValues { (_, c) -> c * other } + } + ) + // endregion + + // region Integer-polynomial relation + public override operator fun Int.plus(other: LabeledPolynomial): LabeledPolynomial = + if (this == 0) other + else + LabeledPolynomial( + other.coefficients + .toMutableMap() + .apply { + val degs = emptyMap() + + val result = this@plus + getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + public override operator fun Int.minus(other: LabeledPolynomial): LabeledPolynomial = + if (this == 0) other + else + LabeledPolynomial( + other.coefficients + .toMutableMap() + .apply { + val degs = emptyMap() + + val result = this@minus - getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + public override operator fun Int.times(other: LabeledPolynomial): LabeledPolynomial = + if (this == 0) zero + else LabeledPolynomial( + other.coefficients + .applyAndRemoveZeros { + mapValues { (_, c) -> this@times * c } + } + ) + // endregion + + // region Constant-constant relation + @JvmName("constantUnaryMinus") + override operator fun C.unaryMinus(): C = ring { -this@unaryMinus } + @JvmName("constantPlus") + override operator fun C.plus(other: C): C = ring { this@plus + other } + @JvmName("constantMinus") + override operator fun C.minus(other: C): C = ring { this@minus - other } + @JvmName("constantTimes") + 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-variable relation + public operator fun C.plus(other: Variable): LabeledPolynomial = + if (isZero()) LabeledPolynomial(mapOf( + mapOf(other to 1U) to ring.one, + )) + else LabeledPolynomial(mapOf( + mapOf(other to 1U) to ring.one, + emptyMap() to this@plus, + )) + public operator fun C.minus(other: Variable): LabeledPolynomial = + if (isZero()) LabeledPolynomial(mapOf( + mapOf(other to 1U) to -ring.one, + )) + else LabeledPolynomial(mapOf( + mapOf(other to 1U) to -ring.one, + emptyMap() to this@minus, + )) + public operator fun C.times(other: Variable): LabeledPolynomial = + if (isZero()) zero + else LabeledPolynomial(mapOf( + mapOf(other to 1U) to this@times, + )) + // endregion + + // region Variable-constant relation + public operator fun Variable.plus(other: C): LabeledPolynomial = + if (other.isZero()) LabeledPolynomial(mapOf( + mapOf(this@plus to 1U) to ring.one, + )) + else LabeledPolynomial(mapOf( + mapOf(this@plus to 1U) to ring.one, + emptyMap() to other, + )) + public operator fun Variable.minus(other: C): LabeledPolynomial = + if (other.isZero()) LabeledPolynomial(mapOf( + mapOf(this@minus to 1U) to -ring.one, + )) + else LabeledPolynomial(mapOf( + mapOf(this@minus to 1U) to -ring.one, + emptyMap() to other, + )) + public operator fun Variable.times(other: C): LabeledPolynomial = + if (other.isZero()) zero + else LabeledPolynomial(mapOf( + mapOf(this@times to 1U) to other, + )) + // endregion + + // region Constant-polynomial relation + override operator fun C.plus(other: LabeledPolynomial): LabeledPolynomial = + if (this.isZero()) other + else with(other.coefficients) { + if (isEmpty()) LabeledPolynomial(mapOf(emptyMap() to this@plus)) + else LabeledPolynomial( + toMutableMap() + .apply { + val degs = emptyMap() + + val result = this@plus + getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + override operator fun C.minus(other: LabeledPolynomial): LabeledPolynomial = + if (this.isZero()) other + else with(other.coefficients) { + if (isEmpty()) LabeledPolynomial(mapOf(emptyMap() to this@minus)) + else LabeledPolynomial( + toMutableMap() + .apply { + forEach { (degs, c) -> if(degs.isNotEmpty()) this[degs] = -c } + + val degs = emptyMap() + + val result = this@minus - getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + override operator fun C.times(other: LabeledPolynomial): LabeledPolynomial = + if (this.isZero()) zero + else LabeledPolynomial( + other.coefficients + .applyAndRemoveZeros { + mapValues { (_, c) -> this@times * c } + } + ) + // endregion + + // region Polynomial-constant relation + /** + * Returns sum of the polynomials. [other] is interpreted as [UnivariatePolynomial]. + */ + override operator fun LabeledPolynomial.plus(other: C): LabeledPolynomial = + if (other.isZero()) this + else with(coefficients) { + if (isEmpty()) LabeledPolynomial(mapOf(emptyMap() to other)) + else LabeledPolynomial( + toMutableMap() + .apply { + val degs = emptyMap() + + val result = getOrElse(degs) { ring.zero } + other + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + /** + * Returns difference of the polynomials. [other] is interpreted as [UnivariatePolynomial]. + */ + override operator fun LabeledPolynomial.minus(other: C): LabeledPolynomial = + if (other.isZero()) this + else with(coefficients) { + if (isEmpty()) LabeledPolynomial(mapOf(emptyMap() to other)) + else LabeledPolynomial( + toMutableMap() + .apply { + forEach { (degs, c) -> if(degs.isNotEmpty()) this[degs] = -c } + + val degs = emptyMap() + + val result = getOrElse(degs) { ring.zero } - other + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + /** + * Returns product of the polynomials. [other] is interpreted as [UnivariatePolynomial]. + */ + override operator fun LabeledPolynomial.times(other: C): LabeledPolynomial = + if (other.isZero()) zero + else LabeledPolynomial( + coefficients + .applyAndRemoveZeros { + mapValues { (_, c) -> c * other } + } + ) + // endregion + + // region Variable-variable relation + public operator fun Variable.plus(other: Variable): LabeledPolynomial = + if (this == other) LabeledPolynomial(mapOf( + mapOf(this to 1U) to ring.one * 2 + )) + else LabeledPolynomial(mapOf( + mapOf(this to 1U) to ring.one, + mapOf(other to 1U) to ring.one, + )) + public operator fun Variable.minus(other: Variable): LabeledPolynomial = + if (this == other) zero + else LabeledPolynomial(mapOf( + mapOf(this to 1U) to ring.one, + mapOf(other to 1U) to -ring.one, + )) + public operator fun Variable.times(other: Variable): LabeledPolynomial = + if (this == other) LabeledPolynomial(mapOf( + mapOf(this to 2U) to ring.one + )) + else LabeledPolynomial(mapOf( + mapOf(this to 1U, other to 1U) to ring.one, + )) + // endregion + + // region Variable-polynomial relation + public operator fun Variable.plus(other: LabeledPolynomial): LabeledPolynomial = + with(other.coefficients) { + if (isEmpty()) LabeledPolynomial(mapOf(mapOf(this@plus to 1u) to ring.one)) + else LabeledPolynomial( + toMutableMap() + .apply { + val degs = mapOf(this@plus to 1U) + + val result = ring.one + getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + public operator fun Variable.minus(other: LabeledPolynomial): LabeledPolynomial = + with(other.coefficients) { + if (isEmpty()) LabeledPolynomial(mapOf(mapOf(this@minus to 1u) to ring.one)) + else LabeledPolynomial( + toMutableMap() + .apply { + forEach { (degs, c) -> if(degs.isNotEmpty()) this[degs] = -c } + + val degs = mapOf(this@minus to 1U) + + val result = ring.one - getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + public operator fun Variable.times(other: LabeledPolynomial): LabeledPolynomial = + LabeledPolynomial( + other.coefficients + .mapKeys { (degs, _) -> degs.toMutableMap().also{ it[this] = if (this in it) it[this]!! + 1U else 1U } } + ) + // endregion + + // region Polynomial-variable relation + public operator fun LabeledPolynomial.plus(other: Variable): LabeledPolynomial = + with(coefficients) { + if (isEmpty()) LabeledPolynomial(mapOf(mapOf(other to 1u) to ring.one)) + else LabeledPolynomial( + toMutableMap() + .apply { + val degs = mapOf(other to 1U) + + val result = ring.one + getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + public operator fun LabeledPolynomial.minus(other: Variable): LabeledPolynomial = + with(coefficients) { + if (isEmpty()) LabeledPolynomial(mapOf(mapOf(other to 1u) to ring.one)) + else LabeledPolynomial( + toMutableMap() + .apply { + val degs = mapOf(other to 1U) + + val result = ring.one - getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + public operator fun LabeledPolynomial.times(other: Variable): LabeledPolynomial = + LabeledPolynomial( + coefficients + .mapKeys { (degs, _) -> degs.toMutableMap().also{ it[other] = if (other in it) it[other]!! + 1U else 1U } } + ) + // endregion + + // region Polynomial-polynomial relation + /** + * Returns negation of the polynomial. + */ + override fun LabeledPolynomial.unaryMinus(): LabeledPolynomial = + LabeledPolynomial( + coefficients.mapValues { -it.value } + ) + /** + * Returns sum of the polynomials. + */ + override operator fun LabeledPolynomial.plus(other: LabeledPolynomial): LabeledPolynomial = + LabeledPolynomial( + coefficients + .applyAndRemoveZeros { + other.coefficients + .mapValuesTo(this) { (key, value) -> if (key in this) this[key]!! + value else value } + } + ) + /** + * Returns difference of the polynomials. + */ + override operator fun LabeledPolynomial.minus(other: LabeledPolynomial): LabeledPolynomial = + LabeledPolynomial( + coefficients + .applyAndRemoveZeros { + other.coefficients + .mapValuesTo(this) { (key, value) -> if (key in this) this[key]!! - value else -value } + } + ) + /** + * Returns product of the polynomials. + */ + override operator fun LabeledPolynomial.times(other: LabeledPolynomial): LabeledPolynomial = + when { + isZero() -> zero + other.isZero() -> zero + else -> LabeledPolynomial( + buildCoefficients { + for ((degs1, c1) in coefficients) for ((degs2, c2) in other.coefficients) { + val degs = degs1.toMutableMap() + degs2.mapValuesTo(degs) { (variable, deg) -> degs.getOrElse(variable) { 0u } + deg } + val c = c1 * c2 + this[degs] = if (degs in this) this[degs]!! + c else c + } + } + ) + } + + override val zero: LabeledPolynomial = LabeledPolynomial(mapOf(emptyMap() to ring.zero)) + override val one: LabeledPolynomial = LabeledPolynomial(mapOf(emptyMap() to ring.one)) + + // TODO: Docs + @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "CovariantEquals") + override fun LabeledPolynomial.equals(other: LabeledPolynomial): Boolean = + when { + this === other -> true + else -> coefficients.size == other.coefficients.size && + coefficients.all { (key, value) -> with(other.coefficients) { key in this && this[key] == value } } + } + // 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. + */ + override val LabeledPolynomial.degree: Int + get() = coefficients.entries.maxOfOrNull { (degs, c) -> if (c.isZero()) -1 else degs.values.sum().toInt() } ?: -1 + /** + * Map that associates variables (that appear in the polynomial in positive exponents) with their most exponents + * in which they are appeared in the polynomial. + * + * As consequence all values in the map are positive integers. Also, if the polynomial is constant, the map is empty. + * And keys of the map is the same as in [variables]. + */ + public val LabeledPolynomial.degrees: Map + get() = + buildMap { + coefficients.entries.forEach { (degs, c) -> + if (c.isNotZero()) degs.mapValuesTo(this) { (variable, deg) -> + max(getOrElse(variable) { 0u }, deg) + } + } + } + /** + * Set of all variables that appear in the polynomial in positive exponents. + */ + public val LabeledPolynomial.variables: Set + get() = + buildSet { + coefficients.entries.forEach { (degs, c) -> if (c.isNotZero()) addAll(degs.keys) } + } + /** + * Count of all variables that appear in the polynomial in positive exponents. + */ + public val LabeledPolynomial.countOfVariables: Int get() = variables.size + + /** + * Checks if the instant is constant polynomial (of degree no more than 0) over considered ring. + */ + override fun LabeledPolynomial.isConstant(): Boolean = + coefficients.all { (degs, c) -> degs.isEmpty() || c.isZero() } + /** + * Checks if the instant is constant non-zero polynomial (of degree no more than 0) over considered ring. + */ + override fun LabeledPolynomial.isNonZeroConstant(): Boolean = + with(coefficients) { + var foundAbsoluteTermAndItIsNotZero = false + for ((degs, c) in this) { + if (degs.isNotEmpty()) if (c.isNotZero()) return@with false + else { + if (c.isZero()) return@with false + else foundAbsoluteTermAndItIsNotZero = true + } + } + foundAbsoluteTermAndItIsNotZero + } + + override fun LabeledPolynomial.asConstantOrNull(): C? = + with(coefficients) { + if(isConstant()) getOrElse(emptyMap()) { ring.zero } + else null + } + +// @Suppress("NOTHING_TO_INLINE") +// public inline fun LabeledPolynomial.substitute(argument: Map): LabeledPolynomial = this.substitute(ring, argument) +// @Suppress("NOTHING_TO_INLINE") +// @JvmName("substitutePolynomial") +// public inline fun LabeledPolynomial.substitute(argument: Map>): LabeledPolynomial = this.substitute(ring, argument) +// +// @Suppress("NOTHING_TO_INLINE") +// public inline fun LabeledPolynomial.asFunction(): (Map) -> LabeledPolynomial = { this.substitute(ring, it) } +// @Suppress("NOTHING_TO_INLINE") +// public inline fun LabeledPolynomial.asFunctionOnConstants(): (Map) -> LabeledPolynomial = { this.substitute(ring, it) } +// @Suppress("NOTHING_TO_INLINE") +// public inline fun LabeledPolynomial.asFunctionOnPolynomials(): (Map>) -> LabeledPolynomial = { this.substitute(ring, it) } +// +// @Suppress("NOTHING_TO_INLINE") +// public inline operator fun LabeledPolynomial.invoke(argument: Map): LabeledPolynomial = this.substitute(ring, argument) +// @Suppress("NOTHING_TO_INLINE") +// @JvmName("invokePolynomial") +// public inline operator fun LabeledPolynomial.invoke(argument: Map>): LabeledPolynomial = this.substitute(ring, argument) + // endregion + + // region Legacy + @Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE") + override inline fun add(left: LabeledPolynomial, right: LabeledPolynomial): LabeledPolynomial = left + right + @Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE") + override inline fun multiply(left: LabeledPolynomial, right: LabeledPolynomial): LabeledPolynomial = left * right + // endregion + + // region Utilities + // TODO: Move to region internal utilities with context receiver + internal fun MutableMap, C>.applyAndRemoveZeros(block: MutableMap, C>.() -> Unit) : MutableMap, C> { + contract { + callsInPlace(block, InvocationKind.EXACTLY_ONCE) + } + block() + for ((degs, c) in this) if (c.isZero()) this.remove(degs) + return this + } + internal fun Map, C>.applyAndRemoveZeros(block: MutableMap, C>.() -> Unit) : Map, C> = + toMutableMap().applyAndRemoveZeros(block) + @OptIn(ExperimentalTypeInference::class) + internal fun buildCoefficients(@BuilderInference builderAction: MutableMap, C>.() -> Unit): Map, C> { + contract { callsInPlace(builderAction, InvocationKind.EXACTLY_ONCE) } + return buildMap { + builderAction() + for ((degs, c) in this) if (c.isZero()) this.remove(degs) + } + } + // endregion +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/NumberedPolynomial.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/NumberedPolynomial.kt new file mode 100644 index 000000000..f1ad9a74f --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/NumberedPolynomial.kt @@ -0,0 +1,689 @@ +package space.kscience.kmath.functions + +import space.kscience.kmath.operations.* +import kotlin.contracts.InvocationKind +import kotlin.contracts.contract +import kotlin.experimental.ExperimentalTypeInference +import kotlin.jvm.JvmName +import kotlin.math.max + + +/** + * Polynomial model without fixation on specific context they are applied to. + * + * @param C the type of constants. + */ +public class NumberedPolynomial +internal constructor( + /** + * Map that collects coefficients of the polynomial. Every monomial `a x_1^{d_1} ... x_n^{d_n}` is represented as + * pair "key-value" in the map, where value is coefficients `a` and + * key is list that associates index of every variable in the monomial with multiplicity of the variable occurring + * in the monomial. For example coefficients of polynomial `5 x_1^2 x_3^3 - 6 x_2` can be represented as + * ``` + * mapOf( + * listOf(2, 0, 3) to 5, + * listOf(0, 1) to (-6), + * ) + * ``` + * and also as + * ``` + * mapOf( + * listOf(2, 0, 3) to 5, + * listOf(0, 1) to (-6), + * listOf(0, 1, 1) to 0, + * ) + * ``` + * It is recommended not to put zero monomials into the map, but is not prohibited. Lists of degrees always do not + * contain any zeros on end, but can contain zeros on start or anywhere in middle. + */ + public val coefficients: Map, C> +) : AbstractPolynomial { + override fun toString(): String = "NumberedPolynomial$coefficients" +} + +// region Internal utilities + +/** + * Represents internal [Polynomial] errors. + */ +internal class NumberedPolynomialError : Error { + constructor(): super() + constructor(message: String): super(message) + constructor(message: String?, cause: Throwable?): super(message, cause) + constructor(cause: Throwable?): super(cause) +} + +/** + * Throws an [PolynomialError] with the given [message]. + */ +internal fun numberedPolynomialError(message: Any): Nothing = throw PolynomialError(message.toString()) + +/** + * Returns the same degrees description of the monomial, but without extra zero degrees on the end. + */ +internal fun List.cleanUp() = subList(0, indexOfLast { it != 0U } + 1) + +// endregion + +// region Constructors and converters +// Waiting for context receivers :( TODO: Replace with context receivers when they will be available + +//context(A) +//@Suppress("FunctionName") +//internal fun > NumberedPolynomial(coefs: Map, C>, toCheckInput: Boolean): NumberedPolynomial { +// if (!toCheckInput) return NumberedPolynomial(coefs) +// +// val fixedCoefs = mutableMapOf, C>() +// +// for (entry in coefs) { +// val key = entry.key.cleanUp() +// val value = entry.value +// fixedCoefs[key] = if (key in fixedCoefs) fixedCoefs[key]!! + value else value +// } +// +// return NumberedPolynomial( +// fixedCoefs +// .filter { it.value.isNotZero() } +// ) +//} +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from end of received +// * lists, sums up proportional monomials, removes zero monomials, and if result is empty map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * @param toCheckInput If it's `true` cleaning of [coefficients] is executed otherwise it is not. +// * +// * @throws PolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//@Suppress("FunctionName") +//internal fun > NumberedPolynomial(pairs: Collection, C>>, toCheckInput: Boolean): NumberedPolynomial { +// if (!toCheckInput) return NumberedPolynomial(pairs.toMap()) +// +// val fixedCoefs = mutableMapOf, C>() +// +// for (entry in pairs) { +// val key = entry.first.cleanUp() +// val value = entry.second +// fixedCoefs[key] = if (key in fixedCoefs) fixedCoefs[key]!! + value else value +// } +// +// return NumberedPolynomial( +// fixedCoefs +// .filter { it.value.isNotZero() } +// ) +//} +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from end of received +// * lists, sums up proportional monomials, removes zero monomials, and if result is empty map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * @param toCheckInput If it's `true` cleaning of [coefficients] is executed otherwise it is not. +// * +// * @throws PolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//@Suppress("FunctionName") +//internal fun > NumberedPolynomial(vararg pairs: Pair, C>, toCheckInput: Boolean): NumberedPolynomial = +// NumberedPolynomial(pairs.toMap(), toCheckInput) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from end of received +// * lists, sums up proportional monomials, removes zero monomials, and if result is empty map adds only element in it. +// * +// * @param coefs Coefficients of the instants. +// * +// * @throws PolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//public fun > NumberedPolynomial(coefs: Map, C>) = NumberedPolynomial(coefs, toCheckInput = true) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from end of received +// * lists, sums up proportional monomials, removes zero monomials, and if result is empty map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * +// * @throws PolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//public fun > NumberedPolynomial(pairs: Collection, C>>) = NumberedPolynomial(pairs, toCheckInput = true) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from end of received +// * lists, sums up proportional monomials, removes zero monomials, and if result is empty map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * +// * @throws PolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(A) +//public fun > NumberedPolynomial(vararg pairs: Pair, C>) = NumberedPolynomial(*pairs, toCheckInput = true) +// +//context(NumberedPolynomialSpace) +//@Suppress("FunctionName") +//internal fun > NumberedPolynomial(coefs: Map, C>, toCheckInput: Boolean): NumberedPolynomial { +// if (!toCheckInput) return NumberedPolynomial(coefs) +// +// val fixedCoefs = mutableMapOf, C>() +// +// for (entry in coefs) { +// val key = entry.key.cleanUp() +// val value = entry.value +// fixedCoefs[key] = if (key in fixedCoefs) fixedCoefs[key]!! + value else value +// } +// +// return NumberedPolynomial( +// fixedCoefs +// .filter { it.value.isNotZero() } +// ) +//} +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from end of received +// * lists, sums up proportional monomials, removes zero monomials, and if result is empty map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * @param toCheckInput If it's `true` cleaning of [coefficients] is executed otherwise it is not. +// * +// * @throws PolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(NumberedPolynomialSpace) +//@Suppress("FunctionName") +//internal fun > NumberedPolynomial(pairs: Collection, C>>, toCheckInput: Boolean): NumberedPolynomial { +// if (!toCheckInput) return NumberedPolynomial(pairs.toMap()) +// +// val fixedCoefs = mutableMapOf, C>() +// +// for (entry in pairs) { +// val key = entry.first.cleanUp() +// val value = entry.second +// fixedCoefs[key] = if (key in fixedCoefs) fixedCoefs[key]!! + value else value +// } +// +// return NumberedPolynomial( +// fixedCoefs +// .filter { it.value.isNotZero() } +// ) +//} +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from end of received +// * lists, sums up proportional monomials, removes zero monomials, and if result is empty map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * @param toCheckInput If it's `true` cleaning of [coefficients] is executed otherwise it is not. +// * +// * @throws PolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(NumberedPolynomialSpace) +//@Suppress("FunctionName") +//internal fun > NumberedPolynomial(vararg pairs: Pair, C>, toCheckInput: Boolean): NumberedPolynomial = +// NumberedPolynomial(pairs.toList(), toCheckInput) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from end of received +// * lists, sums up proportional monomials, removes zero monomials, and if result is empty map adds only element in it. +// * +// * @param coefs Coefficients of the instants. +// * +// * @throws PolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial(coefs: Map, C>) = NumberedPolynomial(coefs, toCheckInput = true) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from end of received +// * lists, sums up proportional monomials, removes zero monomials, and if result is empty map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * +// * @throws PolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial(pairs: Collection, C>>) = NumberedPolynomial(pairs, toCheckInput = true) +///** +// * Gets the coefficients in format of [coefficients] field and cleans it: removes zero degrees from end of received +// * lists, sums up proportional monomials, removes zero monomials, and if result is empty map adds only element in it. +// * +// * @param pairs Collection of pairs that represent monomials. +// * +// * @throws PolynomialError If no coefficient received or if any of degrees in any monomial is negative. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial(vararg pairs: Pair, C>) = NumberedPolynomial(*pairs, toCheckInput = true) + +public fun > C.asNumberedPolynomial() : NumberedPolynomial = NumberedPolynomial(mapOf(emptyList() to this)) + +// endregion + +/** + * Space of polynomials. + * + * @param C the type of operated polynomials. + * @param A the intersection of [Ring] of [C] and [ScaleOperations] of [C]. + * @param ring the [A] instance. + */ +@Suppress("PARAMETER_NAME_CHANGED_ON_OVERRIDE", "INAPPLICABLE_JVM_NAME") +public class NumberedPolynomialSpace>( + 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 + + // 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 NumberedPolynomial.plus(other: Int): NumberedPolynomial = + if (other == 0) this + else + NumberedPolynomial( + coefficients + .toMutableMap() + .apply { + val degs = emptyList() + + val result = getOrElse(degs) { ring.zero } + other + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + public override operator fun NumberedPolynomial.minus(other: Int): NumberedPolynomial = + if (other == 0) this + else + NumberedPolynomial( + coefficients + .toMutableMap() + .apply { + val degs = emptyList() + + val result = getOrElse(degs) { ring.zero } - other + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + public override operator fun NumberedPolynomial.times(other: Int): NumberedPolynomial = + if (other == 0) zero + else NumberedPolynomial( + coefficients + .applyAndRemoveZeros { + mapValues { (_, c) -> c * other } + } + ) + // endregion + + // region Integer-polynomial relation + public override operator fun Int.plus(other: NumberedPolynomial): NumberedPolynomial = + if (this == 0) other + else + NumberedPolynomial( + other.coefficients + .toMutableMap() + .apply { + val degs = emptyList() + + val result = this@plus + getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + public override operator fun Int.minus(other: NumberedPolynomial): NumberedPolynomial = + if (this == 0) other + else + NumberedPolynomial( + other.coefficients + .toMutableMap() + .apply { + val degs = emptyList() + + val result = this@minus - getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + public override operator fun Int.times(other: NumberedPolynomial): NumberedPolynomial = + if (this == 0) zero + else NumberedPolynomial( + other.coefficients + .applyAndRemoveZeros { + mapValues { (_, c) -> this@times * c } + } + ) + // endregion + + // region Constant-constant relation + @JvmName("constantUnaryMinus") + override operator fun C.unaryMinus(): C = ring { -this@unaryMinus } + @JvmName("constantPlus") + override operator fun C.plus(other: C): C = ring { this@plus + other } + @JvmName("constantMinus") + override operator fun C.minus(other: C): C = ring { this@minus - other } + @JvmName("constantTimes") + 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 + override operator fun C.plus(other: NumberedPolynomial): NumberedPolynomial = + if (this.isZero()) other + else with(other.coefficients) { + if (isEmpty()) NumberedPolynomial(mapOf(emptyList() to this@plus)) + else NumberedPolynomial( + toMutableMap() + .apply { + val degs = emptyList() + + val result = this@plus + getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + override operator fun C.minus(other: NumberedPolynomial): NumberedPolynomial = + if (this.isZero()) -other + else with(other.coefficients) { + if (isEmpty()) NumberedPolynomial(mapOf(listOf() to this@minus)) + else NumberedPolynomial( + toMutableMap() + .apply { + forEach { (degs, c) -> if(degs.isNotEmpty()) this[degs] = -c } + + val degs = emptyList() + + val result = this@minus - getOrElse(degs) { ring.zero } + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + override operator fun C.times(other: NumberedPolynomial): NumberedPolynomial = + if (this.isZero()) zero + else NumberedPolynomial( + other.coefficients + .applyAndRemoveZeros { + mapValues { (_, c) -> this@times * c } + } + ) + // endregion + + // region Polynomial-constant relation + /** + * Returns sum of the polynomials. [other] is interpreted as [NumberedPolynomial]. + */ + override operator fun NumberedPolynomial.plus(other: C): NumberedPolynomial = + if (other.isZero()) this + else with(coefficients) { + if (isEmpty()) NumberedPolynomial(mapOf(listOf() to other)) + else NumberedPolynomial( + toMutableMap() + .apply { + val degs = emptyList() + + val result = getOrElse(degs) { ring.zero } + other + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + /** + * Returns difference of the polynomials. [other] is interpreted as [NumberedPolynomial]. + */ + override operator fun NumberedPolynomial.minus(other: C): NumberedPolynomial = + if (other.isZero()) this + else with(coefficients) { + if (isEmpty()) NumberedPolynomial(mapOf(listOf() to other)) + else NumberedPolynomial( + toMutableMap() + .apply { + val degs = emptyList() + + val result = getOrElse(degs) { ring.zero } - other + + if (result.isZero()) remove(degs) + else this[degs] = result + } + ) + } + /** + * Returns product of the polynomials. [other] is interpreted as [NumberedPolynomial]. + */ + override operator fun NumberedPolynomial.times(other: C): NumberedPolynomial = + if (other.isZero()) zero + else NumberedPolynomial( + coefficients + .applyAndRemoveZeros { + mapValues { (_, c) -> c * other } + } + ) + // endregion + + // region Polynomial-polynomial relation + /** + * Returns negation of the polynomial. + */ + override fun NumberedPolynomial.unaryMinus(): NumberedPolynomial = + NumberedPolynomial( + coefficients.mapValues { -it.value } + ) + /** + * Returns sum of the polynomials. + */ + override operator fun NumberedPolynomial.plus(other: NumberedPolynomial): NumberedPolynomial = + NumberedPolynomial( + coefficients + .applyAndRemoveZeros { + other.coefficients + .mapValuesTo(this) { (key, value) -> if (key in this) this[key]!! + value else value } + } + ) + /** + * Returns difference of the polynomials. + */ + override operator fun NumberedPolynomial.minus(other: NumberedPolynomial): NumberedPolynomial = + NumberedPolynomial( + coefficients + .applyAndRemoveZeros { + other.coefficients + .mapValuesTo(this) { (key, value) -> if (key in this) this[key]!! - value else -value } + } + ) + /** + * Returns product of the polynomials. + */ + override operator fun NumberedPolynomial.times(other: NumberedPolynomial): NumberedPolynomial = + when { + isZero() -> zero + other.isZero() -> zero + else -> + NumberedPolynomial( + buildCoefficients { + for ((degs1, c1) in coefficients) for ((degs2, c2) in other.coefficients) { + val degs = + (0..max(degs1.lastIndex, degs2.lastIndex)) + .map { degs1.getOrElse(it) { 0U } + degs2.getOrElse(it) { 0U } } + val c = c1 * c2 + this[degs] = if (degs in this) this[degs]!! + c else c + } + } + ) + } + + public override fun NumberedPolynomial.isZero(): Boolean = coefficients.values.all { it.isZero() } + public override fun NumberedPolynomial.isNotZero(): Boolean = coefficients.values.any { it.isNotZero() } + public override fun NumberedPolynomial.isOne(): Boolean = + with(coefficients) { + var foundAbsoluteTermAndItIsOne = false + for ((degs, c) in this) { + if (degs.isNotEmpty()) if (c.isNotZero()) return@with false + else { + if (c.isNotOne()) return@with false + else foundAbsoluteTermAndItIsOne = true + } + } + foundAbsoluteTermAndItIsOne + } + public override fun NumberedPolynomial.isNotOne(): Boolean = !isOne() + public override fun NumberedPolynomial.isMinusOne(): Boolean = + with(coefficients) { + var foundAbsoluteTermAndItIsMinusOne = false + for ((degs, c) in this) { + if (degs.isNotEmpty()) if (c.isNotZero()) return@with false + else { + if (c.isNotMinusOne()) return@with false + else foundAbsoluteTermAndItIsMinusOne = true + } + } + foundAbsoluteTermAndItIsMinusOne + } + public override fun NumberedPolynomial.isNotMinusOne(): Boolean = !isMinusOne() + + override val zero: NumberedPolynomial = NumberedPolynomial(emptyMap()) + override val one: NumberedPolynomial = + NumberedPolynomial( + mapOf( + listOf() to ring.one // 1 * x_1^0 * x_2^0 * ... + ) + ) + + // TODO: Docs + @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "CovariantEquals") + override fun NumberedPolynomial.equals(other: NumberedPolynomial): Boolean = + when { + this === other -> true + else -> coefficients.size == other.coefficients.size && + coefficients.all { (key, value) -> with(other.coefficients) { key in this && this[key] == value } } + } + // endregion + + // Not sure is it necessary... + // region Polynomial properties + /** + * Count of all variables that appear in the polynomial in positive exponents. + */ + public val NumberedPolynomial.countOfVariables: Int + get() = coefficients.entries.maxOfOrNull { (degs, c) -> if (c.isZero()) 0 else degs.size } ?: 0 + /** + * Degree of the polynomial, [see also](https://en.wikipedia.org/wiki/Degree_of_a_polynomial). If the polynomial is + * zero, degree is -1. + */ + override val NumberedPolynomial.degree: Int + get() = coefficients.entries.maxOfOrNull { (degs, c) -> if (c.isZero()) -1 else degs.sum().toInt() } ?: -1 + /** + * List that associates indices of variables (that appear in the polynomial in positive exponents) with their most + * exponents in which the variables are appeared in the polynomial. + * + * As consequence all values in the list are non-negative integers. Also, if the polynomial is constant, the list is empty. + * And size of the list is [countOfVariables]. + */ + public val NumberedPolynomial.degrees: List + get() = + buildList(countOfVariables) { + repeat(countOfVariables) { add(0U) } + coefficients.entries.forEach { (degs, c) -> + if (c.isNotZero()) degs.forEachIndexed { index, deg -> + this[index] = max(this[index], deg) + } + } + } + + /** + * Checks if the instant is constant polynomial (of degree no more than 0) over considered ring. + */ + override fun NumberedPolynomial.isConstant(): Boolean = + coefficients.all { (degs, c) -> degs.isEmpty() || c.isZero() } + /** + * Checks if the instant is constant non-zero polynomial (of degree no more than 0) over considered ring. + */ + override fun NumberedPolynomial.isNonZeroConstant(): Boolean = + with(coefficients) { + var foundAbsoluteTermAndItIsNotZero = false + for ((degs, c) in this) { + if (degs.isNotEmpty()) if (c.isNotZero()) return@with false + else { + if (c.isZero()) return@with false + else foundAbsoluteTermAndItIsNotZero = true + } + } + foundAbsoluteTermAndItIsNotZero + } + + override fun NumberedPolynomial.asConstantOrNull(): C? = + with(coefficients) { + if(isConstant()) getOrElse(emptyList()) { ring.zero } + else null + } + + @Suppress("NOTHING_TO_INLINE") + public inline fun NumberedPolynomial.substitute(argument: Map): NumberedPolynomial = this.substitute(ring, argument) + @Suppress("NOTHING_TO_INLINE") + @JvmName("substitutePolynomial") + public inline fun NumberedPolynomial.substitute(argument: Map>): NumberedPolynomial = this.substitute(ring, argument) + + @Suppress("NOTHING_TO_INLINE") + public inline fun NumberedPolynomial.asFunction(): (Map) -> NumberedPolynomial = { this.substitute(ring, it) } + @Suppress("NOTHING_TO_INLINE") + public inline fun NumberedPolynomial.asFunctionOnConstants(): (Map) -> NumberedPolynomial = { this.substitute(ring, it) } + @Suppress("NOTHING_TO_INLINE") + public inline fun NumberedPolynomial.asFunctionOnPolynomials(): (Map>) -> NumberedPolynomial = { this.substitute(ring, it) } + + @Suppress("NOTHING_TO_INLINE") + public inline operator fun NumberedPolynomial.invoke(argument: Map): NumberedPolynomial = this.substitute(ring, argument) + @Suppress("NOTHING_TO_INLINE") + @JvmName("invokePolynomial") + public inline operator fun NumberedPolynomial.invoke(argument: Map>): NumberedPolynomial = this.substitute(ring, argument) + // endregion + + // region Legacy + @Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE") + override inline fun add(left: NumberedPolynomial, right: NumberedPolynomial): NumberedPolynomial = left + right + @Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE") + override inline fun multiply(left: NumberedPolynomial, right: NumberedPolynomial): NumberedPolynomial = left * right + // endregion + + // region Utilities + // TODO: Move to region internal utilities with context receiver + internal fun MutableMap, C>.applyAndRemoveZeros(block: MutableMap, C>.() -> Unit) : MutableMap, C> { + contract { + callsInPlace(block, InvocationKind.EXACTLY_ONCE) + } + block() + for ((degs, c) in this) if (c.isZero()) this.remove(degs) + return this + } + internal fun Map, C>.applyAndRemoveZeros(block: MutableMap, C>.() -> Unit) : Map, C> = + toMutableMap().applyAndRemoveZeros(block) + @OptIn(ExperimentalTypeInference::class) + internal fun buildCoefficients(@BuilderInference builderAction: MutableMap, C>.() -> Unit): Map, C> { + contract { callsInPlace(builderAction, InvocationKind.EXACTLY_ONCE) } + return buildMap { + builderAction() + for ((degs, c) in this) if (c.isZero()) this.remove(degs) + } + } + // endregion +} \ No newline at end of file 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 30280a396..99d6b0659 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,43 +5,38 @@ package space.kscience.kmath.functions -import space.kscience.kmath.functions.AbstractPolynomialSpace.Companion.optimizedAddMultiplied -import space.kscience.kmath.functions.AbstractPolynomialSpace.Companion.optimizedMultiply import space.kscience.kmath.operations.* import kotlin.jvm.JvmName import kotlin.math.max import kotlin.math.min /** - * Polynomial coefficients model without fixation on specific context they are applied to. + * Polynomial model without fixation on specific context they are applied to. * * @param coefficients constant is the leftmost coefficient. */ 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 -// ) -// } } +// region Internal utilities + /** * Represents internal [Polynomial] errors. */ -internal class PolynomialError(message: String): Error(message) +internal class PolynomialError : Error { + constructor(): super() + constructor(message: String): super(message) + constructor(message: String?, cause: Throwable?): super(message, cause) + constructor(cause: Throwable?): super(cause) +} + +/** + * Throws an [PolynomialError] with the given [message]. + */ +internal fun polynomialError(message: Any): Nothing = throw PolynomialError(message.toString()) + +// endregion // region Constructors and converters @@ -66,16 +61,16 @@ public fun T.asPolynomial() : Polynomial = Polynomial(listOf(this)) // endregion /** - * Space of polynomials constructed over ring. + * Space of univariate 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]. */ -@Suppress("INAPPLICABLE_JVM_NAME") // KT-31420 +@Suppress("INAPPLICABLE_JVM_NAME") // TODO: KT-31420 public open class PolynomialSpace>( public val ring: A, -) : AbstractPolynomialSpace>{ +) : AbstractPolynomialSpace> { // region Constant-integer relation @JvmName("constantIntPlus") public override operator fun C.plus(other: Int): C = ring { optimizedAddMultiplied(this@plus, one, other) } @@ -102,8 +97,14 @@ public open class PolynomialSpace>( coefficients .toMutableList() .apply { - if (isEmpty()) this[0] = ring.zero + other - else this[0] = this[0]!! + other + val result = getOrElse(0) { ring.zero } + other + val isResultZero = result.isZero() + + when { + size == 0 && !isResultZero -> add(result) + size > 1 || !isResultZero -> this[0] = result + else -> clear() + } } ) public override operator fun Polynomial.minus(other: Int): Polynomial = @@ -113,8 +114,14 @@ public open class PolynomialSpace>( coefficients .toMutableList() .apply { - if (isEmpty()) this[0] = ring.zero - other - else this[0] = this[0]!! - other + val result = getOrElse(0) { ring.zero } - other + val isResultZero = result.isZero() + + when { + size == 0 && !isResultZero -> add(result) + size > 1 || !isResultZero -> this[0] = result + else -> clear() + } } ) public override operator fun Polynomial.times(other: Int): Polynomial = @@ -134,8 +141,14 @@ public open class PolynomialSpace>( other.coefficients .toMutableList() .apply { - if (isEmpty()) this[0] = ring.zero + this@plus - else this[0] = this[0]!! + this@plus + val result = this@plus + getOrElse(0) { ring.zero } + val isResultZero = result.isZero() + + when { + size == 0 && !isResultZero -> add(result) + size > 1 || !isResultZero -> this[0] = result + else -> clear() + } } ) public override operator fun Int.minus(other: Polynomial): Polynomial = @@ -145,8 +158,16 @@ public open class PolynomialSpace>( other.coefficients .toMutableList() .apply { - if (isEmpty()) this[0] = ring.zero - this@minus - else this[0] = this[0]!! - this@minus + forEachIndexed { index, c -> if (index != 0) this[index] = -c } + + val result = this@minus - getOrElse(0) { ring.zero } + val isResultZero = result.isZero() + + when { + size == 0 && !isResultZero -> add(result) + size > 1 || !isResultZero -> this[0] = result + else -> clear() + } } ) public override operator fun Int.times(other: Polynomial): Polynomial = @@ -183,12 +204,20 @@ public open class PolynomialSpace>( // region Constant-polynomial relation public override operator fun C.plus(other: Polynomial): Polynomial = - with(other.coefficients) { + if (this.isZero()) other + else with(other.coefficients) { if (isEmpty()) Polynomial(listOf(this@plus)) else Polynomial( toMutableList() .apply { - this[0] += this@plus + val result = if (size == 0) this@plus else this@plus + get(0) + val isResultZero = result.isZero() + + when { + size == 0 && !isResultZero -> add(result) + size > 1 || !isResultZero -> this[0] = result + else -> clear() + } } ) } @@ -196,12 +225,22 @@ public open class PolynomialSpace>( // listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1) // ) public override operator fun C.minus(other: Polynomial): Polynomial = - with(other.coefficients) { + if (this.isZero()) other + else with(other.coefficients) { if (isEmpty()) Polynomial(listOf(-this@minus)) else Polynomial( toMutableList() .apply { - this[0] -= this@minus + forEachIndexed { index, c -> if (index != 0) this[index] = -c } + + val result = if (size == 0) this@minus else this@minus - get(0) + val isResultZero = result.isZero() + + when { + size == 0 && !isResultZero -> add(result) + size > 1 || !isResultZero -> this[0] = result + else -> clear() + } } ) } @@ -209,9 +248,10 @@ public open class PolynomialSpace>( // listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1) // ) public override operator fun C.times(other: Polynomial): Polynomial = - Polynomial( + if (this.isZero()) other + else Polynomial( other.coefficients -// .subList(0, other.degree + 1) + .subList(0, other.degree + 1) .map { it * this } ) // endregion @@ -224,7 +264,14 @@ public open class PolynomialSpace>( else Polynomial( toMutableList() .apply { - this[0] += other + val result = if (size == 0) other else get(0) + other + val isResultZero = result.isZero() + + when { + size == 0 && !isResultZero -> add(result) + size > 1 || !isResultZero -> this[0] = result + else -> clear() + } } ) } @@ -232,12 +279,20 @@ public open class PolynomialSpace>( // listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1) // ) public override operator fun Polynomial.minus(other: C): Polynomial = - with(coefficients) { + if (other.isZero()) this + else with(coefficients) { if (isEmpty()) Polynomial(listOf(other)) else Polynomial( toMutableList() .apply { - this[0] -= other + val result = if (size == 0) other else get(0) - other + val isResultZero = result.isZero() + + when { + size == 0 && !isResultZero -> add(result) + size > 1 || !isResultZero -> this[0] = result + else -> clear() + } } ) } @@ -245,9 +300,10 @@ public open class PolynomialSpace>( // listOf(coefficients[0] - other) + coefficients.subList(1, degree + 1) // ) public override operator fun Polynomial.times(other: C): Polynomial = - Polynomial( + if (other.isZero()) this + else Polynomial( coefficients -// .subList(0, degree + 1) + .subList(0, degree + 1) .map { it * other } ) // endregion @@ -283,8 +339,8 @@ public open class PolynomialSpace>( val thisDegree = degree val otherDegree = other.degree return when { - thisDegree == -1 -> this - otherDegree == -1 -> other + thisDegree == -1 -> zero + otherDegree == -1 -> zero else -> Polynomial( (0..(thisDegree + otherDegree)) @@ -293,6 +349,7 @@ public open class PolynomialSpace>( .map { coefficients[it] * other.coefficients[d - it] } .reduce { acc, rational -> acc + rational } } + .run { subList(0, indexOfLast { it.isNotZero() } + 1) } ) } } @@ -321,8 +378,8 @@ public open class PolynomialSpace>( } // 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? = @@ -354,8 +411,8 @@ public open class PolynomialSpace>( 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 + // endregion } /** diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Variable.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Variable.kt new file mode 100644 index 000000000..410604fd3 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Variable.kt @@ -0,0 +1,38 @@ +package space.kscience.kmath.functions + +import kotlin.reflect.KProperty + + +/** + * Represents class of labeled variables like usual + * `x`, `y`, `z`, `a`, `b`, `n`, `m`, etc. + * + * Variables does not contain any information about field (or ring, ets.) they are considered in + * and therefore about coefficient. + * + * @property name Is the label or name of variable. For `x` it is `"x"`, for `n` – `"n"`, etc. + */ +public data class Variable (val name: String) : Comparable { + /** + * Represents the variable as a string. + * + * @return Only name of the variable. + */ + override fun toString(): String = name + /** + * Compares two variables. + * Comparison is realised by comparison of variables' names. + * + * Used in [LabeledPolynomial] and [LabeledRationalFunction] to sort monomials in + * [LabeledPolynomial.toString] and [LabeledRationalFunction.toString] in lexicographic order. + * + * @see Comparable.compareTo + * @sample LabeledPolynomial.monomialComparator + * @return Only name of the variable. + */ + override fun compareTo(other: Variable): Int = name.compareTo(other.name) + + public companion object { + public operator fun getValue(thisRef: Any?, property: KProperty<*>) : Variable = Variable(property.name) + } +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/algebraicStub.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/algebraicStub.kt new file mode 100644 index 000000000..9e5043b8c --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/algebraicStub.kt @@ -0,0 +1,51 @@ +/* + * 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.Group + + +// 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/labeledPolynomialUtil.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledPolynomialUtil.kt new file mode 100644 index 000000000..62ac31b64 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledPolynomialUtil.kt @@ -0,0 +1,490 @@ +package space.kscience.kmath.functions + +import space.kscience.kmath.operations.* +import kotlin.contracts.* + + +// TODO: Docs + +// region Sort of legacy + +//// region Constants +// +//// TODO: Reuse underlying ring extensions +// +//context(LabeledPolynomialSpace) +//@Suppress("NOTHING_TO_INLINE") +//fun > numberConstant(value: Int): C = ring { number(value) } +// +//context(LabeledPolynomialSpace) +//fun > power(arg: C, pow: UInt): C = ring { power(arg, pow) } +// +//context(LabeledPolynomialSpace) +//fun > multiplyWithPower(base: C, arg: C, pow: UInt): C = ring { multiplyWithPower(base, arg, pow) } +// +//// endregion + +//// region Variables +// +//context(LabeledPolynomialSpace) +//fun > power(arg: Variable, pow: UInt): LabeledPolynomial = +// if (pow == 0U) one +// else LabeledPolynomial(mapOf( +// mapOf(arg to pow) to ring.one +// )) +// +//// endregion + +//// region Polynomials +// +//context(LabeledPolynomialSpace) +//fun > number(value: Int): LabeledPolynomial = ring { LabeledPolynomial(mapOf(emptyMap() to number(value))) } +// +//context(LabeledPolynomialSpace) +//fun > multiplyWithPower(base: LabeledPolynomial, arg: LabeledPolynomial, pow: UInt): LabeledPolynomial = +// when { +// arg.isZero() && pow > 0U -> base +// arg.isOne() -> base +// arg.isMinusOne() -> if (pow % 2U == 0U) base else -base +// else -> multiplyWithPowerInternalLogic(base, arg, pow) +// } +// +//// Trivial but slow as duck +//context(LabeledPolynomialSpace) +//internal tailrec fun > multiplyWithPowerInternalLogic(base: LabeledPolynomial, arg: LabeledPolynomial, exponent: UInt): LabeledPolynomial = +// when { +// exponent == 0U -> base +// exponent == 1U -> base * arg +// exponent % 2U == 0U -> multiplyWithPowerInternalLogic(base, arg * arg, exponent / 2U) +// exponent % 2U == 1U -> multiplyWithPowerInternalLogic(base * arg, arg * arg, exponent / 2U) +// else -> error("Error in raising ring instant by unsigned integer: got reminder by division by 2 different from 0 and 1") +// } +// +//// endregion + +// endregion + +// region Utilities + +// TODO: Docs +@OptIn(ExperimentalContracts::class) +public inline fun , R> A.labeledPolynomial(block: LabeledPolynomialSpace.() -> R): R { + contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) } + return LabeledPolynomialSpace(this).block() +} + +// endregion + +//// region String representations +// +///** +// * Represents the polynomial as a [String] with names of variables substituted with names from [names]. +// * Consider that monomials are sorted in lexicographic order. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial.represent(names: Map = emptyMap()): String = +// coefficients.entries +// .sortedWith { o1, o2 -> LabeledPolynomial.monomialComparator.compare(o1.key, o2.key) } +// .asSequence() +// .map { (degs, t) -> +// if (degs.isEmpty()) "$t" +// else { +// when { +// t.isOne() -> "" +// t.isMinusOne() -> "-" +// else -> "$t " +// } + +// degs +// .toSortedMap() +// .filter { it.value > 0U } +// .map { (variable, deg) -> +// val variableName = names.getOrDefault(variable, variable.toString()) +// when (deg) { +// 1U -> variableName +// else -> "$variableName^$deg" +// } +// } +// .joinToString(separator = " ") { it } +// } +// } +// .joinToString(separator = " + ") { it } +// .ifEmpty { "0" } +// +///** +// * Represents the polynomial as a [String] naming variables by [namer]. +// * Consider that monomials are sorted in lexicographic order. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial.represent(namer: (Variable) -> String): String = +// coefficients.entries +// .sortedWith { o1, o2 -> LabeledPolynomial.monomialComparator.compare(o1.key, o2.key) } +// .asSequence() +// .map { (degs, t) -> +// if (degs.isEmpty()) "$t" +// else { +// when { +// t.isOne() -> "" +// t.isMinusOne() -> "-" +// else -> "$t " +// } + +// degs +// .toSortedMap() +// .filter { it.value > 0U } +// .map { (variable, deg) -> +// when (deg) { +// 1U -> namer(variable) +// else -> "${namer(variable)}^$deg" +// } +// } +// .joinToString(separator = " ") { it } +// } +// } +// .joinToString(separator = " + ") { it } +// .ifEmpty { "0" } +// +///** +// * Represents the polynomial as a [String] with names of variables substituted with names from [names] and with +// * brackets around the string if needed (i.e. when there are at least two addends in the representation). +// * Consider that monomials are sorted in lexicographic order. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial.representWithBrackets(names: Map = emptyMap()): String = +// with(represent(names)) { if (coefficients.count() == 1) this else "($this)" } +// +///** +// * Represents the polynomial as a [String] naming variables by [namer] and with brackets around the string if needed +// * (i.e. when there are at least two addends in the representation). +// * Consider that monomials are sorted in lexicographic order. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial.representWithBrackets(namer: (Variable) -> String): String = +// with(represent(namer)) { if (coefficients.count() == 1) this else "($this)" } +// +///** +// * Represents the polynomial as a [String] with names of variables substituted with names from [names]. +// * Consider that monomials are sorted in **reversed** lexicographic order. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial.representReversed(names: Map = emptyMap()): String = +// coefficients.entries +// .sortedWith { o1, o2 -> -LabeledPolynomial.monomialComparator.compare(o1.key, o2.key) } +// .asSequence() +// .map { (degs, t) -> +// if (degs.isEmpty()) "$t" +// else { +// when { +// t.isOne() -> "" +// t.isMinusOne() -> "-" +// else -> "$t " +// } + +// degs +// .toSortedMap() +// .filter { it.value > 0U } +// .map { (variable, deg) -> +// val variableName = names.getOrDefault(variable, variable.toString()) +// when (deg) { +// 1U -> variableName +// else -> "$variableName^$deg" +// } +// } +// .joinToString(separator = " ") { it } +// } +// } +// .joinToString(separator = " + ") { it } +// .ifEmpty { "0" } +// +///** +// * Represents the polynomial as a [String] naming variables by [namer]. +// * Consider that monomials are sorted in **reversed** lexicographic order. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial.representReversed(namer: (Variable) -> String): String = +// coefficients.entries +// .sortedWith { o1, o2 -> -LabeledPolynomial.monomialComparator.compare(o1.key, o2.key) } +// .asSequence() +// .map { (degs, t) -> +// if (degs.isEmpty()) "$t" +// else { +// when { +// t.isOne() -> "" +// t.isMinusOne() -> "-" +// else -> "$t " +// } + +// degs +// .toSortedMap() +// .filter { it.value > 0U } +// .map { (variable, deg) -> +// when (deg) { +// 1U -> namer(variable) +// else -> "${namer(variable)}^$deg" +// } +// } +// .joinToString(separator = " ") { it } +// } +// } +// .joinToString(separator = " + ") { it } +// .ifEmpty { "0" } +// +///** +// * Represents the polynomial as a [String] with names of variables substituted with names from [names] and with +// * brackets around the string if needed (i.e. when there are at least two addends in the representation). +// * Consider that monomials are sorted in **reversed** lexicographic order. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial.representReversedWithBrackets(names: Map = emptyMap()): String = +// with(representReversed(names)) { if (coefficients.count() == 1) this else "($this)" } +// +///** +// * Represents the polynomial as a [String] naming variables by [namer] and with brackets around the string if needed +// * (i.e. when there are at least two addends in the representation). +// * Consider that monomials are sorted in **reversed** lexicographic order. +// */ +//context(LabeledPolynomialSpace) +//fun > LabeledPolynomial.representReversedWithBrackets(namer: (Variable) -> String): String = +// with(representReversed(namer)) { if (coefficients.count() == 1) this else "($this)" } +// +//// endregion + +// region Operator extensions + +//// region Field case +// +//operator fun > Polynomial.div(other: Polynomial): Polynomial { +// if (other.isZero()) throw ArithmeticException("/ by zero") +// if (isZero()) return this +// +// fun Map, T>.leadingTerm() = +// this +// .asSequence() +// .map { Pair(it.key, it.value) } +// .reduce { (accDegs, accC), (listDegs, listC) -> +// for (i in 0..accDegs.lastIndex) { +// if (accDegs[i] > listDegs.getOrElse(i) { 0 }) return@reduce accDegs to accC +// if (accDegs[i] < listDegs.getOrElse(i) { 0 }) return@reduce listDegs to listC +// } +// if (accDegs.size < listDegs.size) listDegs to listC else accDegs to accC +// } +// +// var thisCoefficients = coefficients.toMutableMap() +// val otherCoefficients = other.coefficients +// val quotientCoefficients = HashMap, T>() +// +// var (thisLeadingTermDegs, thisLeadingTermC) = thisCoefficients.leadingTerm() +// val (otherLeadingTermDegs, otherLeadingTermC) = otherCoefficients.leadingTerm() +// +// while ( +// thisLeadingTermDegs.size >= otherLeadingTermDegs.size && +// (0..otherLeadingTermDegs.lastIndex).all { thisLeadingTermDegs[it] >= otherLeadingTermDegs[it] } +// ) { +// val multiplierDegs = +// thisLeadingTermDegs +// .mapIndexed { index, deg -> deg - otherLeadingTermDegs.getOrElse(index) { 0 } } +// .cleanUp() +// val multiplierC = thisLeadingTermC / otherLeadingTermC +// +// quotientCoefficients[multiplierDegs] = multiplierC +// +// for ((degs, t) in otherCoefficients) { +// val productDegs = +// (0..max(degs.lastIndex, multiplierDegs.lastIndex)) +// .map { degs.getOrElse(it) { 0 } + multiplierDegs.getOrElse(it) { 0 } } +// .cleanUp() +// val productC = t * multiplierC +// thisCoefficients[productDegs] = +// if (productDegs in thisCoefficients) thisCoefficients[productDegs]!! - productC else -productC +// } +// +// thisCoefficients = thisCoefficients.filterValues { it.isNotZero() }.toMutableMap() +// +// if (thisCoefficients.isEmpty()) +// return Polynomial(quotientCoefficients, toCheckInput = false) +// +// val t = thisCoefficients.leadingTerm() +// thisLeadingTermDegs = t.first +// thisLeadingTermC = t.second +// } +// +// return Polynomial(quotientCoefficients, toCheckInput = false) +//} +// +//operator fun > Polynomial.div(other: T): Polynomial = +// if (other.isZero()) throw ArithmeticException("/ by zero") +// else +// Polynomial( +// coefficients +// .mapValues { it.value / other }, +// toCheckInput = false +// ) +// +//operator fun > Polynomial.rem(other: Polynomial): Polynomial { +// if (other.isZero()) throw ArithmeticException("/ by zero") +// if (isZero()) return this +// +// fun Map, T>.leadingTerm() = +// this +// .asSequence() +// .map { Pair(it.key, it.value) } +// .reduce { (accDegs, accC), (listDegs, listC) -> +// for (i in 0..accDegs.lastIndex) { +// if (accDegs[i] > listDegs.getOrElse(i) { 0 }) return@reduce accDegs to accC +// if (accDegs[i] < listDegs.getOrElse(i) { 0 }) return@reduce listDegs to listC +// } +// if (accDegs.size < listDegs.size) listDegs to listC else accDegs to accC +// } +// +// var thisCoefficients = coefficients.toMutableMap() +// val otherCoefficients = other.coefficients +// +// var (thisLeadingTermDegs, thisLeadingTermC) = thisCoefficients.leadingTerm() +// val (otherLeadingTermDegs, otherLeadingTermC) = otherCoefficients.leadingTerm() +// +// while ( +// thisLeadingTermDegs.size >= otherLeadingTermDegs.size && +// (0..otherLeadingTermDegs.lastIndex).all { thisLeadingTermDegs[it] >= otherLeadingTermDegs[it] } +// ) { +// val multiplierDegs = +// thisLeadingTermDegs +// .mapIndexed { index, deg -> deg - otherLeadingTermDegs.getOrElse(index) { 0 } } +// .cleanUp() +// val multiplierC = thisLeadingTermC / otherLeadingTermC +// +// for ((degs, t) in otherCoefficients) { +// val productDegs = +// (0..max(degs.lastIndex, multiplierDegs.lastIndex)) +// .map { degs.getOrElse(it) { 0 } + multiplierDegs.getOrElse(it) { 0 } } +// .cleanUp() +// val productC = t * multiplierC +// thisCoefficients[productDegs] = +// if (productDegs in thisCoefficients) thisCoefficients[productDegs]!! - productC else -productC +// } +// +// thisCoefficients = thisCoefficients.filterValues { it.isNotZero() }.toMutableMap() +// +// if (thisCoefficients.isEmpty()) +// return Polynomial(thisCoefficients, toCheckInput = false) +// +// val t = thisCoefficients.leadingTerm() +// thisLeadingTermDegs = t.first +// thisLeadingTermC = t.second +// } +// +// return Polynomial(thisCoefficients, toCheckInput = false) +//} +// +//infix fun > Polynomial.divrem(other: Polynomial): Polynomial.Companion.DividingResult { +// if (other.isZero()) throw ArithmeticException("/ by zero") +// if (isZero()) return Polynomial.Companion.DividingResult(this, this) +// +// fun Map, T>.leadingTerm() = +// this +// .asSequence() +// .map { Pair(it.key, it.value) } +// .reduce { (accDegs, accC), (listDegs, listC) -> +// for (i in 0..accDegs.lastIndex) { +// if (accDegs[i] > listDegs.getOrElse(i) { 0 }) return@reduce accDegs to accC +// if (accDegs[i] < listDegs.getOrElse(i) { 0 }) return@reduce listDegs to listC +// } +// if (accDegs.size < listDegs.size) listDegs to listC else accDegs to accC +// } +// +// var thisCoefficients = coefficients.toMutableMap() +// val otherCoefficients = other.coefficients +// val quotientCoefficients = HashMap, T>() +// +// var (thisLeadingTermDegs, thisLeadingTermC) = thisCoefficients.leadingTerm() +// val (otherLeadingTermDegs, otherLeadingTermC) = otherCoefficients.leadingTerm() +// +// while ( +// thisLeadingTermDegs.size >= otherLeadingTermDegs.size && +// (0..otherLeadingTermDegs.lastIndex).all { thisLeadingTermDegs[it] >= otherLeadingTermDegs[it] } +// ) { +// val multiplierDegs = +// thisLeadingTermDegs +// .mapIndexed { index, deg -> deg - otherLeadingTermDegs.getOrElse(index) { 0 } } +// .cleanUp() +// val multiplierC = thisLeadingTermC / otherLeadingTermC +// +// quotientCoefficients[multiplierDegs] = multiplierC +// +// for ((degs, t) in otherCoefficients) { +// val productDegs = +// (0..max(degs.lastIndex, multiplierDegs.lastIndex)) +// .map { degs.getOrElse(it) { 0 } + multiplierDegs.getOrElse(it) { 0 } } +// .cleanUp() +// val productC = t * multiplierC +// thisCoefficients[productDegs] = +// if (productDegs in thisCoefficients) thisCoefficients[productDegs]!! - productC else -productC +// } +// +// thisCoefficients = thisCoefficients.filterValues { it.isNotZero() }.toMutableMap() +// +// if (thisCoefficients.isEmpty()) +// return Polynomial.Companion.DividingResult( +// Polynomial(quotientCoefficients, toCheckInput = false), +// Polynomial(thisCoefficients, toCheckInput = false) +// ) +// +// val t = thisCoefficients.leadingTerm() +// thisLeadingTermDegs = t.first +// thisLeadingTermC = t.second +// } +// +// return Polynomial.Companion.DividingResult( +// Polynomial(quotientCoefficients, toCheckInput = false), +// Polynomial(thisCoefficients, toCheckInput = false) +// ) +//} +// +//// endregion + +// endregion + +//// region Polynomial substitution and functional representation +// +//public fun LabeledPolynomial.substitute(ring: Ring, args: Map): LabeledPolynomial = ring { +// if (coefficients.isEmpty()) return this@substitute +// LabeledPolynomial( +// buildMap { +// coefficients.forEach { (degs, c) -> +// val newDegs = degs.filterKeys { it !in args } +// val newC = degs.entries.asSequence().filter { it.key in args }.fold(c) { acc, (variable, deg) -> +// multiplyWithPower(acc, args[variable]!!, deg) +// } +// this[newDegs] = if (newDegs in this) this[newDegs]!! + newC else newC +// } +// } +// ) +//} +// +//// TODO: Replace with optimisation: the [result] may be unboxed, and all operations may be performed as soon as +//// possible on it +//@JvmName("substitutePolynomial") +//fun LabeledPolynomial.substitute(ring: Ring, arg: Map>) : LabeledPolynomial = +// ring.labeledPolynomial { +// if (coefficients.isEmpty()) return zero +// coefficients +// .asSequence() +// .map { (degs, c) -> +// degs.entries +// .asSequence() +// .filter { it.key in arg } +// .fold(LabeledPolynomial(mapOf(degs.filterKeys { it !in arg } to c))) { acc, (index, deg) -> +// multiplyWithPower(acc, arg[index]!!, deg) +// } +// } +// .reduce { acc, polynomial -> acc + polynomial } // TODO: Rewrite. Might be slow. +// } +// +//// TODO: Substitute rational function +// +//fun > LabeledPolynomial.asFunctionOver(ring: A): (Map) -> LabeledPolynomial = +// { substitute(ring, it) } +// +//fun > LabeledPolynomial.asPolynomialFunctionOver(ring: A): (Map>) -> LabeledPolynomial = +// { substitute(ring, it) } +// +//// endregion + +//// region Algebraic derivative and antiderivative +//// TODO +//// endregion \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedPolynomialUtil.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedPolynomialUtil.kt new file mode 100644 index 000000000..d4053442d --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedPolynomialUtil.kt @@ -0,0 +1,605 @@ +package space.kscience.kmath.functions + +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.operations.* +import kotlin.contracts.* +import kotlin.jvm.JvmName + + +// TODO: Docs + +// region Sort of legacy + +//// region Constants +// +//// TODO: Reuse underlying ring extensions +// +//context(NumberedPolynomialSpace) +//@Suppress("NOTHING_TO_INLINE") +//public fun > numberConstant(value: Int): C = ring { number(value) } +// +//context(NumberedPolynomialSpace) +//public fun > multiplyWithPower(base: C, arg: C, pow: UInt): C = ring { multiplyWithPower(base, arg, pow) } +// +//// endregion + +//// region Polynomials +// +//context(NumberedPolynomialSpace) +//public fun > number(value: Int): NumberedPolynomial = ring { NumberedPolynomial(mapOf(emptyList() to number(value))) } +// +//context(NumberedPolynomialSpace) +//public fun > multiplyWithPower(base: NumberedPolynomial, arg: NumberedPolynomial, pow: UInt): NumberedPolynomial = +// when { +// arg.isZero() && pow > 0U -> base +// arg.isOne() -> base +// arg.isMinusOne() -> if (pow % 2U == 0U) base else -base +// else -> multiplyWithPowerInternalLogic(base, arg, pow) +// } +// +//// Trivial but slow as duck +//context(NumberedPolynomialSpace) +//internal tailrec fun > multiplyWithPowerInternalLogic(base: NumberedPolynomial, arg: NumberedPolynomial, exponent: UInt): NumberedPolynomial = +// when { +// exponent == 0U -> base +// exponent == 1U -> base * arg +// exponent % 2U == 0U -> multiplyWithPowerInternalLogic(base, arg * arg, exponent / 2U) +// exponent % 2U == 1U -> multiplyWithPowerInternalLogic(base * arg, arg * arg, exponent / 2U) +// else -> error("Error in raising ring instant by unsigned integer: got reminder by division by 2 different from 0 and 1") +// } +// +//// endregion + +// endregion + +// region Utilities + +// TODO: Docs +@OptIn(ExperimentalContracts::class) +public inline fun , R> A.numberedPolynomial(block: NumberedPolynomialSpace.() -> R): R { + contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) } + return NumberedPolynomialSpace(this).block() +} + +// endregion + +//// region String representations +// +///** +// * Represents the polynomial as a [String] where name of variable with index `i` is [withVariableName] + `"_${i+1}"`. +// * Consider that monomials are sorted in lexicographic order. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial.represent(withVariableName: String = NumberedPolynomial.defaultVariableName): String = +// coefficients.entries +// .sortedWith { o1, o2 -> NumberedPolynomial.monomialComparator.compare(o1.key, o2.key) } +// .asSequence() +// .map { (degs, t) -> +// if (degs.isEmpty()) "$t" +// else { +// when { +// t.isOne() -> "" +// t.isMinusOne() -> "-" +// else -> "$t " +// } + +// degs +// .mapIndexed { index, deg -> +// when (deg) { +// 0U -> "" +// 1U -> "${withVariableName}_${index+1}" +// else -> "${withVariableName}_${index+1}^$deg" +// } +// } +// .filter { it.isNotEmpty() } +// .joinToString(separator = " ") { it } +// } +// } +// .joinToString(separator = " + ") { it } +// .ifEmpty { "0" } +// +///** +// * Represents the polynomial as a [String] naming variables by [namer]. +// * Consider that monomials are sorted in lexicographic order. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial.represent(namer: (Int) -> String): String = +// coefficients.entries +// .sortedWith { o1, o2 -> NumberedPolynomial.monomialComparator.compare(o1.key, o2.key) } +// .asSequence() +// .map { (degs, t) -> +// if (degs.isEmpty()) "$t" +// else { +// when { +// t.isOne() -> "" +// t.isMinusOne() -> "-" +// else -> "$t " +// } + +// degs +// .mapIndexed { index, deg -> +// when (deg) { +// 0U -> "" +// 1U -> namer(index) +// else -> "${namer(index)}^$deg" +// } +// } +// .filter { it.isNotEmpty() } +// .joinToString(separator = " ") { it } +// } +// } +// .joinToString(separator = " + ") { it } +// .ifEmpty { "0" } +// +///** +// * Represents the polynomial as a [String] where name of variable with index `i` is [withVariableName] + `"_${i+1}"` +// * and with brackets around the string if needed (i.e. when there are at least two addends in the representation). +// * Consider that monomials are sorted in lexicographic order. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial.representWithBrackets(withVariableName: String = NumberedPolynomial.defaultVariableName): String = +// with(represent(withVariableName)) { if (coefficients.count() == 1) this else "($this)" } +// +///** +// * Represents the polynomial as a [String] naming variables by [namer] and with brackets around the string if needed +// * (i.e. when there are at least two addends in the representation). +// * Consider that monomials are sorted in lexicographic order. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial.representWithBrackets(namer: (Int) -> String): String = +// with(represent(namer)) { if (coefficients.count() == 1) this else "($this)" } +// +///** +// * Represents the polynomial as a [String] where name of variable with index `i` is [withVariableName] + `"_${i+1}"`. +// * Consider that monomials are sorted in **reversed** lexicographic order. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial.representReversed(withVariableName: String = NumberedPolynomial.defaultVariableName): String = +// coefficients.entries +// .sortedWith { o1, o2 -> -NumberedPolynomial.monomialComparator.compare(o1.key, o2.key) } +// .asSequence() +// .map { (degs, t) -> +// if (degs.isEmpty()) "$t" +// else { +// when { +// t.isOne() -> "" +// t.isMinusOne() -> "-" +// else -> "$t " +// } + +// degs +// .mapIndexed { index, deg -> +// when (deg) { +// 0U -> "" +// 1U -> "${withVariableName}_${index+1}" +// else -> "${withVariableName}_${index+1}^$deg" +// } +// } +// .filter { it.isNotEmpty() } +// .joinToString(separator = " ") { it } +// } +// } +// .joinToString(separator = " + ") { it } +// .ifEmpty { "0" } +// +///** +// * Represents the polynomial as a [String] naming variables by [namer]. +// * Consider that monomials are sorted in **reversed** lexicographic order. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial.representReversed(namer: (Int) -> String): String = +// coefficients.entries +// .sortedWith { o1, o2 -> -NumberedPolynomial.monomialComparator.compare(o1.key, o2.key) } +// .asSequence() +// .map { (degs, t) -> +// if (degs.isEmpty()) "$t" +// else { +// when { +// t.isOne() -> "" +// t.isMinusOne() -> "-" +// else -> "$t " +// } + +// degs +// .mapIndexed { index, deg -> +// when (deg) { +// 0U -> "" +// 1U -> namer(index) +// else -> "${namer(index)}^$deg" +// } +// } +// .filter { it.isNotEmpty() } +// .joinToString(separator = " ") { it } +// } +// } +// .joinToString(separator = " + ") { it } +// .ifEmpty { "0" } +// +///** +// * Represents the polynomial as a [String] where name of variable with index `i` is [withVariableName] + `"_${i+1}"` +// * and with brackets around the string if needed (i.e. when there are at least two addends in the representation). +// * Consider that monomials are sorted in **reversed** lexicographic order. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial.representReversedWithBrackets(withVariableName: String = NumberedPolynomial.defaultVariableName): String = +// with(representReversed(withVariableName)) { if (coefficients.count() == 1) this else "($this)" } +// +///** +// * Represents the polynomial as a [String] naming variables by [namer] and with brackets around the string if needed +// * (i.e. when there are at least two addends in the representation). +// * Consider that monomials are sorted in **reversed** lexicographic order. +// */ +//context(NumberedPolynomialSpace) +//public fun > NumberedPolynomial.representReversedWithBrackets(namer: (Int) -> String): String = +// with(representReversed(namer)) { if (coefficients.count() == 1) this else "($this)" } +// +//// endregion + +//// region Polynomial substitution and functional representation +// +//public fun NumberedPolynomial.substitute(ring: Ring, args: Map): NumberedPolynomial = ring { +// if (coefficients.isEmpty()) return this@substitute +// NumberedPolynomial( +// buildMap { +// coefficients.forEach { (degs, c) -> +// val newDegs = degs.mapIndexed { index, deg -> if (index in args) 0U else deg }.cleanUp() +// val newC = degs.foldIndexed(c) { index, acc, deg -> +// if (index in args) multiplyWithPower(acc, args[index]!!, deg) +// else acc +// } +// this[newDegs] = if (newDegs in this) this[newDegs]!! + newC else newC +// } +// } +// ) +//} +// +//// TODO: Replace with optimisation: the [result] may be unboxed, and all operations may be performed as soon as +//// possible on it +//@JvmName("substitutePolynomial") +//public fun NumberedPolynomial.substitute(ring: Ring, arg: Map>) : NumberedPolynomial = +// ring.numberedPolynomialSpace { +// if (coefficients.isEmpty()) return zero +// coefficients +// .asSequence() +// .map { (degs, c) -> +// degs.foldIndexed( +// NumberedPolynomial( +// degs.mapIndexed { index, deg -> if (index in arg) 0U else deg } to c +// ) +// ) { index, acc, deg -> if (index in arg) multiplyWithPower(acc, arg[index]!!, deg) else acc } +// } +// .reduce { acc, polynomial -> acc + polynomial } // TODO: Rewrite. Might be slow. +// } +// +//// TODO: Substitute rational function +// +//public fun > NumberedPolynomial.asFunctionOver(ring: A): (Map) -> NumberedPolynomial = +// { substitute(ring, it) } +// +//public fun > NumberedPolynomial.asPolynomialFunctionOver(ring: A): (Map>) -> NumberedPolynomial = +// { substitute(ring, it) } +// +//// endregion + +// region Operator extensions + +//// region Field case +// +//operator fun > Polynomial.div(other: Polynomial): Polynomial { +// if (other.isZero()) throw ArithmeticException("/ by zero") +// if (isZero()) return this +// +// fun Map, T>.leadingTerm() = +// this +// .asSequence() +// .map { Pair(it.key, it.value) } +// .reduce { (accDegs, accC), (listDegs, listC) -> +// for (i in 0..accDegs.lastIndex) { +// if (accDegs[i] > listDegs.getOrElse(i) { 0 }) return@reduce accDegs to accC +// if (accDegs[i] < listDegs.getOrElse(i) { 0 }) return@reduce listDegs to listC +// } +// if (accDegs.size < listDegs.size) listDegs to listC else accDegs to accC +// } +// +// var thisCoefficients = coefficients.toMutableMap() +// val otherCoefficients = other.coefficients +// val quotientCoefficients = HashMap, T>() +// +// var (thisLeadingTermDegs, thisLeadingTermC) = thisCoefficients.leadingTerm() +// val (otherLeadingTermDegs, otherLeadingTermC) = otherCoefficients.leadingTerm() +// +// while ( +// thisLeadingTermDegs.size >= otherLeadingTermDegs.size && +// (0..otherLeadingTermDegs.lastIndex).all { thisLeadingTermDegs[it] >= otherLeadingTermDegs[it] } +// ) { +// val multiplierDegs = +// thisLeadingTermDegs +// .mapIndexed { index, deg -> deg - otherLeadingTermDegs.getOrElse(index) { 0 } } +// .cleanUp() +// val multiplierC = thisLeadingTermC / otherLeadingTermC +// +// quotientCoefficients[multiplierDegs] = multiplierC +// +// for ((degs, t) in otherCoefficients) { +// val productDegs = +// (0..max(degs.lastIndex, multiplierDegs.lastIndex)) +// .map { degs.getOrElse(it) { 0 } + multiplierDegs.getOrElse(it) { 0 } } +// .cleanUp() +// val productC = t * multiplierC +// thisCoefficients[productDegs] = +// if (productDegs in thisCoefficients) thisCoefficients[productDegs]!! - productC else -productC +// } +// +// thisCoefficients = thisCoefficients.filterValues { it.isNotZero() }.toMutableMap() +// +// if (thisCoefficients.isEmpty()) +// return Polynomial(quotientCoefficients, toCheckInput = false) +// +// val t = thisCoefficients.leadingTerm() +// thisLeadingTermDegs = t.first +// thisLeadingTermC = t.second +// } +// +// return Polynomial(quotientCoefficients, toCheckInput = false) +//} +// +//operator fun > Polynomial.div(other: T): Polynomial = +// if (other.isZero()) throw ArithmeticException("/ by zero") +// else +// Polynomial( +// coefficients +// .mapValues { it.value / other }, +// toCheckInput = false +// ) +// +//operator fun > Polynomial.rem(other: Polynomial): Polynomial { +// if (other.isZero()) throw ArithmeticException("/ by zero") +// if (isZero()) return this +// +// fun Map, T>.leadingTerm() = +// this +// .asSequence() +// .map { Pair(it.key, it.value) } +// .reduce { (accDegs, accC), (listDegs, listC) -> +// for (i in 0..accDegs.lastIndex) { +// if (accDegs[i] > listDegs.getOrElse(i) { 0 }) return@reduce accDegs to accC +// if (accDegs[i] < listDegs.getOrElse(i) { 0 }) return@reduce listDegs to listC +// } +// if (accDegs.size < listDegs.size) listDegs to listC else accDegs to accC +// } +// +// var thisCoefficients = coefficients.toMutableMap() +// val otherCoefficients = other.coefficients +// +// var (thisLeadingTermDegs, thisLeadingTermC) = thisCoefficients.leadingTerm() +// val (otherLeadingTermDegs, otherLeadingTermC) = otherCoefficients.leadingTerm() +// +// while ( +// thisLeadingTermDegs.size >= otherLeadingTermDegs.size && +// (0..otherLeadingTermDegs.lastIndex).all { thisLeadingTermDegs[it] >= otherLeadingTermDegs[it] } +// ) { +// val multiplierDegs = +// thisLeadingTermDegs +// .mapIndexed { index, deg -> deg - otherLeadingTermDegs.getOrElse(index) { 0 } } +// .cleanUp() +// val multiplierC = thisLeadingTermC / otherLeadingTermC +// +// for ((degs, t) in otherCoefficients) { +// val productDegs = +// (0..max(degs.lastIndex, multiplierDegs.lastIndex)) +// .map { degs.getOrElse(it) { 0 } + multiplierDegs.getOrElse(it) { 0 } } +// .cleanUp() +// val productC = t * multiplierC +// thisCoefficients[productDegs] = +// if (productDegs in thisCoefficients) thisCoefficients[productDegs]!! - productC else -productC +// } +// +// thisCoefficients = thisCoefficients.filterValues { it.isNotZero() }.toMutableMap() +// +// if (thisCoefficients.isEmpty()) +// return Polynomial(thisCoefficients, toCheckInput = false) +// +// val t = thisCoefficients.leadingTerm() +// thisLeadingTermDegs = t.first +// thisLeadingTermC = t.second +// } +// +// return Polynomial(thisCoefficients, toCheckInput = false) +//} +// +//infix fun > Polynomial.divrem(other: Polynomial): Polynomial.Companion.DividingResult { +// if (other.isZero()) throw ArithmeticException("/ by zero") +// if (isZero()) return Polynomial.Companion.DividingResult(this, this) +// +// fun Map, T>.leadingTerm() = +// this +// .asSequence() +// .map { Pair(it.key, it.value) } +// .reduce { (accDegs, accC), (listDegs, listC) -> +// for (i in 0..accDegs.lastIndex) { +// if (accDegs[i] > listDegs.getOrElse(i) { 0 }) return@reduce accDegs to accC +// if (accDegs[i] < listDegs.getOrElse(i) { 0 }) return@reduce listDegs to listC +// } +// if (accDegs.size < listDegs.size) listDegs to listC else accDegs to accC +// } +// +// var thisCoefficients = coefficients.toMutableMap() +// val otherCoefficients = other.coefficients +// val quotientCoefficients = HashMap, T>() +// +// var (thisLeadingTermDegs, thisLeadingTermC) = thisCoefficients.leadingTerm() +// val (otherLeadingTermDegs, otherLeadingTermC) = otherCoefficients.leadingTerm() +// +// while ( +// thisLeadingTermDegs.size >= otherLeadingTermDegs.size && +// (0..otherLeadingTermDegs.lastIndex).all { thisLeadingTermDegs[it] >= otherLeadingTermDegs[it] } +// ) { +// val multiplierDegs = +// thisLeadingTermDegs +// .mapIndexed { index, deg -> deg - otherLeadingTermDegs.getOrElse(index) { 0 } } +// .cleanUp() +// val multiplierC = thisLeadingTermC / otherLeadingTermC +// +// quotientCoefficients[multiplierDegs] = multiplierC +// +// for ((degs, t) in otherCoefficients) { +// val productDegs = +// (0..max(degs.lastIndex, multiplierDegs.lastIndex)) +// .map { degs.getOrElse(it) { 0 } + multiplierDegs.getOrElse(it) { 0 } } +// .cleanUp() +// val productC = t * multiplierC +// thisCoefficients[productDegs] = +// if (productDegs in thisCoefficients) thisCoefficients[productDegs]!! - productC else -productC +// } +// +// thisCoefficients = thisCoefficients.filterValues { it.isNotZero() }.toMutableMap() +// +// if (thisCoefficients.isEmpty()) +// return Polynomial.Companion.DividingResult( +// Polynomial(quotientCoefficients, toCheckInput = false), +// Polynomial(thisCoefficients, toCheckInput = false) +// ) +// +// val t = thisCoefficients.leadingTerm() +// thisLeadingTermDegs = t.first +// thisLeadingTermC = t.second +// } +// +// return Polynomial.Companion.DividingResult( +// Polynomial(quotientCoefficients, toCheckInput = false), +// Polynomial(thisCoefficients, toCheckInput = false) +// ) +//} +// +//// endregion + +// 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 NumberedPolynomial.substitute(args: Map): NumberedPolynomial = Double.algebra { + val acc = LinkedHashMap, Double>(coefficients.size) + for ((degs, c) in coefficients) { + val newDegs = degs.mapIndexed { index, deg -> if (index !in args) deg else 0u }.cleanUp() + val newC = args.entries.fold(c) { product, (variable, substitutor) -> + val deg = degs.getOrElse(variable) { 0u } + if (deg == 0u) product else product * substitutor.pow(deg.toInt()) + } + if (newDegs !in acc) acc[newDegs] = c + else acc[newDegs] = acc[newDegs]!! + c + } + return NumberedPolynomial(acc) +} + +/** + * 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 NumberedPolynomial.substitute(ring: Ring, args: Map): NumberedPolynomial = ring { + val acc = LinkedHashMap, C>(coefficients.size) + for ((degs, c) in coefficients) { + val newDegs = degs.mapIndexed { index, deg -> if (index !in args) deg else 0u }.cleanUp() + val newC = args.entries.fold(c) { product, (variable, substitutor) -> + val deg = degs.getOrElse(variable) { 0u } + if (deg == 0u) product else product * power(substitutor, deg) + } + if (newDegs !in acc) acc[newDegs] = c + else acc[newDegs] = acc[newDegs]!! + c + } + return NumberedPolynomial(acc) +} + +// TODO: (Waiting for hero) Replace with optimisation: the [result] may be unboxed, and all operations may be performed +// as soon as possible on it +@JvmName("substitutePolynomial") +public fun NumberedPolynomial.substitute(ring: Ring, args: Map>) : NumberedPolynomial = TODO() /*ring.numberedPolynomial { + val acc = LinkedHashMap, NumberedPolynomial>(coefficients.size) + for ((degs, c) in coefficients) { + val newDegs = degs.mapIndexed { index, deg -> if (index !in args) deg else 0u }.cleanUp() + val newC = args.entries.fold(c.asNumberedPolynomial()) { product, (variable, substitutor) -> + val deg = degs.getOrElse(variable) { 0u } + if (deg == 0u) product else product * power(substitutor, deg) + } + if (newDegs !in acc) acc[newDegs] = c.asNumberedPolynomial() + else acc[newDegs] = acc[newDegs]!! + c + } +}*/ + +/** + * Represent the polynomial as a regular context-less function. + */ +public fun > NumberedPolynomial.asFunction(ring: A): (Map) -> NumberedPolynomial = { substitute(ring, it) } + +/** + * Represent the polynomial as a regular context-less function. + */ +public fun > NumberedPolynomial.asPolynomialFunctionOver(ring: A): (Map>) -> NumberedPolynomial = { substitute(ring, it) } + +// endregion + +// region Algebraic derivative and antiderivative + +/** + * Returns algebraic derivative of received polynomial. + */ +@UnstableKMathAPI +public fun NumberedPolynomial.derivativeBy( + algebra: A, + variable: Int, +): Polynomial where A : Ring, A : NumericAlgebra = algebra { + TODO() +} + +/** + * Returns algebraic derivative of received polynomial. + */ +@UnstableKMathAPI +public fun NumberedPolynomial.derivativeBy( + algebra: A, + variables: IntArray, +): Polynomial where A : Ring, A : NumericAlgebra = algebra { + TODO() +} + +/** + * Returns algebraic derivative of received polynomial. + */ +@UnstableKMathAPI +public fun NumberedPolynomial.derivativeBy( + algebra: A, + variables: Collection, +): Polynomial where A : Ring, A : NumericAlgebra = derivativeBy(algebra, variables.toIntArray()) + +/** + * Create a polynomial witch represents indefinite integral version of this polynomial + */ +@UnstableKMathAPI +public fun NumberedPolynomial.antiderivativeBy( + algebra: A, + variable: Int, +): Polynomial where A : Field, A : NumericAlgebra = algebra { + TODO() +} + +/** + * Create a polynomial witch represents indefinite integral version of this polynomial + */ +@UnstableKMathAPI +public fun NumberedPolynomial.antiderivativeBy( + algebra: A, + variables: IntArray, +): Polynomial where A : Field, A : NumericAlgebra = algebra { + TODO() +} + +/** + * Create a polynomial witch represents indefinite integral version of this polynomial + */ +@UnstableKMathAPI +public fun NumberedPolynomial.antiderivativeBy( + algebra: A, + variables: Collection, +): Polynomial where A : Field, A : NumericAlgebra = antiderivativeBy(algebra, variables.toIntArray()) + +// endregion \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt index 1a3eb7874..4d99b3a45 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt @@ -121,7 +121,7 @@ public fun Polynomial.antiderivative( ): Polynomial where A : Field, A : NumericAlgebra = algebra { val integratedCoefficients = buildList(coefficients.size + 1) { add(zero) - coefficients.forEachIndexed{ index, t -> add(t / (number(index) + one)) } + coefficients.forEachIndexed{ index, t -> add(t / number(index + 1)) } } Polynomial(integratedCoefficients) } @@ -136,4 +136,6 @@ public fun > Polynomial.integrate( ): C = algebra { val integral = antiderivative(algebra) integral.substitute(algebra, range.endInclusive) - integral.substitute(algebra, range.start) -} \ No newline at end of file +} + +// endregion \ No newline at end of file