kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

Added support for all polynomials. But standard utilities still are not fully implemented.

Browse Source
This commit is contained in:
Gleb Minaev 2022-03-13 03:27:00 +03:00
parent ffea8cc223
commit 843d63c76a
10 changed files with 2911 additions and 95 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -17,6 +17,9 @@ public interface AbstractPolynomial<C>
/** /**
* Abstraction of ring of polynomials of type [P] over ring of constants of type [C]. * 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") @Suppress("INAPPLICABLE_JVM_NAME", "PARAMETER_NAME_CHANGED_ON_OVERRIDE")
public interface AbstractPolynomialSpace<C, P: AbstractPolynomial<C>> : Ring<P> { public interface AbstractPolynomialSpace<C, P: AbstractPolynomial<C>> : Ring<P> {
@ -307,48 +310,4 @@ public interface AbstractPolynomialSpace<C, P: AbstractPolynomial<C>> : Ring<P>
override fun add(left: P, right: P): P = left + right override fun add(left: P, right: P): P = left + right
override fun multiply(left: P, right: P): P = left * right override fun multiply(left: P, right: P): P = left * right
// endregion // endregion
public companion object {
// TODO: All of this should be moved to algebraic structures' place for utilities
// TODO: Move receiver to context receiver
/**
* Multiplication of element and integer.
*
* @receiver the multiplicand.
* @param other the multiplier.
* @return the difference.
*/
internal tailrec fun <C> Group<C>.optimizedMultiply(arg: C, other: Int): C =
when {
other == 0 -> zero
other == 1 -> arg
other == -1 -> -arg
other % 2 == 0 -> optimizedMultiply(arg + arg, other / 2)
other % 2 == 1 -> optimizedAddMultiplied(arg, arg + arg, other / 2)
other % 2 == -1 -> optimizedAddMultiplied(-arg, arg + arg, other / 2)
else -> error("Error in multiplication group instant by integer: got reminder by division by 2 different from 0, 1 and -1")
}
// TODO: Move receiver to context receiver
/**
* Adds product of [arg] and [multiplier] to [base].
*
* @receiver the algebra to provide multiplication.
* @param base the augend.
* @param arg the multiplicand.
* @param multiplier the multiplier.
* @return sum of the augend [base] and product of the multiplicand [arg] and the multiplier [multiplier].
* @author Gleb Minaev
*/
internal tailrec fun <C> Group<C>.optimizedAddMultiplied(base: C, arg: C, multiplier: Int): C =
when {
multiplier == 0 -> base
multiplier == 1 -> base + arg
multiplier == -1 -> base - arg
multiplier % 2 == 0 -> optimizedAddMultiplied(base, arg + arg, multiplier / 2)
multiplier % 2 == 1 -> optimizedAddMultiplied(base + arg, arg + arg, multiplier / 2)
multiplier % 2 == -1 -> optimizedAddMultiplied(base + arg, arg + arg, multiplier / 2)
else -> error("Error in multiplication group instant by integer: got reminder by division by 2 different from 0, 1 and -1")
}
}
} }

2
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/AbstractRationalFunction.kt
View File

@ -5,8 +5,6 @@
package space.kscience.kmath.functions 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 space.kscience.kmath.operations.*
import kotlin.js.JsName import kotlin.js.JsName
import kotlin.jvm.JvmName import kotlin.jvm.JvmName

927
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt Normal file
View File

@ -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<C>
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<Map<Variable, UInt>, C>
) : AbstractPolynomial<C> {
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<Variable, UInt>.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 <C, A: Ring<C>> LabeledPolynomial(coefs: Map<Map<Variable, UInt>, C>, toCheckInput: Boolean = true): LabeledPolynomial<C> {
// if (!toCheckInput) return LabeledPolynomial<C>(coefs)
//
// // Map for cleaned coefficients.
// val fixedCoefs = mutableMapOf<Map<Variable, UInt>, 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<C>(
// 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 <C, A: Ring<C>> LabeledPolynomial(pairs: Collection<Pair<Map<Variable, UInt>, C>>, toCheckInput: Boolean): LabeledPolynomial<C> {
// if (!toCheckInput) return LabeledPolynomial<C>(pairs.toMap())
//
// // Map for cleaned coefficients.
// val fixedCoefs = mutableMapOf<Map<Variable, UInt>, 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<C>(
// 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 <C, A: Ring<C>> LabeledPolynomial(vararg pairs: Pair<Map<Variable, UInt>, C>, toCheckInput: Boolean): LabeledPolynomial<C> {
// if (!toCheckInput) return LabeledPolynomial<C>(pairs.toMap())
//
// // Map for cleaned coefficients.
// val fixedCoefs = mutableMapOf<Map<Variable, UInt>, 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<C>(
// 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 <C, A: Ring<C>> LabeledPolynomial(coefs: Map<Map<Variable, UInt>, C>): LabeledPolynomial<C> = 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 <C, A: Ring<C>> LabeledPolynomial(pairs: Collection<Pair<Map<Variable, UInt>, C>>): LabeledPolynomial<C> = 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 <C, A: Ring<C>> LabeledPolynomial(vararg pairs: Pair<Map<Variable, UInt>, C>): LabeledPolynomial<C> = 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<C, A>)
//internal fun <C, A: Ring<C>> LabeledPolynomial(coefs: Map<Map<Variable, UInt>, C>, toCheckInput: Boolean = true): LabeledPolynomial<C> {
// if (!toCheckInput) return LabeledPolynomial<C>(coefs)
//
// // Map for cleaned coefficients.
// val fixedCoefs = mutableMapOf<Map<Variable, UInt>, 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<C>(
// 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<C, A>)
//internal fun <C, A: Ring<C>> LabeledPolynomial(pairs: Collection<Pair<Map<Variable, UInt>, C>>, toCheckInput: Boolean): LabeledPolynomial<C> {
// if (!toCheckInput) return LabeledPolynomial<C>(pairs.toMap())
//
// // Map for cleaned coefficients.
// val fixedCoefs = mutableMapOf<Map<Variable, UInt>, 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<C>(
// 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<C, A>)
//internal fun <C, A: Ring<C>> LabeledPolynomial(vararg pairs: Pair<Map<Variable, UInt>, C>, toCheckInput: Boolean): LabeledPolynomial<C> {
// if (!toCheckInput) return LabeledPolynomial<C>(pairs.toMap())
//
// // Map for cleaned coefficients.
// val fixedCoefs = mutableMapOf<Map<Variable, UInt>, 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<C>(
// 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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial(coefs: Map<Map<Variable, UInt>, C>): LabeledPolynomial<C> = 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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial(pairs: Collection<Pair<Map<Variable, UInt>, C>>): LabeledPolynomial<C> = 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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial(vararg pairs: Pair<Map<Variable, UInt>, C>): LabeledPolynomial<C> = LabeledPolynomial(*pairs, toCheckInput = true)
//
//fun <C> C.asLabeledPolynomial() : LabeledPolynomial<C> = LabeledPolynomial<C>(mapOf(emptyMap<Variable, UInt>() to this))
//
//context(A)
//fun <C, A: Ring<C>> Variable.asLabeledPolynomial() : LabeledPolynomial<C> = LabeledPolynomial<C>(mapOf(mapOf<Variable, UInt>(this to 1U) to one))
//
//context(LabeledPolynomialSpace<C, A>)
//fun <C, A: Ring<C>> Variable.asLabeledPolynomial() : LabeledPolynomial<C> = LabeledPolynomial<C>(mapOf(mapOf<Variable, UInt>(this to 1U) to ring.one))
//
//context(A)
//fun <C, A: Ring<C>> Variable.asLabeledPolynomial(c: C) : LabeledPolynomial<C> =
// if(c.isZero()) LabeledPolynomial<C>(emptyMap())
// else LabeledPolynomial<C>(mapOf(mapOf<Variable, UInt>(this to 1U) to c))
//
//context(LabeledPolynomialSpace<C, A>)
//fun <C, A: Ring<C>> Variable.asLabeledPolynomial(c: C) : LabeledPolynomial<C> =
// if(c.isZero()) zero
// else LabeledPolynomial<C>(mapOf(mapOf<Variable, UInt>(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<C, A : Ring<C>>(
public val ring: A,
) : AbstractPolynomialSpace<C, LabeledPolynomial<C>> {
// region Constant-integer relation
@JvmName("constantIntPlus")
public override operator fun C.plus(other: Int): C = ring { optimizedAddMultiplied(this@plus, one, other) }
@JvmName("constantIntMinus")
public override operator fun C.minus(other: Int): C = ring { optimizedAddMultiplied(this@minus, one, -other) }
@JvmName("constantIntTimes")
public override operator fun C.times(other: Int): C = ring { optimizedMultiply(this@times, other) }
// endregion
// 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<C> =
if (other == 0) LabeledPolynomial<C>(mapOf(
mapOf(this@plus to 1U) to ring.one,
))
else LabeledPolynomial<C>(mapOf(
mapOf(this@plus to 1U) to ring.one,
emptyMap<Variable, UInt>() to ring.one * other,
))
public operator fun Variable.minus(other: Int): LabeledPolynomial<C> =
if (other == 0) LabeledPolynomial<C>(mapOf(
mapOf(this@minus to 1U) to -ring.one,
))
else LabeledPolynomial<C>(mapOf(
mapOf(this@minus to 1U) to -ring.one,
emptyMap<Variable, UInt>() to ring.one * other,
))
public operator fun Variable.times(other: Int): LabeledPolynomial<C> =
if (other == 0) zero
else LabeledPolynomial<C>(mapOf(
mapOf(this to 1U) to ring.one * other,
))
// endregion
// region Integer-variable relation
public operator fun Int.plus(other: Variable): LabeledPolynomial<C> =
if (this == 0) LabeledPolynomial<C>(mapOf(
mapOf(other to 1U) to ring.one,
))
else LabeledPolynomial<C>(mapOf(
mapOf(other to 1U) to ring.one,
emptyMap<Variable, UInt>() to ring.one * this@plus,
))
public operator fun Int.minus(other: Variable): LabeledPolynomial<C> =
if (this == 0) LabeledPolynomial<C>(mapOf(
mapOf(other to 1U) to -ring.one,
))
else LabeledPolynomial<C>(mapOf(
mapOf(other to 1U) to -ring.one,
emptyMap<Variable, UInt>() to ring.one * this@minus,
))
public operator fun Int.times(other: Variable): LabeledPolynomial<C> =
if (this == 0) zero
else LabeledPolynomial<C>(mapOf(
mapOf(other to 1U) to ring.one * this@times,
))
// endregion
// region Polynomial-integer relation
public override operator fun LabeledPolynomial<C>.plus(other: Int): LabeledPolynomial<C> =
if (other == 0) this
else
LabeledPolynomial(
coefficients
.toMutableMap()
.apply {
val degs = emptyMap<Variable, UInt>()
val result = getOrElse(degs) { ring.zero } + other
if (result.isZero()) remove(degs)
else this[degs] = result
}
)
public override operator fun LabeledPolynomial<C>.minus(other: Int): LabeledPolynomial<C> =
if (other == 0) this
else
LabeledPolynomial(
coefficients
.toMutableMap()
.apply {
val degs = emptyMap<Variable, UInt>()
val result = getOrElse(degs) { ring.zero } - other
if (result.isZero()) remove(degs)
else this[degs] = result
}
)
public override operator fun LabeledPolynomial<C>.times(other: Int): LabeledPolynomial<C> =
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<C>): LabeledPolynomial<C> =
if (this == 0) other
else
LabeledPolynomial(
other.coefficients
.toMutableMap()
.apply {
val degs = emptyMap<Variable, UInt>()
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<C>): LabeledPolynomial<C> =
if (this == 0) other
else
LabeledPolynomial(
other.coefficients
.toMutableMap()
.apply {
val degs = emptyMap<Variable, UInt>()
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<C>): LabeledPolynomial<C> =
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<C> =
if (isZero()) LabeledPolynomial<C>(mapOf(
mapOf(other to 1U) to ring.one,
))
else LabeledPolynomial<C>(mapOf(
mapOf(other to 1U) to ring.one,
emptyMap<Variable, UInt>() to this@plus,
))
public operator fun C.minus(other: Variable): LabeledPolynomial<C> =
if (isZero()) LabeledPolynomial<C>(mapOf(
mapOf(other to 1U) to -ring.one,
))
else LabeledPolynomial<C>(mapOf(
mapOf(other to 1U) to -ring.one,
emptyMap<Variable, UInt>() to this@minus,
))
public operator fun C.times(other: Variable): LabeledPolynomial<C> =
if (isZero()) zero
else LabeledPolynomial<C>(mapOf(
mapOf(other to 1U) to this@times,
))
// endregion
// region Variable-constant relation
public operator fun Variable.plus(other: C): LabeledPolynomial<C> =
if (other.isZero()) LabeledPolynomial<C>(mapOf(
mapOf(this@plus to 1U) to ring.one,
))
else LabeledPolynomial<C>(mapOf(
mapOf(this@plus to 1U) to ring.one,
emptyMap<Variable, UInt>() to other,
))
public operator fun Variable.minus(other: C): LabeledPolynomial<C> =
if (other.isZero()) LabeledPolynomial<C>(mapOf(
mapOf(this@minus to 1U) to -ring.one,
))
else LabeledPolynomial<C>(mapOf(
mapOf(this@minus to 1U) to -ring.one,
emptyMap<Variable, UInt>() to other,
))
public operator fun Variable.times(other: C): LabeledPolynomial<C> =
if (other.isZero()) zero
else LabeledPolynomial<C>(mapOf(
mapOf(this@times to 1U) to other,
))
// endregion
// region Constant-polynomial relation
override operator fun C.plus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
if (this.isZero()) other
else with(other.coefficients) {
if (isEmpty()) LabeledPolynomial<C>(mapOf(emptyMap<Variable, UInt>() to this@plus))
else LabeledPolynomial<C>(
toMutableMap()
.apply {
val degs = emptyMap<Variable, UInt>()
val result = this@plus + getOrElse(degs) { ring.zero }
if (result.isZero()) remove(degs)
else this[degs] = result
}
)
}
override operator fun C.minus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
if (this.isZero()) other
else with(other.coefficients) {
if (isEmpty()) LabeledPolynomial<C>(mapOf(emptyMap<Variable, UInt>() to this@minus))
else LabeledPolynomial<C>(
toMutableMap()
.apply {
forEach { (degs, c) -> if(degs.isNotEmpty()) this[degs] = -c }
val degs = emptyMap<Variable, UInt>()
val result = this@minus - getOrElse(degs) { ring.zero }
if (result.isZero()) remove(degs)
else this[degs] = result
}
)
}
override operator fun C.times(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
if (this.isZero()) zero
else LabeledPolynomial<C>(
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<C>.plus(other: C): LabeledPolynomial<C> =
if (other.isZero()) this
else with(coefficients) {
if (isEmpty()) LabeledPolynomial<C>(mapOf(emptyMap<Variable, UInt>() to other))
else LabeledPolynomial<C>(
toMutableMap()
.apply {
val degs = emptyMap<Variable, UInt>()
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<C>.minus(other: C): LabeledPolynomial<C> =
if (other.isZero()) this
else with(coefficients) {
if (isEmpty()) LabeledPolynomial<C>(mapOf(emptyMap<Variable, UInt>() to other))
else LabeledPolynomial<C>(
toMutableMap()
.apply {
forEach { (degs, c) -> if(degs.isNotEmpty()) this[degs] = -c }
val degs = emptyMap<Variable, UInt>()
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<C>.times(other: C): LabeledPolynomial<C> =
if (other.isZero()) zero
else LabeledPolynomial<C>(
coefficients
.applyAndRemoveZeros {
mapValues { (_, c) -> c * other }
}
)
// endregion
// region Variable-variable relation
public operator fun Variable.plus(other: Variable): LabeledPolynomial<C> =
if (this == other) LabeledPolynomial<C>(mapOf(
mapOf(this to 1U) to ring.one * 2
))
else LabeledPolynomial<C>(mapOf(
mapOf(this to 1U) to ring.one,
mapOf(other to 1U) to ring.one,
))
public operator fun Variable.minus(other: Variable): LabeledPolynomial<C> =
if (this == other) zero
else LabeledPolynomial<C>(mapOf(
mapOf(this to 1U) to ring.one,
mapOf(other to 1U) to -ring.one,
))
public operator fun Variable.times(other: Variable): LabeledPolynomial<C> =
if (this == other) LabeledPolynomial<C>(mapOf(
mapOf(this to 2U) to ring.one
))
else LabeledPolynomial<C>(mapOf(
mapOf(this to 1U, other to 1U) to ring.one,
))
// endregion
// region Variable-polynomial relation
public operator fun Variable.plus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
with(other.coefficients) {
if (isEmpty()) LabeledPolynomial<C>(mapOf(mapOf(this@plus to 1u) to ring.one))
else LabeledPolynomial<C>(
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<C>): LabeledPolynomial<C> =
with(other.coefficients) {
if (isEmpty()) LabeledPolynomial<C>(mapOf(mapOf(this@minus to 1u) to ring.one))
else LabeledPolynomial<C>(
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<C>): LabeledPolynomial<C> =
LabeledPolynomial<C>(
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<C>.plus(other: Variable): LabeledPolynomial<C> =
with(coefficients) {
if (isEmpty()) LabeledPolynomial<C>(mapOf(mapOf(other to 1u) to ring.one))
else LabeledPolynomial<C>(
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<C>.minus(other: Variable): LabeledPolynomial<C> =
with(coefficients) {
if (isEmpty()) LabeledPolynomial<C>(mapOf(mapOf(other to 1u) to ring.one))
else LabeledPolynomial<C>(
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<C>.times(other: Variable): LabeledPolynomial<C> =
LabeledPolynomial<C>(
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<C>.unaryMinus(): LabeledPolynomial<C> =
LabeledPolynomial<C>(
coefficients.mapValues { -it.value }
)
/**
* Returns sum of the polynomials.
*/
override operator fun LabeledPolynomial<C>.plus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
LabeledPolynomial<C>(
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<C>.minus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
LabeledPolynomial<C>(
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<C>.times(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
when {
isZero() -> zero
other.isZero() -> zero
else -> LabeledPolynomial<C>(
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<C> = LabeledPolynomial<C>(mapOf(emptyMap<Variable, UInt>() to ring.zero))
override val one: LabeledPolynomial<C> = LabeledPolynomial<C>(mapOf(emptyMap<Variable, UInt>() to ring.one))
// TODO: Docs
@Suppress("EXTENSION_SHADOWED_BY_MEMBER", "CovariantEquals")
override fun LabeledPolynomial<C>.equals(other: LabeledPolynomial<C>): 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<C>.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<C>.degrees: Map<Variable, UInt>
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<C>.variables: Set<Variable>
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<C>.countOfVariables: Int get() = variables.size
/**
* Checks if the instant is constant polynomial (of degree no more than 0) over considered ring.
*/
override fun LabeledPolynomial<C>.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<C>.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<C>.asConstantOrNull(): C? =
with(coefficients) {
if(isConstant()) getOrElse(emptyMap()) { ring.zero }
else null
}
// @Suppress("NOTHING_TO_INLINE")
// public inline fun LabeledPolynomial<C>.substitute(argument: Map<Variable, C>): LabeledPolynomial<C> = this.substitute(ring, argument)
// @Suppress("NOTHING_TO_INLINE")
// @JvmName("substitutePolynomial")
// public inline fun LabeledPolynomial<C>.substitute(argument: Map<Variable, LabeledPolynomial<C>>): LabeledPolynomial<C> = this.substitute(ring, argument)
//
// @Suppress("NOTHING_TO_INLINE")
// public inline fun LabeledPolynomial<C>.asFunction(): (Map<Variable, C>) -> LabeledPolynomial<C> = { this.substitute(ring, it) }
// @Suppress("NOTHING_TO_INLINE")
// public inline fun LabeledPolynomial<C>.asFunctionOnConstants(): (Map<Variable, C>) -> LabeledPolynomial<C> = { this.substitute(ring, it) }
// @Suppress("NOTHING_TO_INLINE")
// public inline fun LabeledPolynomial<C>.asFunctionOnPolynomials(): (Map<Variable, LabeledPolynomial<C>>) -> LabeledPolynomial<C> = { this.substitute(ring, it) }
//
// @Suppress("NOTHING_TO_INLINE")
// public inline operator fun LabeledPolynomial<C>.invoke(argument: Map<Variable, C>): LabeledPolynomial<C> = this.substitute(ring, argument)
// @Suppress("NOTHING_TO_INLINE")
// @JvmName("invokePolynomial")
// public inline operator fun LabeledPolynomial<C>.invoke(argument: Map<Variable, LabeledPolynomial<C>>): LabeledPolynomial<C> = this.substitute(ring, argument)
// endregion
// region Legacy
@Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
override inline fun add(left: LabeledPolynomial<C>, right: LabeledPolynomial<C>): LabeledPolynomial<C> = left + right
@Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
override inline fun multiply(left: LabeledPolynomial<C>, right: LabeledPolynomial<C>): LabeledPolynomial<C> = left * right
// endregion
// region Utilities
// TODO: Move to region internal utilities with context receiver
internal fun MutableMap<Map<Variable, UInt>, C>.applyAndRemoveZeros(block: MutableMap<Map<Variable, UInt>, C>.() -> Unit) : MutableMap<Map<Variable, UInt>, C> {
contract {
callsInPlace(block, InvocationKind.EXACTLY_ONCE)
}
block()
for ((degs, c) in this) if (c.isZero()) this.remove(degs)
return this
}
internal fun Map<Map<Variable, UInt>, C>.applyAndRemoveZeros(block: MutableMap<Map<Variable, UInt>, C>.() -> Unit) : Map<Map<Variable, UInt>, C> =
toMutableMap().applyAndRemoveZeros(block)
@OptIn(ExperimentalTypeInference::class)
internal fun buildCoefficients(@BuilderInference builderAction: MutableMap<Map<Variable, UInt>, C>.() -> Unit): Map<Map<Variable, UInt>, C> {
contract { callsInPlace(builderAction, InvocationKind.EXACTLY_ONCE) }
return buildMap {
builderAction()
for ((degs, c) in this) if (c.isZero()) this.remove(degs)
}
}
// endregion
}

689
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/NumberedPolynomial.kt Normal file
View File

@ -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<C>
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<List<UInt>, C>
) : AbstractPolynomial<C> {
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<UInt>.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 <C, A: Ring<C>> NumberedPolynomial(coefs: Map<List<UInt>, C>, toCheckInput: Boolean): NumberedPolynomial<C> {
// if (!toCheckInput) return NumberedPolynomial<C>(coefs)
//
// val fixedCoefs = mutableMapOf<List<UInt>, 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<C>(
// 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 <C, A: Ring<C>> NumberedPolynomial(pairs: Collection<Pair<List<UInt>, C>>, toCheckInput: Boolean): NumberedPolynomial<C> {
// if (!toCheckInput) return NumberedPolynomial(pairs.toMap())
//
// val fixedCoefs = mutableMapOf<List<UInt>, 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<C>(
// 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 <C, A: Ring<C>> NumberedPolynomial(vararg pairs: Pair<List<UInt>, C>, toCheckInput: Boolean): NumberedPolynomial<C> =
// 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 <C, A: Ring<C>> NumberedPolynomial(coefs: Map<List<UInt>, 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 <C, A: Ring<C>> NumberedPolynomial(pairs: Collection<Pair<List<UInt>, 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 <C, A: Ring<C>> NumberedPolynomial(vararg pairs: Pair<List<UInt>, C>) = NumberedPolynomial(*pairs, toCheckInput = true)
//
//context(NumberedPolynomialSpace<C, A>)
//@Suppress("FunctionName")
//internal fun <C, A: Ring<C>> NumberedPolynomial(coefs: Map<List<UInt>, C>, toCheckInput: Boolean): NumberedPolynomial<C> {
// if (!toCheckInput) return NumberedPolynomial(coefs)
//
// val fixedCoefs = mutableMapOf<List<UInt>, 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<C>(
// 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<C, A>)
//@Suppress("FunctionName")
//internal fun <C, A: Ring<C>> NumberedPolynomial(pairs: Collection<Pair<List<UInt>, C>>, toCheckInput: Boolean): NumberedPolynomial<C> {
// if (!toCheckInput) return NumberedPolynomial(pairs.toMap())
//
// val fixedCoefs = mutableMapOf<List<UInt>, 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<C>(
// 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<C, A>)
//@Suppress("FunctionName")
//internal fun <C, A: Ring<C>> NumberedPolynomial(vararg pairs: Pair<List<UInt>, C>, toCheckInput: Boolean): NumberedPolynomial<C> =
// 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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial(coefs: Map<List<UInt>, 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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial(pairs: Collection<Pair<List<UInt>, 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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial(vararg pairs: Pair<List<UInt>, C>) = NumberedPolynomial(*pairs, toCheckInput = true)
public fun <C, A: Ring<C>> C.asNumberedPolynomial() : NumberedPolynomial<C> = NumberedPolynomial<C>(mapOf(emptyList<UInt>() 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<C, A : Ring<C>>(
public val ring: A,
) : AbstractPolynomialSpace<C, NumberedPolynomial<C>> {
// region Constant-integer relation
@JvmName("constantIntPlus")
public override operator fun C.plus(other: Int): C = ring { optimizedAddMultiplied(this@plus, one, other) }
@JvmName("constantIntMinus")
public override operator fun C.minus(other: Int): C = ring { optimizedAddMultiplied(this@minus, one, -other) }
@JvmName("constantIntTimes")
public override operator fun C.times(other: Int): C = ring { optimizedMultiply(this@times, other) }
// endregion
// 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<C>.plus(other: Int): NumberedPolynomial<C> =
if (other == 0) this
else
NumberedPolynomial(
coefficients
.toMutableMap()
.apply {
val degs = emptyList<UInt>()
val result = getOrElse(degs) { ring.zero } + other
if (result.isZero()) remove(degs)
else this[degs] = result
}
)
public override operator fun NumberedPolynomial<C>.minus(other: Int): NumberedPolynomial<C> =
if (other == 0) this
else
NumberedPolynomial(
coefficients
.toMutableMap()
.apply {
val degs = emptyList<UInt>()
val result = getOrElse(degs) { ring.zero } - other
if (result.isZero()) remove(degs)
else this[degs] = result
}
)
public override operator fun NumberedPolynomial<C>.times(other: Int): NumberedPolynomial<C> =
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<C>): NumberedPolynomial<C> =
if (this == 0) other
else
NumberedPolynomial(
other.coefficients
.toMutableMap()
.apply {
val degs = emptyList<UInt>()
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<C>): NumberedPolynomial<C> =
if (this == 0) other
else
NumberedPolynomial(
other.coefficients
.toMutableMap()
.apply {
val degs = emptyList<UInt>()
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<C>): NumberedPolynomial<C> =
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<C>): NumberedPolynomial<C> =
if (this.isZero()) other
else with(other.coefficients) {
if (isEmpty()) NumberedPolynomial<C>(mapOf(emptyList<UInt>() to this@plus))
else NumberedPolynomial<C>(
toMutableMap()
.apply {
val degs = emptyList<UInt>()
val result = this@plus + getOrElse(degs) { ring.zero }
if (result.isZero()) remove(degs)
else this[degs] = result
}
)
}
override operator fun C.minus(other: NumberedPolynomial<C>): NumberedPolynomial<C> =
if (this.isZero()) -other
else with(other.coefficients) {
if (isEmpty()) NumberedPolynomial<C>(mapOf(listOf<UInt>() to this@minus))
else NumberedPolynomial<C>(
toMutableMap()
.apply {
forEach { (degs, c) -> if(degs.isNotEmpty()) this[degs] = -c }
val degs = emptyList<UInt>()
val result = this@minus - getOrElse(degs) { ring.zero }
if (result.isZero()) remove(degs)
else this[degs] = result
}
)
}
override operator fun C.times(other: NumberedPolynomial<C>): NumberedPolynomial<C> =
if (this.isZero()) zero
else NumberedPolynomial<C>(
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<C>.plus(other: C): NumberedPolynomial<C> =
if (other.isZero()) this
else with(coefficients) {
if (isEmpty()) NumberedPolynomial<C>(mapOf(listOf<UInt>() to other))
else NumberedPolynomial<C>(
toMutableMap()
.apply {
val degs = emptyList<UInt>()
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<C>.minus(other: C): NumberedPolynomial<C> =
if (other.isZero()) this
else with(coefficients) {
if (isEmpty()) NumberedPolynomial<C>(mapOf(listOf<UInt>() to other))
else NumberedPolynomial<C>(
toMutableMap()
.apply {
val degs = emptyList<UInt>()
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<C>.times(other: C): NumberedPolynomial<C> =
if (other.isZero()) zero
else NumberedPolynomial<C>(
coefficients
.applyAndRemoveZeros {
mapValues { (_, c) -> c * other }
}
)
// endregion
// region Polynomial-polynomial relation
/**
* Returns negation of the polynomial.
*/
override fun NumberedPolynomial<C>.unaryMinus(): NumberedPolynomial<C> =
NumberedPolynomial<C>(
coefficients.mapValues { -it.value }
)
/**
* Returns sum of the polynomials.
*/
override operator fun NumberedPolynomial<C>.plus(other: NumberedPolynomial<C>): NumberedPolynomial<C> =
NumberedPolynomial<C>(
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<C>.minus(other: NumberedPolynomial<C>): NumberedPolynomial<C> =
NumberedPolynomial<C>(
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<C>.times(other: NumberedPolynomial<C>): NumberedPolynomial<C> =
when {
isZero() -> zero
other.isZero() -> zero
else ->
NumberedPolynomial<C>(
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<C>.isZero(): Boolean = coefficients.values.all { it.isZero() }
public override fun NumberedPolynomial<C>.isNotZero(): Boolean = coefficients.values.any { it.isNotZero() }
public override fun NumberedPolynomial<C>.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<C>.isNotOne(): Boolean = !isOne()
public override fun NumberedPolynomial<C>.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<C>.isNotMinusOne(): Boolean = !isMinusOne()
override val zero: NumberedPolynomial<C> = NumberedPolynomial<C>(emptyMap())
override val one: NumberedPolynomial<C> =
NumberedPolynomial<C>(
mapOf(
listOf<UInt>() to ring.one // 1 * x_1^0 * x_2^0 * ...
)
)
// TODO: Docs
@Suppress("EXTENSION_SHADOWED_BY_MEMBER", "CovariantEquals")
override fun NumberedPolynomial<C>.equals(other: NumberedPolynomial<C>): 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<C>.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<C>.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<C>.degrees: List<UInt>
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<C>.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<C>.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<C>.asConstantOrNull(): C? =
with(coefficients) {
if(isConstant()) getOrElse(emptyList()) { ring.zero }
else null
}
@Suppress("NOTHING_TO_INLINE")
public inline fun NumberedPolynomial<C>.substitute(argument: Map<Int, C>): NumberedPolynomial<C> = this.substitute(ring, argument)
@Suppress("NOTHING_TO_INLINE")
@JvmName("substitutePolynomial")
public inline fun NumberedPolynomial<C>.substitute(argument: Map<Int, NumberedPolynomial<C>>): NumberedPolynomial<C> = this.substitute(ring, argument)
@Suppress("NOTHING_TO_INLINE")
public inline fun NumberedPolynomial<C>.asFunction(): (Map<Int, C>) -> NumberedPolynomial<C> = { this.substitute(ring, it) }
@Suppress("NOTHING_TO_INLINE")
public inline fun NumberedPolynomial<C>.asFunctionOnConstants(): (Map<Int, C>) -> NumberedPolynomial<C> = { this.substitute(ring, it) }
@Suppress("NOTHING_TO_INLINE")
public inline fun NumberedPolynomial<C>.asFunctionOnPolynomials(): (Map<Int, NumberedPolynomial<C>>) -> NumberedPolynomial<C> = { this.substitute(ring, it) }
@Suppress("NOTHING_TO_INLINE")
public inline operator fun NumberedPolynomial<C>.invoke(argument: Map<Int, C>): NumberedPolynomial<C> = this.substitute(ring, argument)
@Suppress("NOTHING_TO_INLINE")
@JvmName("invokePolynomial")
public inline operator fun NumberedPolynomial<C>.invoke(argument: Map<Int, NumberedPolynomial<C>>): NumberedPolynomial<C> = this.substitute(ring, argument)
// endregion
// region Legacy
@Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
override inline fun add(left: NumberedPolynomial<C>, right: NumberedPolynomial<C>): NumberedPolynomial<C> = left + right
@Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
override inline fun multiply(left: NumberedPolynomial<C>, right: NumberedPolynomial<C>): NumberedPolynomial<C> = left * right
// endregion
// region Utilities
// TODO: Move to region internal utilities with context receiver
internal fun MutableMap<List<UInt>, C>.applyAndRemoveZeros(block: MutableMap<List<UInt>, C>.() -> Unit) : MutableMap<List<UInt>, C> {
contract {
callsInPlace(block, InvocationKind.EXACTLY_ONCE)
}
block()
for ((degs, c) in this) if (c.isZero()) this.remove(degs)
return this
}
internal fun Map<List<UInt>, C>.applyAndRemoveZeros(block: MutableMap<List<UInt>, C>.() -> Unit) : Map<List<UInt>, C> =
toMutableMap().applyAndRemoveZeros(block)
@OptIn(ExperimentalTypeInference::class)
internal fun buildCoefficients(@BuilderInference builderAction: MutableMap<List<UInt>, C>.() -> Unit): Map<List<UInt>, C> {
contract { callsInPlace(builderAction, InvocationKind.EXACTLY_ONCE) }
return buildMap {
builderAction()
for ((degs, c) in this) if (c.isZero()) this.remove(degs)
}
}
// endregion
}

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

@ -5,43 +5,38 @@
package space.kscience.kmath.functions 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 space.kscience.kmath.operations.*
import kotlin.jvm.JvmName import kotlin.jvm.JvmName
import kotlin.math.max import kotlin.math.max
import kotlin.math.min 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. * @param coefficients constant is the leftmost coefficient.
*/ */
public class Polynomial<T>(public val coefficients: List<T>) : AbstractPolynomial<T> { public class Polynomial<T>(public val coefficients: List<T>) : AbstractPolynomial<T> {
override fun toString(): String = "Polynomial$coefficients" override fun toString(): String = "Polynomial$coefficients"
// public companion object {
// /**
// * Default name of variables used in string representations.
// *
// * @see Polynomial.toString
// */
// public var defaultVariableName: String = "x"
//
// /**
// * Represents result of division with remainder.
// */
// public data class DividingResult<C>(
// val quotient: Polynomial<C>,
// val reminder: Polynomial<C>
// )
// }
} }
// region Internal utilities
/** /**
* Represents internal [Polynomial] errors. * 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 // region Constructors and converters
@ -66,16 +61,16 @@ public fun <T> T.asPolynomial() : Polynomial<T> = Polynomial(listOf(this))
// endregion // 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 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 A type of underlying ring of constants. It's [Ring] of [C].
* @param ring underlying ring of constants of type [A]. * @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<C, A : Ring<C>>( public open class PolynomialSpace<C, A : Ring<C>>(
public val ring: A, public val ring: A,
) : AbstractPolynomialSpace<C, Polynomial<C>>{ ) : AbstractPolynomialSpace<C, Polynomial<C>> {
// region Constant-integer relation // region Constant-integer relation
@JvmName("constantIntPlus") @JvmName("constantIntPlus")
public override operator fun C.plus(other: Int): C = ring { optimizedAddMultiplied(this@plus, one, other) } public override operator fun C.plus(other: Int): C = ring { optimizedAddMultiplied(this@plus, one, other) }
@ -102,8 +97,14 @@ public open class PolynomialSpace<C, A : Ring<C>>(
coefficients coefficients
.toMutableList() .toMutableList()
.apply { .apply {
if (isEmpty()) this[0] = ring.zero + other val result = getOrElse(0) { ring.zero } + other
else this[0] = this[0]!! + other val isResultZero = result.isZero()
when {
size == 0 && !isResultZero -> add(result)
size > 1 || !isResultZero -> this[0] = result
else -> clear()
}
} }
) )
public override operator fun Polynomial<C>.minus(other: Int): Polynomial<C> = public override operator fun Polynomial<C>.minus(other: Int): Polynomial<C> =
@ -113,8 +114,14 @@ public open class PolynomialSpace<C, A : Ring<C>>(
coefficients coefficients
.toMutableList() .toMutableList()
.apply { .apply {
if (isEmpty()) this[0] = ring.zero - other val result = getOrElse(0) { ring.zero } - other
else this[0] = this[0]!! - other val isResultZero = result.isZero()
when {
size == 0 && !isResultZero -> add(result)
size > 1 || !isResultZero -> this[0] = result
else -> clear()
}
} }
) )
public override operator fun Polynomial<C>.times(other: Int): Polynomial<C> = public override operator fun Polynomial<C>.times(other: Int): Polynomial<C> =
@ -134,8 +141,14 @@ public open class PolynomialSpace<C, A : Ring<C>>(
other.coefficients other.coefficients
.toMutableList() .toMutableList()
.apply { .apply {
if (isEmpty()) this[0] = ring.zero + this@plus val result = this@plus + getOrElse(0) { ring.zero }
else this[0] = this[0]!! + this@plus 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<C>): Polynomial<C> = public override operator fun Int.minus(other: Polynomial<C>): Polynomial<C> =
@ -145,8 +158,16 @@ public open class PolynomialSpace<C, A : Ring<C>>(
other.coefficients other.coefficients
.toMutableList() .toMutableList()
.apply { .apply {
if (isEmpty()) this[0] = ring.zero - this@minus forEachIndexed { index, c -> if (index != 0) this[index] = -c }
else this[0] = this[0]!! - this@minus
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<C>): Polynomial<C> = public override operator fun Int.times(other: Polynomial<C>): Polynomial<C> =
@ -183,12 +204,20 @@ public open class PolynomialSpace<C, A : Ring<C>>(
// region Constant-polynomial relation // region Constant-polynomial relation
public override operator fun C.plus(other: Polynomial<C>): Polynomial<C> = public override operator fun C.plus(other: Polynomial<C>): Polynomial<C> =
with(other.coefficients) { if (this.isZero()) other
else with(other.coefficients) {
if (isEmpty()) Polynomial(listOf(this@plus)) if (isEmpty()) Polynomial(listOf(this@plus))
else Polynomial( else Polynomial(
toMutableList() toMutableList()
.apply { .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<C, A : Ring<C>>(
// listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1) // listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1)
// ) // )
public override operator fun C.minus(other: Polynomial<C>): Polynomial<C> = public override operator fun C.minus(other: Polynomial<C>): Polynomial<C> =
with(other.coefficients) { if (this.isZero()) other
else with(other.coefficients) {
if (isEmpty()) Polynomial(listOf(-this@minus)) if (isEmpty()) Polynomial(listOf(-this@minus))
else Polynomial( else Polynomial(
toMutableList() toMutableList()
.apply { .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<C, A : Ring<C>>(
// listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1) // listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1)
// ) // )
public override operator fun C.times(other: Polynomial<C>): Polynomial<C> = public override operator fun C.times(other: Polynomial<C>): Polynomial<C> =
Polynomial( if (this.isZero()) other
else Polynomial(
other.coefficients other.coefficients
// .subList(0, other.degree + 1) .subList(0, other.degree + 1)
.map { it * this } .map { it * this }
) )
// endregion // endregion
@ -224,7 +264,14 @@ public open class PolynomialSpace<C, A : Ring<C>>(
else Polynomial( else Polynomial(
toMutableList() toMutableList()
.apply { .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<C, A : Ring<C>>(
// listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1) // listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1)
// ) // )
public override operator fun Polynomial<C>.minus(other: C): Polynomial<C> = public override operator fun Polynomial<C>.minus(other: C): Polynomial<C> =
with(coefficients) { if (other.isZero()) this
else with(coefficients) {
if (isEmpty()) Polynomial(listOf(other)) if (isEmpty()) Polynomial(listOf(other))
else Polynomial( else Polynomial(
toMutableList() toMutableList()
.apply { .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<C, A : Ring<C>>(
// listOf(coefficients[0] - other) + coefficients.subList(1, degree + 1) // listOf(coefficients[0] - other) + coefficients.subList(1, degree + 1)
// ) // )
public override operator fun Polynomial<C>.times(other: C): Polynomial<C> = public override operator fun Polynomial<C>.times(other: C): Polynomial<C> =
Polynomial( if (other.isZero()) this
else Polynomial(
coefficients coefficients
// .subList(0, degree + 1) .subList(0, degree + 1)
.map { it * other } .map { it * other }
) )
// endregion // endregion
@ -283,8 +339,8 @@ public open class PolynomialSpace<C, A : Ring<C>>(
val thisDegree = degree val thisDegree = degree
val otherDegree = other.degree val otherDegree = other.degree
return when { return when {
thisDegree == -1 -> this thisDegree == -1 -> zero
otherDegree == -1 -> other otherDegree == -1 -> zero
else -> else ->
Polynomial( Polynomial(
(0..(thisDegree + otherDegree)) (0..(thisDegree + otherDegree))
@ -293,6 +349,7 @@ public open class PolynomialSpace<C, A : Ring<C>>(
.map { coefficients[it] * other.coefficients[d - it] } .map { coefficients[it] * other.coefficients[d - it] }
.reduce { acc, rational -> acc + rational } .reduce { acc, rational -> acc + rational }
} }
.run { subList(0, indexOfLast { it.isNotZero() } + 1) }
) )
} }
} }
@ -321,8 +378,8 @@ public open class PolynomialSpace<C, A : Ring<C>>(
} }
// endregion // endregion
// Not sure is it necessary...
// region Polynomial properties // region Polynomial properties
public override val Polynomial<C>.degree: Int get() = coefficients.indexOfLast { it != ring.zero } public override val Polynomial<C>.degree: Int get() = coefficients.indexOfLast { it != ring.zero }
public override fun Polynomial<C>.asConstantOrNull(): C? = public override fun Polynomial<C>.asConstantOrNull(): C? =
@ -354,8 +411,8 @@ public open class PolynomialSpace<C, A : Ring<C>>(
public inline operator fun Polynomial<C>.invoke(argument: C): C = this.substitute(ring, argument) public inline operator fun Polynomial<C>.invoke(argument: C): C = this.substitute(ring, argument)
@Suppress("NOTHING_TO_INLINE") @Suppress("NOTHING_TO_INLINE")
public inline operator fun Polynomial<C>.invoke(argument: Polynomial<C>): Polynomial<C> = this.substitute(ring, argument) public inline operator fun Polynomial<C>.invoke(argument: Polynomial<C>): Polynomial<C> = this.substitute(ring, argument)
// endregion
// endregion
} }
/** /**

38
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Variable.kt Normal file
View File

@ -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` &ndash; `"n"`, etc.
*/
public data class Variable (val name: String) : Comparable<Variable> {
/**
* 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)
}
}

51
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/algebraicStub.kt Normal file
View File

@ -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 <C> Group<C>.optimizedMultiply(arg: C, other: Int): C =
when {
other == 0 -> zero
other == 1 -> arg
other == -1 -> -arg
other % 2 == 0 -> optimizedMultiply(arg + arg, other / 2)
other % 2 == 1 -> optimizedAddMultiplied(arg, arg + arg, other / 2)
other % 2 == -1 -> optimizedAddMultiplied(-arg, arg + arg, other / 2)
else -> error("Error in multiplication group instant by integer: got reminder by division by 2 different from 0, 1 and -1")
}
// TODO: Move receiver to context receiver
/**
* Adds product of [arg] and [multiplier] to [base].
*
* @receiver the algebra to provide multiplication.
* @param base the augend.
* @param arg the multiplicand.
* @param multiplier the multiplier.
* @return sum of the augend [base] and product of the multiplicand [arg] and the multiplier [multiplier].
* @author Gleb Minaev
*/
internal tailrec fun <C> Group<C>.optimizedAddMultiplied(base: C, arg: C, multiplier: Int): C =
when {
multiplier == 0 -> base
multiplier == 1 -> base + arg
multiplier == -1 -> base - arg
multiplier % 2 == 0 -> optimizedAddMultiplied(base, arg + arg, multiplier / 2)
multiplier % 2 == 1 -> optimizedAddMultiplied(base + arg, arg + arg, multiplier / 2)
multiplier % 2 == -1 -> optimizedAddMultiplied(base + arg, arg + arg, multiplier / 2)
else -> error("Error in multiplication group instant by integer: got reminder by division by 2 different from 0, 1 and -1")
}

490
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/labeledPolynomialUtil.kt Normal file
View File

@ -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<C, A>)
//@Suppress("NOTHING_TO_INLINE")
//fun <C, A: Ring<C>> numberConstant(value: Int): C = ring { number<C>(value) }
//
//context(LabeledPolynomialSpace<C, A>)
//fun <C, A: Ring<C>> power(arg: C, pow: UInt): C = ring { power(arg, pow) }
//
//context(LabeledPolynomialSpace<C, A>)
//fun <C, A: Ring<C>> multiplyWithPower(base: C, arg: C, pow: UInt): C = ring { multiplyWithPower<C>(base, arg, pow) }
//
//// endregion
//// region Variables
//
//context(LabeledPolynomialSpace<C, A>)
//fun <C, A: Ring<C>> power(arg: Variable, pow: UInt): LabeledPolynomial<C> =
// if (pow == 0U) one
// else LabeledPolynomial<C>(mapOf(
// mapOf(arg to pow) to ring.one
// ))
//
//// endregion
//// region Polynomials
//
//context(LabeledPolynomialSpace<C, A>)
//fun <C, A: Ring<C>> number(value: Int): LabeledPolynomial<C> = ring { LabeledPolynomial<C>(mapOf(emptyMap<Variable, UInt>() to number<C>(value))) }
//
//context(LabeledPolynomialSpace<C, A>)
//fun <C, A: Ring<C>> multiplyWithPower(base: LabeledPolynomial<C>, arg: LabeledPolynomial<C>, pow: UInt): LabeledPolynomial<C> =
// 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<C, A>)
//internal tailrec fun <C, A: Ring<C>> multiplyWithPowerInternalLogic(base: LabeledPolynomial<C>, arg: LabeledPolynomial<C>, exponent: UInt): LabeledPolynomial<C> =
// 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 <C, A : Ring<C>, R> A.labeledPolynomial(block: LabeledPolynomialSpace<C, A>.() -> 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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial<C>.represent(names: Map<Variable, String> = 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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial<C>.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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial<C>.representWithBrackets(names: Map<Variable, String> = 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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial<C>.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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial<C>.representReversed(names: Map<Variable, String> = 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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial<C>.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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial<C>.representReversedWithBrackets(names: Map<Variable, String> = 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<C, A>)
//fun <C, A: Ring<C>> LabeledPolynomial<C>.representReversedWithBrackets(namer: (Variable) -> String): String =
// with(representReversed(namer)) { if (coefficients.count() == 1) this else "($this)" }
//
//// endregion
// region Operator extensions
//// region Field case
//
//operator fun <T: Field<T>> Polynomial<T>.div(other: Polynomial<T>): Polynomial<T> {
// if (other.isZero()) throw ArithmeticException("/ by zero")
// if (isZero()) return this
//
// fun Map<List<Int>, 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<List<Int>, 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 <T: Field<T>> Polynomial<T>.div(other: T): Polynomial<T> =
// if (other.isZero()) throw ArithmeticException("/ by zero")
// else
// Polynomial(
// coefficients
// .mapValues { it.value / other },
// toCheckInput = false
// )
//
//operator fun <T: Field<T>> Polynomial<T>.rem(other: Polynomial<T>): Polynomial<T> {
// if (other.isZero()) throw ArithmeticException("/ by zero")
// if (isZero()) return this
//
// fun Map<List<Int>, 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 <T: Field<T>> Polynomial<T>.divrem(other: Polynomial<T>): Polynomial.Companion.DividingResult<T> {
// if (other.isZero()) throw ArithmeticException("/ by zero")
// if (isZero()) return Polynomial.Companion.DividingResult(this, this)
//
// fun Map<List<Int>, 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<List<Int>, 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 <C> LabeledPolynomial<C>.substitute(ring: Ring<C>, args: Map<Variable, C>): LabeledPolynomial<C> = ring {
// if (coefficients.isEmpty()) return this@substitute
// LabeledPolynomial<C>(
// 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 <C> LabeledPolynomial<C>.substitute(ring: Ring<C>, arg: Map<Variable, LabeledPolynomial<C>>) : LabeledPolynomial<C> =
// 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 <C, A : Ring<C>> LabeledPolynomial<C>.asFunctionOver(ring: A): (Map<Variable, C>) -> LabeledPolynomial<C> =
// { substitute(ring, it) }
//
//fun <C, A : Ring<C>> LabeledPolynomial<C>.asPolynomialFunctionOver(ring: A): (Map<Variable, LabeledPolynomial<C>>) -> LabeledPolynomial<C> =
// { substitute(ring, it) }
//
//// endregion
//// region Algebraic derivative and antiderivative
//// TODO
//// endregion

605
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/numberedPolynomialUtil.kt Normal file
View File

@ -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<C, A>)
//@Suppress("NOTHING_TO_INLINE")
//public fun <C, A: Ring<C>> numberConstant(value: Int): C = ring { number<C>(value) }
//
//context(NumberedPolynomialSpace<C, A>)
//public fun <C, A: Ring<C>> multiplyWithPower(base: C, arg: C, pow: UInt): C = ring { multiplyWithPower<C>(base, arg, pow) }
//
//// endregion
//// region Polynomials
//
//context(NumberedPolynomialSpace<C, A>)
//public fun <C, A: Ring<C>> number(value: Int): NumberedPolynomial<C> = ring { NumberedPolynomial<C>(mapOf(emptyList<UInt>() to number<C>(value))) }
//
//context(NumberedPolynomialSpace<C, A>)
//public fun <C, A: Ring<C>> multiplyWithPower(base: NumberedPolynomial<C>, arg: NumberedPolynomial<C>, pow: UInt): NumberedPolynomial<C> =
// 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<C, A>)
//internal tailrec fun <C, A: Ring<C>> multiplyWithPowerInternalLogic(base: NumberedPolynomial<C>, arg: NumberedPolynomial<C>, exponent: UInt): NumberedPolynomial<C> =
// 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 <C, A : Ring<C>, R> A.numberedPolynomial(block: NumberedPolynomialSpace<C, A>.() -> 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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial<C>.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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial<C>.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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial<C>.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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial<C>.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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial<C>.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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial<C>.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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial<C>.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<C, A>)
//public fun <C, A: Ring<C>> NumberedPolynomial<C>.representReversedWithBrackets(namer: (Int) -> String): String =
// with(representReversed(namer)) { if (coefficients.count() == 1) this else "($this)" }
//
//// endregion
//// region Polynomial substitution and functional representation
//
//public fun <C> NumberedPolynomial<C>.substitute(ring: Ring<C>, args: Map<Int, C>): NumberedPolynomial<C> = ring {
// if (coefficients.isEmpty()) return this@substitute
// NumberedPolynomial<C>(
// 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 <C> NumberedPolynomial<C>.substitute(ring: Ring<C>, arg: Map<Int, NumberedPolynomial<C>>) : NumberedPolynomial<C> =
// 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 <C, A : Ring<C>> NumberedPolynomial<C>.asFunctionOver(ring: A): (Map<Int, C>) -> NumberedPolynomial<C> =
// { substitute(ring, it) }
//
//public fun <C, A : Ring<C>> NumberedPolynomial<C>.asPolynomialFunctionOver(ring: A): (Map<Int, NumberedPolynomial<C>>) -> NumberedPolynomial<C> =
// { substitute(ring, it) }
//
//// endregion
// region Operator extensions
//// region Field case
//
//operator fun <T: Field<T>> Polynomial<T>.div(other: Polynomial<T>): Polynomial<T> {
// if (other.isZero()) throw ArithmeticException("/ by zero")
// if (isZero()) return this
//
// fun Map<List<Int>, 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<List<Int>, 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 <T: Field<T>> Polynomial<T>.div(other: T): Polynomial<T> =
// if (other.isZero()) throw ArithmeticException("/ by zero")
// else
// Polynomial(
// coefficients
// .mapValues { it.value / other },
// toCheckInput = false
// )
//
//operator fun <T: Field<T>> Polynomial<T>.rem(other: Polynomial<T>): Polynomial<T> {
// if (other.isZero()) throw ArithmeticException("/ by zero")
// if (isZero()) return this
//
// fun Map<List<Int>, 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 <T: Field<T>> Polynomial<T>.divrem(other: Polynomial<T>): Polynomial.Companion.DividingResult<T> {
// if (other.isZero()) throw ArithmeticException("/ by zero")
// if (isZero()) return Polynomial.Companion.DividingResult(this, this)
//
// fun Map<List<Int>, 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<List<Int>, 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<Double>.substitute(args: Map<Int, Double>): NumberedPolynomial<Double> = Double.algebra {
val acc = LinkedHashMap<List<UInt>, 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<Double>(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 <C> NumberedPolynomial<C>.substitute(ring: Ring<C>, args: Map<Int, C>): NumberedPolynomial<C> = ring {
val acc = LinkedHashMap<List<UInt>, 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<C>(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 <C> NumberedPolynomial<C>.substitute(ring: Ring<C>, args: Map<Int, NumberedPolynomial<C>>) : NumberedPolynomial<C> = TODO() /*ring.numberedPolynomial {
val acc = LinkedHashMap<List<UInt>, NumberedPolynomial<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.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 <C, A : Ring<C>> NumberedPolynomial<C>.asFunction(ring: A): (Map<Int, C>) -> NumberedPolynomial<C> = { substitute(ring, it) }
/**
* Represent the polynomial as a regular context-less function.
*/
public fun <C, A : Ring<C>> NumberedPolynomial<C>.asPolynomialFunctionOver(ring: A): (Map<Int, NumberedPolynomial<C>>) -> NumberedPolynomial<C> = { substitute(ring, it) }
// endregion
// region Algebraic derivative and antiderivative
/**
* Returns algebraic derivative of received polynomial.
*/
@UnstableKMathAPI
public fun <C, A> NumberedPolynomial<C>.derivativeBy(
algebra: A,
variable: Int,
): Polynomial<C> where A : Ring<C>, A : NumericAlgebra<C> = algebra {
TODO()
}
/**
* Returns algebraic derivative of received polynomial.
*/
@UnstableKMathAPI
public fun <C, A> NumberedPolynomial<C>.derivativeBy(
algebra: A,
variables: IntArray,
): Polynomial<C> where A : Ring<C>, A : NumericAlgebra<C> = algebra {
TODO()
}
/**
* Returns algebraic derivative of received polynomial.
*/
@UnstableKMathAPI
public fun <C, A> NumberedPolynomial<C>.derivativeBy(
algebra: A,
variables: Collection<Int>,
): Polynomial<C> where A : Ring<C>, A : NumericAlgebra<C> = derivativeBy(algebra, variables.toIntArray())
/**
* Create a polynomial witch represents indefinite integral version of this polynomial
*/
@UnstableKMathAPI
public fun <C, A> NumberedPolynomial<C>.antiderivativeBy(
algebra: A,
variable: Int,
): Polynomial<C> where A : Field<C>, A : NumericAlgebra<C> = algebra {
TODO()
}
/**
* Create a polynomial witch represents indefinite integral version of this polynomial
*/
@UnstableKMathAPI
public fun <C, A> NumberedPolynomial<C>.antiderivativeBy(
algebra: A,
variables: IntArray,
): Polynomial<C> where A : Field<C>, A : NumericAlgebra<C> = algebra {
TODO()
}
/**
* Create a polynomial witch represents indefinite integral version of this polynomial
*/
@UnstableKMathAPI
public fun <C, A> NumberedPolynomial<C>.antiderivativeBy(
algebra: A,
variables: Collection<Int>,
): Polynomial<C> where A : Field<C>, A : NumericAlgebra<C> = antiderivativeBy(algebra, variables.toIntArray())
// endregion

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

@ -121,7 +121,7 @@ public fun <C, A> Polynomial<C>.antiderivative(
): Polynomial<C> where A : Field<C>, A : NumericAlgebra<C> = algebra { ): Polynomial<C> where A : Field<C>, A : NumericAlgebra<C> = algebra {
val integratedCoefficients = buildList(coefficients.size + 1) { val integratedCoefficients = buildList(coefficients.size + 1) {
add(zero) add(zero)
coefficients.forEachIndexed{ index, t -> add(t / (number(index) + one)) } coefficients.forEachIndexed{ index, t -> add(t / number(index + 1)) }
} }
Polynomial(integratedCoefficients) Polynomial(integratedCoefficients)
} }
@ -137,3 +137,5 @@ public fun <C : Comparable<C>> Polynomial<C>.integrate(
val integral = antiderivative(algebra) val integral = antiderivative(algebra)
integral.substitute(algebra, range.endInclusive) - integral.substitute(algebra, range.start) integral.substitute(algebra, range.endInclusive) - integral.substitute(algebra, range.start)
} }
// endregion