forked from kscience/kmath
Restructured Polynomial
This commit is contained in:
parent
c80f70fe0f
commit
2483c56f1c
@ -0,0 +1,354 @@
|
|||||||
|
/*
|
||||||
|
* Copyright 2018-2021 KMath contributors.
|
||||||
|
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
|
||||||
|
*/
|
||||||
|
|
||||||
|
package space.kscience.kmath.functions
|
||||||
|
|
||||||
|
import space.kscience.kmath.operations.*
|
||||||
|
import kotlin.js.JsName
|
||||||
|
import kotlin.jvm.JvmName
|
||||||
|
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Abstraction of polynomials.
|
||||||
|
*/
|
||||||
|
public interface AbstractPolynomial<C>
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Abstraction of ring of polynomials of type [P] over ring of constants of type [C].
|
||||||
|
*/
|
||||||
|
@Suppress("INAPPLICABLE_JVM_NAME", "PARAMETER_NAME_CHANGED_ON_OVERRIDE")
|
||||||
|
public interface AbstractPolynomialSpace<C, P: AbstractPolynomial<C>> : Ring<P> {
|
||||||
|
// region Constant-integer relation
|
||||||
|
/**
|
||||||
|
* Returns sum of the constant and the integer represented as constant (member of underlying ring).
|
||||||
|
*
|
||||||
|
* The operation is equivalent to adding [other] copies of unit of underlying ring to [this].
|
||||||
|
*/
|
||||||
|
@JvmName("constantIntPlus")
|
||||||
|
public operator fun C.plus(other: Int): C
|
||||||
|
/**
|
||||||
|
* Returns difference between the constant and the integer represented as constant (member of underlying ring).
|
||||||
|
*
|
||||||
|
* The operation is equivalent to subtraction [other] copies of unit of underlying ring from [this].
|
||||||
|
*/
|
||||||
|
@JvmName("constantIntMinus")
|
||||||
|
public operator fun C.minus(other: Int): C
|
||||||
|
/**
|
||||||
|
* Returns product of the constant and the integer represented as constant (member of underlying ring).
|
||||||
|
*
|
||||||
|
* The operation is equivalent to sum of [other] copies of [this].
|
||||||
|
*/
|
||||||
|
@JvmName("constantIntTimes")
|
||||||
|
public operator fun C.times(other: Int): C
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Integer-constant relation
|
||||||
|
/**
|
||||||
|
* Returns sum of the integer represented as constant (member of underlying ring) and the constant.
|
||||||
|
*
|
||||||
|
* The operation is equivalent to adding [this] copies of unit of underlying ring to [other].
|
||||||
|
*/
|
||||||
|
@JvmName("intConstantPlus")
|
||||||
|
public operator fun Int.plus(other: C): C
|
||||||
|
/**
|
||||||
|
* Returns difference between the integer represented as constant (member of underlying ring) and the constant.
|
||||||
|
*
|
||||||
|
* The operation is equivalent to subtraction [this] copies of unit of underlying ring from [other].
|
||||||
|
*/
|
||||||
|
@JvmName("intConstantMinus")
|
||||||
|
public operator fun Int.minus(other: C): C
|
||||||
|
/**
|
||||||
|
* Returns product of the integer represented as constant (member of underlying ring) and the constant.
|
||||||
|
*
|
||||||
|
* The operation is equivalent to sum of [this] copies of [other].
|
||||||
|
*/
|
||||||
|
@JvmName("intConstantTimes")
|
||||||
|
public operator fun Int.times(other: C): C
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Polynomial-integer relation
|
||||||
|
/**
|
||||||
|
* Returns sum of the constant and the integer represented as polynomial.
|
||||||
|
*
|
||||||
|
* The operation is equivalent to adding [other] copies of unit polynomial to [this].
|
||||||
|
*/
|
||||||
|
public operator fun P.plus(other: Int): P = optimizedAddMultiplied(this, one, other)
|
||||||
|
/**
|
||||||
|
* Returns difference between the constant and the integer represented as polynomial.
|
||||||
|
*
|
||||||
|
* The operation is equivalent to subtraction [other] copies of unit polynomial from [this].
|
||||||
|
*/
|
||||||
|
public operator fun P.minus(other: Int): P = optimizedAddMultiplied(this, one, -other)
|
||||||
|
/**
|
||||||
|
* Returns product of the constant and the integer represented as polynomial.
|
||||||
|
*
|
||||||
|
* The operation is equivalent to sum of [other] copies of [this].
|
||||||
|
*/
|
||||||
|
public operator fun P.times(other: Int): P = optimizedMultiply(this, other)
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Integer-polynomial relation
|
||||||
|
/**
|
||||||
|
* Returns sum of the integer represented as polynomial and the constant.
|
||||||
|
*
|
||||||
|
* The operation is equivalent to adding [this] copies of unit polynomial to [other].
|
||||||
|
*/
|
||||||
|
public operator fun Int.plus(other: P): P = optimizedAddMultiplied(other, one, this)
|
||||||
|
/**
|
||||||
|
* Returns difference between the integer represented as polynomial and the constant.
|
||||||
|
*
|
||||||
|
* The operation is equivalent to subtraction [this] copies of unit polynomial from [other].
|
||||||
|
*/
|
||||||
|
public operator fun Int.minus(other: P): P = optimizedAddMultiplied(-other, one, this)
|
||||||
|
/**
|
||||||
|
* Returns product of the integer represented as polynomial and the constant.
|
||||||
|
*
|
||||||
|
* The operation is equivalent to sum of [this] copies of [other].
|
||||||
|
*/
|
||||||
|
public operator fun Int.times(other: P): P = optimizedMultiply(other, this)
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Constant-constant relation
|
||||||
|
/**
|
||||||
|
* Returns the same constant.
|
||||||
|
*/
|
||||||
|
@JvmName("constantUnaryPlus")
|
||||||
|
@JsName("constantUnaryPlus")
|
||||||
|
public operator fun C.unaryPlus(): C = this
|
||||||
|
/**
|
||||||
|
* Returns negation of the constant.
|
||||||
|
*/
|
||||||
|
@JvmName("constantUnaryMinus")
|
||||||
|
@JsName("constantUnaryMinus")
|
||||||
|
public operator fun C.unaryMinus(): C
|
||||||
|
/**
|
||||||
|
* Returns sum of the constants.
|
||||||
|
*/
|
||||||
|
@JvmName("constantPlus")
|
||||||
|
@JsName("constantPlus")
|
||||||
|
public operator fun C.plus(other: C): C
|
||||||
|
/**
|
||||||
|
* Returns difference of the constants.
|
||||||
|
*/
|
||||||
|
@JvmName("constantMinus")
|
||||||
|
@JsName("constantMinus")
|
||||||
|
public operator fun C.minus(other: C): C
|
||||||
|
/**
|
||||||
|
* Returns product of the constants.
|
||||||
|
*/
|
||||||
|
@JvmName("constantTimes")
|
||||||
|
@JsName("constantTimes")
|
||||||
|
public operator fun C.times(other: C): C
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Check if the instant is zero constant.
|
||||||
|
*/
|
||||||
|
@JvmName("constantIsZero")
|
||||||
|
public fun C.isZero(): Boolean
|
||||||
|
/**
|
||||||
|
* Check if the instant is NOT zero constant.
|
||||||
|
*/
|
||||||
|
@JvmName("constantIsNotZero")
|
||||||
|
public fun C.isNotZero(): Boolean
|
||||||
|
/**
|
||||||
|
* Check if the instant is unit constant.
|
||||||
|
*/
|
||||||
|
@JvmName("constantIsOne")
|
||||||
|
public fun C.isOne(): Boolean
|
||||||
|
/**
|
||||||
|
* Check if the instant is NOT unit constant.
|
||||||
|
*/
|
||||||
|
@JvmName("constantIsNotOne")
|
||||||
|
public fun C.isNotOne(): Boolean
|
||||||
|
/**
|
||||||
|
* Check if the instant is minus unit constant.
|
||||||
|
*/
|
||||||
|
@JvmName("constantIsMinusOne")
|
||||||
|
public fun C.isMinusOne(): Boolean
|
||||||
|
/**
|
||||||
|
* Check if the instant is NOT minus unit constant.
|
||||||
|
*/
|
||||||
|
@JvmName("constantIsNotMinusOne")
|
||||||
|
public fun C.isNotMinusOne(): Boolean
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Constant-polynomial relation
|
||||||
|
/**
|
||||||
|
* Returns sum of the constant represented as polynomial and the polynomial.
|
||||||
|
*/
|
||||||
|
public operator fun C.plus(other: P): P
|
||||||
|
/**
|
||||||
|
* Returns difference between the constant represented as polynomial and the polynomial.
|
||||||
|
*/
|
||||||
|
public operator fun C.minus(other: P): P
|
||||||
|
/**
|
||||||
|
* Returns product of the constant represented as polynomial and the polynomial.
|
||||||
|
*/
|
||||||
|
public operator fun C.times(other: P): P
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Polynomial-constant relation
|
||||||
|
/**
|
||||||
|
* Returns sum of the constant represented as polynomial and the polynomial.
|
||||||
|
*/
|
||||||
|
public operator fun P.plus(other: C): P
|
||||||
|
/**
|
||||||
|
* Returns difference between the constant represented as polynomial and the polynomial.
|
||||||
|
*/
|
||||||
|
public operator fun P.minus(other: C): P
|
||||||
|
/**
|
||||||
|
* Returns product of the constant represented as polynomial and the polynomial.
|
||||||
|
*/
|
||||||
|
public operator fun P.times(other: C): P
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Polynomial-polynomial relation
|
||||||
|
/**
|
||||||
|
* Returns the same polynomial.
|
||||||
|
*/
|
||||||
|
public override operator fun P.unaryPlus(): P = this
|
||||||
|
/**
|
||||||
|
* Returns negation of the polynomial.
|
||||||
|
*/
|
||||||
|
public override operator fun P.unaryMinus(): P
|
||||||
|
/**
|
||||||
|
* Returns sum of the polynomials.
|
||||||
|
*/
|
||||||
|
public override operator fun P.plus(other: P): P
|
||||||
|
/**
|
||||||
|
* Returns difference of the polynomials.
|
||||||
|
*/
|
||||||
|
public override operator fun P.minus(other: P): P
|
||||||
|
/**
|
||||||
|
* Returns product of the polynomials.
|
||||||
|
*/
|
||||||
|
public override operator fun P.times(other: P): P
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Check if the instant is zero polynomial.
|
||||||
|
*/
|
||||||
|
public fun P.isZero(): Boolean = this == zero
|
||||||
|
/**
|
||||||
|
* Check if the instant is NOT zero polynomial.
|
||||||
|
*/
|
||||||
|
public fun P.isNotZero(): Boolean = this != zero
|
||||||
|
/**
|
||||||
|
* Check if the instant is unit polynomial.
|
||||||
|
*/
|
||||||
|
public fun P.isOne(): Boolean = this == one
|
||||||
|
/**
|
||||||
|
* Check if the instant is NOT unit polynomial.
|
||||||
|
*/
|
||||||
|
public fun P.isNotOne(): Boolean = this != one
|
||||||
|
/**
|
||||||
|
* Check if the instant is minus unit polynomial.
|
||||||
|
*/
|
||||||
|
public fun P.isMinusOne(): Boolean = this == -one
|
||||||
|
/**
|
||||||
|
* Check if the instant is NOT minus unit polynomial.
|
||||||
|
*/
|
||||||
|
public fun P.isNotMinusOne(): Boolean = this != -one
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Instance of zero polynomial (zero of the polynomial ring).
|
||||||
|
*/
|
||||||
|
override val zero: P
|
||||||
|
/**
|
||||||
|
* Instance of unit polynomial (unit of the polynomial ring).
|
||||||
|
*/
|
||||||
|
override val one: P
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Checks equality of the polynomials.
|
||||||
|
*/
|
||||||
|
@Suppress("EXTENSION_SHADOWED_BY_MEMBER", "CovariantEquals")
|
||||||
|
public fun P.equals(other: P): Boolean
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// Not sure is it necessary...
|
||||||
|
// region Polynomial properties
|
||||||
|
/**
|
||||||
|
* Degree of the polynomial, [see also](https://en.wikipedia.org/wiki/Degree_of_a_polynomial). If the polynomial is
|
||||||
|
* zero, degree is -1.
|
||||||
|
*/
|
||||||
|
public val P.degree: Int
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Checks if the instant is constant polynomial (of degree no more than 0) over considered ring.
|
||||||
|
*/
|
||||||
|
public fun P.isConstant(): Boolean = degree <= 0
|
||||||
|
/**
|
||||||
|
* Checks if the instant is **not** constant polynomial (of degree no more than 0) over considered ring.
|
||||||
|
*/
|
||||||
|
public fun P.isNotConstant(): Boolean = !isConstant()
|
||||||
|
/**
|
||||||
|
* Checks if the instant is constant non-zero polynomial (of degree no more than 0) over considered ring.
|
||||||
|
*/
|
||||||
|
public fun P.isNonZeroConstant(): Boolean = degree == 0
|
||||||
|
/**
|
||||||
|
* Checks if the instant is **not** constant non-zero polynomial (of degree no more than 0) over considered ring.
|
||||||
|
*/
|
||||||
|
public fun P.isNotNonZeroConstant(): Boolean = !isNonZeroConstant()
|
||||||
|
/**
|
||||||
|
* If polynomial is a constant polynomial represents and returns it as constant.
|
||||||
|
* Otherwise, (when the polynomial is not constant polynomial) returns `null`.
|
||||||
|
*/
|
||||||
|
public fun P.asConstantOrNull(): C?
|
||||||
|
/**
|
||||||
|
* If polynomial is a constant polynomial represents and returns it as constant.
|
||||||
|
* Otherwise, (when the polynomial is not constant polynomial) raises corresponding exception.
|
||||||
|
*/
|
||||||
|
public fun P.asConstant(): C = asConstantOrNull() ?: error("Can not represent non-constant polynomial as a constant")
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Legacy of Ring interface
|
||||||
|
override fun add(left: P, right: P): P = left + right
|
||||||
|
override fun multiply(left: P, right: P): P = left * right
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
public companion object {
|
||||||
|
// TODO: All of this should be moved to algebraic structures' place for utilities
|
||||||
|
// TODO: Move receiver to context receiver
|
||||||
|
/**
|
||||||
|
* Multiplication of element and integer.
|
||||||
|
*
|
||||||
|
* @receiver the multiplicand.
|
||||||
|
* @param other the multiplier.
|
||||||
|
* @return the difference.
|
||||||
|
*/
|
||||||
|
internal tailrec fun <C> Group<C>.optimizedMultiply(arg: C, other: Int): C =
|
||||||
|
when {
|
||||||
|
other == 0 -> zero
|
||||||
|
other == 1 -> arg
|
||||||
|
other == -1 -> -arg
|
||||||
|
other % 2 == 0 -> optimizedMultiply(arg + arg, other / 2)
|
||||||
|
other % 2 == 1 -> optimizedAddMultiplied(arg, arg + arg, other / 2)
|
||||||
|
other % 2 == -1 -> optimizedAddMultiplied(-arg, arg + arg, other / 2)
|
||||||
|
else -> error("Error in multiplication group instant by integer: got reminder by division by 2 different from 0, 1 and -1")
|
||||||
|
}
|
||||||
|
|
||||||
|
// TODO: Move receiver to context receiver
|
||||||
|
/**
|
||||||
|
* Adds product of [arg] and [multiplier] to [base].
|
||||||
|
*
|
||||||
|
* @receiver the algebra to provide multiplication.
|
||||||
|
* @param base the augend.
|
||||||
|
* @param arg the multiplicand.
|
||||||
|
* @param multiplier the multiplier.
|
||||||
|
* @return sum of the augend [base] and product of the multiplicand [arg] and the multiplier [multiplier].
|
||||||
|
* @author Gleb Minaev
|
||||||
|
*/
|
||||||
|
internal tailrec fun <C> Group<C>.optimizedAddMultiplied(base: C, arg: C, multiplier: Int): C =
|
||||||
|
when {
|
||||||
|
multiplier == 0 -> base
|
||||||
|
multiplier == 1 -> base + arg
|
||||||
|
multiplier == -1 -> base - arg
|
||||||
|
multiplier % 2 == 0 -> optimizedAddMultiplied(base, arg + arg, multiplier / 2)
|
||||||
|
multiplier % 2 == 1 -> optimizedAddMultiplied(base + arg, arg + arg, multiplier / 2)
|
||||||
|
multiplier % 2 == -1 -> optimizedAddMultiplied(base + arg, arg + arg, multiplier / 2)
|
||||||
|
else -> error("Error in multiplication group instant by integer: got reminder by division by 2 different from 0, 1 and -1")
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
@ -117,16 +117,16 @@ public fun <T : Comparable<T>> PiecewisePolynomial(
|
|||||||
* Return a value of polynomial function with given [ring] a given [arg] or null if argument is outside piecewise
|
* Return a value of polynomial function with given [ring] a given [arg] or null if argument is outside piecewise
|
||||||
* definition.
|
* definition.
|
||||||
*/
|
*/
|
||||||
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.value(ring: C, arg: T): T? =
|
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.substitute(ring: C, arg: T): T? =
|
||||||
findPiece(arg)?.value(ring, arg)
|
findPiece(arg)?.substitute(ring, arg)
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Convert this polynomial to a function returning nullable value (null if argument is outside piecewise range).
|
* Convert this polynomial to a function returning nullable value (null if argument is outside piecewise range).
|
||||||
*/
|
*/
|
||||||
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { value(ring, it) }
|
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { substitute(ring, it) }
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Convert this polynomial to a function using [defaultValue] for arguments outside the piecewise range.
|
* Convert this polynomial to a function using [defaultValue] for arguments outside the piecewise range.
|
||||||
*/
|
*/
|
||||||
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C, defaultValue: T): (T) -> T =
|
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C, defaultValue: T): (T) -> T =
|
||||||
{ value(ring, it) ?: defaultValue }
|
{ substitute(ring, it) ?: defaultValue }
|
||||||
|
@ -5,128 +5,371 @@
|
|||||||
|
|
||||||
package space.kscience.kmath.functions
|
package space.kscience.kmath.functions
|
||||||
|
|
||||||
import space.kscience.kmath.misc.UnstableKMathAPI
|
import space.kscience.kmath.functions.AbstractPolynomialSpace.Companion.optimizedAddMultiplied
|
||||||
|
import space.kscience.kmath.functions.AbstractPolynomialSpace.Companion.optimizedMultiply
|
||||||
import space.kscience.kmath.operations.*
|
import space.kscience.kmath.operations.*
|
||||||
import kotlin.contracts.InvocationKind
|
import kotlin.jvm.JvmName
|
||||||
import kotlin.contracts.contract
|
|
||||||
import kotlin.math.max
|
import kotlin.math.max
|
||||||
import kotlin.math.pow
|
import kotlin.math.min
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Polynomial coefficients model without fixation on specific context they are applied to.
|
* Polynomial coefficients 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<out T>(public val coefficients: List<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>
|
||||||
|
// )
|
||||||
|
// }
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Returns a [Polynomial] instance with given [coefficients].
|
* Represents internal [Polynomial] errors.
|
||||||
|
*/
|
||||||
|
internal class PolynomialError(message: String): Error(message)
|
||||||
|
|
||||||
|
// region Constructors and converters
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Returns a [Polynomial] instance with given [coefficients]. The collection of coefficients will be reversed if
|
||||||
|
* [reverse] parameter is true.
|
||||||
*/
|
*/
|
||||||
@Suppress("FunctionName")
|
@Suppress("FunctionName")
|
||||||
public fun <T> Polynomial(vararg coefficients: T): Polynomial<T> = Polynomial(coefficients.toList())
|
public fun <T> Polynomial(coefficients: List<T>, reverse: Boolean = false): Polynomial<T> =
|
||||||
|
Polynomial(with(coefficients) { if (reverse) reversed() else this })
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Evaluates the value of the given double polynomial for given double argument.
|
* Returns a [Polynomial] instance with given [coefficients]. The collection of coefficients will be reversed if
|
||||||
|
* [reverse] parameter is true.
|
||||||
*/
|
*/
|
||||||
public fun Polynomial<Double>.value(arg: Double): Double = coefficients.reduceIndexed { index, acc, c ->
|
@Suppress("FunctionName")
|
||||||
acc + c * arg.pow(index)
|
public fun <T> Polynomial(vararg coefficients: T, reverse: Boolean = false): Polynomial<T> =
|
||||||
}
|
Polynomial(with(coefficients) { if (reverse) reversed() else toList() })
|
||||||
|
|
||||||
|
public fun <T> T.asPolynomial() : Polynomial<T> = Polynomial(listOf(this))
|
||||||
|
|
||||||
|
// endregion
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Evaluates the value of the given polynomial for given argument.
|
* Space of polynomials constructed over ring.
|
||||||
* https://en.wikipedia.org/wiki/Horner%27s_method
|
|
||||||
*/
|
|
||||||
public fun <T, C : Ring<T>> Polynomial<T>.value(ring: C, arg: T): T = ring {
|
|
||||||
if (coefficients.isEmpty()) return@ring zero
|
|
||||||
var result: T = coefficients.last()
|
|
||||||
for (j in coefficients.size - 2 downTo 0) {
|
|
||||||
result = (arg * result) + coefficients[j]
|
|
||||||
}
|
|
||||||
return result
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Represent the polynomial as a regular context-less function.
|
|
||||||
*/
|
|
||||||
public fun <T, C : Ring<T>> Polynomial<T>.asFunction(ring: C): (T) -> T = { value(ring, it) }
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Create a polynomial witch represents differentiated version of this polynomial
|
|
||||||
*/
|
|
||||||
@UnstableKMathAPI
|
|
||||||
public fun <T, A> Polynomial<T>.differentiate(
|
|
||||||
algebra: A,
|
|
||||||
): Polynomial<T> where A : Ring<T>, A : NumericAlgebra<T> = algebra {
|
|
||||||
Polynomial(coefficients.drop(1).mapIndexed { index, t -> number(index) * t })
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Create a polynomial witch represents indefinite integral version of this polynomial
|
|
||||||
*/
|
|
||||||
@UnstableKMathAPI
|
|
||||||
public fun <T, A> Polynomial<T>.integrate(
|
|
||||||
algebra: A,
|
|
||||||
): Polynomial<T> where A : Field<T>, A : NumericAlgebra<T> = algebra {
|
|
||||||
val integratedCoefficients = buildList(coefficients.size + 1) {
|
|
||||||
add(zero)
|
|
||||||
coefficients.forEachIndexed{ index, t -> add(t / (number(index) + one)) }
|
|
||||||
}
|
|
||||||
Polynomial(integratedCoefficients)
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Compute a definite integral of a given polynomial in a [range]
|
|
||||||
*/
|
|
||||||
@UnstableKMathAPI
|
|
||||||
public fun <T : Comparable<T>> Polynomial<T>.integrate(
|
|
||||||
algebra: Field<T>,
|
|
||||||
range: ClosedRange<T>,
|
|
||||||
): T = algebra {
|
|
||||||
val integral = integrate(algebra)
|
|
||||||
integral.value(algebra, range.endInclusive) - integral.value(algebra, range.start)
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Space of polynomials.
|
|
||||||
*
|
*
|
||||||
* @param T the type of operated polynomials.
|
* @param C the type of constants. Polynomials have them as a coefficients in their terms.
|
||||||
* @param C the intersection of [Ring] of [T] and [ScaleOperations] of [T].
|
* @param A type of underlying ring of constants. It's [Ring] of [C].
|
||||||
* @param ring the [C] instance.
|
* @param ring underlying ring of constants of type [A].
|
||||||
*/
|
*/
|
||||||
public class PolynomialSpace<T, C>(
|
@Suppress("INAPPLICABLE_JVM_NAME") // KT-31420
|
||||||
private val ring: C,
|
public open class PolynomialSpace<C, A : Ring<C>>(
|
||||||
) : Group<Polynomial<T>>, ScaleOperations<Polynomial<T>> where C : Ring<T>, C : ScaleOperations<T> {
|
public val ring: A,
|
||||||
override val zero: Polynomial<T> = Polynomial(emptyList())
|
) : AbstractPolynomialSpace<C, Polynomial<C>>{
|
||||||
|
// region Constant-integer relation
|
||||||
|
@JvmName("constantIntPlus")
|
||||||
|
public override operator fun C.plus(other: Int): C = ring { optimizedAddMultiplied(this@plus, one, other) }
|
||||||
|
@JvmName("constantIntMinus")
|
||||||
|
public override operator fun C.minus(other: Int): C = ring { optimizedAddMultiplied(this@minus, one, -other) }
|
||||||
|
@JvmName("constantIntTimes")
|
||||||
|
public override operator fun C.times(other: Int): C = ring { optimizedMultiply(this@times, other) }
|
||||||
|
// endregion
|
||||||
|
|
||||||
override fun Polynomial<T>.unaryMinus(): Polynomial<T> = ring {
|
// region Integer-constant relation
|
||||||
|
@JvmName("intConstantPlus")
|
||||||
|
public override operator fun Int.plus(other: C): C = ring { optimizedAddMultiplied(other, one, this@plus) }
|
||||||
|
@JvmName("intConstantMinus")
|
||||||
|
public override operator fun Int.minus(other: C): C = ring { optimizedAddMultiplied(-other, one, this@minus) }
|
||||||
|
@JvmName("intConstantTimes")
|
||||||
|
public override operator fun Int.times(other: C): C = ring { optimizedMultiply(other, this@times) }
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Polynomial-integer relation
|
||||||
|
public override operator fun Polynomial<C>.plus(other: Int): Polynomial<C> =
|
||||||
|
if (other == 0) this
|
||||||
|
else
|
||||||
|
Polynomial(
|
||||||
|
coefficients
|
||||||
|
.toMutableList()
|
||||||
|
.apply {
|
||||||
|
if (isEmpty()) this[0] = ring.zero + other
|
||||||
|
else this[0] = this[0]!! + other
|
||||||
|
}
|
||||||
|
)
|
||||||
|
public override operator fun Polynomial<C>.minus(other: Int): Polynomial<C> =
|
||||||
|
if (other == 0) this
|
||||||
|
else
|
||||||
|
Polynomial(
|
||||||
|
coefficients
|
||||||
|
.toMutableList()
|
||||||
|
.apply {
|
||||||
|
if (isEmpty()) this[0] = ring.zero - other
|
||||||
|
else this[0] = this[0]!! - other
|
||||||
|
}
|
||||||
|
)
|
||||||
|
public override operator fun Polynomial<C>.times(other: Int): Polynomial<C> =
|
||||||
|
if (other == 0) zero
|
||||||
|
else Polynomial(
|
||||||
|
coefficients
|
||||||
|
.subList(0, degree + 1)
|
||||||
|
.map { it * other }
|
||||||
|
)
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Integer-polynomial relation
|
||||||
|
public override operator fun Int.plus(other: Polynomial<C>): Polynomial<C> =
|
||||||
|
if (this == 0) other
|
||||||
|
else
|
||||||
|
Polynomial(
|
||||||
|
other.coefficients
|
||||||
|
.toMutableList()
|
||||||
|
.apply {
|
||||||
|
if (isEmpty()) this[0] = ring.zero + this@plus
|
||||||
|
else this[0] = this[0]!! + this@plus
|
||||||
|
}
|
||||||
|
)
|
||||||
|
public override operator fun Int.minus(other: Polynomial<C>): Polynomial<C> =
|
||||||
|
if (this == 0) other
|
||||||
|
else
|
||||||
|
Polynomial(
|
||||||
|
other.coefficients
|
||||||
|
.toMutableList()
|
||||||
|
.apply {
|
||||||
|
if (isEmpty()) this[0] = ring.zero - this@minus
|
||||||
|
else this[0] = this[0]!! - this@minus
|
||||||
|
}
|
||||||
|
)
|
||||||
|
public override operator fun Int.times(other: Polynomial<C>): Polynomial<C> =
|
||||||
|
if (this == 0) zero
|
||||||
|
else Polynomial(
|
||||||
|
other.coefficients
|
||||||
|
.subList(0, other.degree + 1)
|
||||||
|
.map { it * this }
|
||||||
|
)
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Constant-constant relation
|
||||||
|
@JvmName("constantUnaryMinus")
|
||||||
|
public override operator fun C.unaryMinus(): C = ring { -this@unaryMinus }
|
||||||
|
@JvmName("constantPlus")
|
||||||
|
public override operator fun C.plus(other: C): C = ring { this@plus + other }
|
||||||
|
@JvmName("constantMinus")
|
||||||
|
public override operator fun C.minus(other: C): C = ring { this@minus - other }
|
||||||
|
@JvmName("constantTimes")
|
||||||
|
public override operator fun C.times(other: C): C = ring { this@times * other }
|
||||||
|
@JvmName("constantIsZero")
|
||||||
|
public override fun C.isZero(): Boolean = ring { this == zero }
|
||||||
|
@JvmName("constantIsNotZero")
|
||||||
|
public override fun C.isNotZero(): Boolean = ring { this != zero }
|
||||||
|
@JvmName("constantIsOne")
|
||||||
|
public override fun C.isOne(): Boolean = ring { this == one }
|
||||||
|
@JvmName("constantIsNotOne")
|
||||||
|
public override fun C.isNotOne(): Boolean = ring { this != one }
|
||||||
|
@JvmName("constantIsMinusOne")
|
||||||
|
public override fun C.isMinusOne(): Boolean = ring { this == -one }
|
||||||
|
@JvmName("constantIsNotMinusOne")
|
||||||
|
public override fun C.isNotMinusOne(): Boolean = ring { this != -one }
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Constant-polynomial relation
|
||||||
|
public override operator fun C.plus(other: Polynomial<C>): Polynomial<C> =
|
||||||
|
with(other.coefficients) {
|
||||||
|
if (isEmpty()) Polynomial(listOf(this@plus))
|
||||||
|
else Polynomial(
|
||||||
|
toMutableList()
|
||||||
|
.apply {
|
||||||
|
this[0] += this@plus
|
||||||
|
}
|
||||||
|
)
|
||||||
|
}
|
||||||
|
// if (degree == -1) UnivariatePolynomial(other) else UnivariatePolynomial(
|
||||||
|
// listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1)
|
||||||
|
// )
|
||||||
|
public override operator fun C.minus(other: Polynomial<C>): Polynomial<C> =
|
||||||
|
with(other.coefficients) {
|
||||||
|
if (isEmpty()) Polynomial(listOf(-this@minus))
|
||||||
|
else Polynomial(
|
||||||
|
toMutableList()
|
||||||
|
.apply {
|
||||||
|
this[0] -= this@minus
|
||||||
|
}
|
||||||
|
)
|
||||||
|
}
|
||||||
|
// if (degree == -1) UnivariatePolynomial(other) else UnivariatePolynomial(
|
||||||
|
// listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1)
|
||||||
|
// )
|
||||||
|
public override operator fun C.times(other: Polynomial<C>): Polynomial<C> =
|
||||||
|
Polynomial(
|
||||||
|
other.coefficients
|
||||||
|
// .subList(0, other.degree + 1)
|
||||||
|
.map { it * this }
|
||||||
|
)
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Polynomial-constant relation
|
||||||
|
public override operator fun Polynomial<C>.plus(other: C): Polynomial<C> =
|
||||||
|
if (other.isZero()) this
|
||||||
|
else with(coefficients) {
|
||||||
|
if (isEmpty()) Polynomial(listOf(other))
|
||||||
|
else Polynomial(
|
||||||
|
toMutableList()
|
||||||
|
.apply {
|
||||||
|
this[0] += other
|
||||||
|
}
|
||||||
|
)
|
||||||
|
}
|
||||||
|
// if (degree == -1) UnivariatePolynomial(other) else UnivariatePolynomial(
|
||||||
|
// listOf(coefficients[0] + other) + coefficients.subList(1, degree + 1)
|
||||||
|
// )
|
||||||
|
public override operator fun Polynomial<C>.minus(other: C): Polynomial<C> =
|
||||||
|
with(coefficients) {
|
||||||
|
if (isEmpty()) Polynomial(listOf(other))
|
||||||
|
else Polynomial(
|
||||||
|
toMutableList()
|
||||||
|
.apply {
|
||||||
|
this[0] -= other
|
||||||
|
}
|
||||||
|
)
|
||||||
|
}
|
||||||
|
// if (degree == -1) UnivariatePolynomial(-other) else UnivariatePolynomial(
|
||||||
|
// listOf(coefficients[0] - other) + coefficients.subList(1, degree + 1)
|
||||||
|
// )
|
||||||
|
public override operator fun Polynomial<C>.times(other: C): Polynomial<C> =
|
||||||
|
Polynomial(
|
||||||
|
coefficients
|
||||||
|
// .subList(0, degree + 1)
|
||||||
|
.map { it * other }
|
||||||
|
)
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Polynomial-polynomial relation
|
||||||
|
public override operator fun Polynomial<C>.unaryMinus(): Polynomial<C> =
|
||||||
Polynomial(coefficients.map { -it })
|
Polynomial(coefficients.map { -it })
|
||||||
|
public override operator fun Polynomial<C>.plus(other: Polynomial<C>): Polynomial<C> =
|
||||||
|
Polynomial(
|
||||||
|
(0..max(degree, other.degree))
|
||||||
|
.map {
|
||||||
|
when {
|
||||||
|
it > degree -> other.coefficients[it]
|
||||||
|
it > other.degree -> coefficients[it]
|
||||||
|
else -> coefficients[it] + other.coefficients[it]
|
||||||
}
|
}
|
||||||
|
}
|
||||||
override fun add(left: Polynomial<T>, right: Polynomial<T>): Polynomial<T> {
|
.ifEmpty { listOf(ring.zero) }
|
||||||
val dim = max(left.coefficients.size, right.coefficients.size)
|
)
|
||||||
|
public override operator fun Polynomial<C>.minus(other: Polynomial<C>): Polynomial<C> =
|
||||||
return ring {
|
Polynomial(
|
||||||
Polynomial(List(dim) { index ->
|
(0..max(degree, other.degree))
|
||||||
left.coefficients.getOrElse(index) { zero } + right.coefficients.getOrElse(index) { zero }
|
.map {
|
||||||
})
|
when {
|
||||||
|
it > degree -> -other.coefficients[it]
|
||||||
|
it > other.degree -> coefficients[it]
|
||||||
|
else -> coefficients[it] - other.coefficients[it]
|
||||||
|
}
|
||||||
|
}
|
||||||
|
.ifEmpty { listOf(ring.zero) }
|
||||||
|
)
|
||||||
|
public override operator fun Polynomial<C>.times(other: Polynomial<C>): Polynomial<C> {
|
||||||
|
val thisDegree = degree
|
||||||
|
val otherDegree = other.degree
|
||||||
|
return when {
|
||||||
|
thisDegree == -1 -> this
|
||||||
|
otherDegree == -1 -> other
|
||||||
|
else ->
|
||||||
|
Polynomial(
|
||||||
|
(0..(thisDegree + otherDegree))
|
||||||
|
.map { d ->
|
||||||
|
(max(0, d - otherDegree)..(min(thisDegree, d)))
|
||||||
|
.map { coefficients[it] * other.coefficients[d - it] }
|
||||||
|
.reduce { acc, rational -> acc + rational }
|
||||||
|
}
|
||||||
|
)
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
override fun scale(a: Polynomial<T>, value: Double): Polynomial<T> =
|
public override fun Polynomial<C>.isZero(): Boolean = coefficients.all { it.isZero() }
|
||||||
ring { Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * value }) }
|
public override fun Polynomial<C>.isNotZero(): Boolean = coefficients.any { it.isNotZero() }
|
||||||
|
public override fun Polynomial<C>.isOne(): Boolean =
|
||||||
|
with(coefficients) { isNotEmpty() && asSequence().withIndex().any { (index, c) -> if (index == 0) c.isOne() else c.isZero() } } // TODO: It's better to write new methods like `anyIndexed`. But what's better way to do it?
|
||||||
|
public override fun Polynomial<C>.isNotOne(): Boolean = !isOne()
|
||||||
|
public override fun Polynomial<C>.isMinusOne(): Boolean =
|
||||||
|
with(coefficients) { isNotEmpty() && asSequence().withIndex().any { (index, c) -> if (index == 0) c.isMinusOne() else c.isZero() } } // TODO: It's better to write new methods like `anyIndexed`. But what's better way to do it?
|
||||||
|
public override fun Polynomial<C>.isNotMinusOne(): Boolean = !isMinusOne()
|
||||||
|
|
||||||
|
override val zero: Polynomial<C> = Polynomial(emptyList())
|
||||||
|
override val one: Polynomial<C> = Polynomial(listOf(ring.one))
|
||||||
|
|
||||||
|
@Suppress("EXTENSION_SHADOWED_BY_MEMBER", "CovariantEquals")
|
||||||
|
public override fun Polynomial<C>.equals(other: Polynomial<C>): Boolean =
|
||||||
|
when {
|
||||||
|
this === other -> true
|
||||||
|
else -> {
|
||||||
|
if (this.degree == other.degree)
|
||||||
|
(0..degree).all { coefficients[it] == other.coefficients[it] }
|
||||||
|
else false
|
||||||
|
}
|
||||||
|
}
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// Not sure is it necessary...
|
||||||
|
// region Polynomial properties
|
||||||
|
public override val Polynomial<C>.degree: Int get() = coefficients.indexOfLast { it != ring.zero }
|
||||||
|
|
||||||
|
public override fun Polynomial<C>.asConstantOrNull(): C? =
|
||||||
|
with(coefficients) {
|
||||||
|
when {
|
||||||
|
isEmpty() -> ring.zero
|
||||||
|
degree > 0 -> null
|
||||||
|
else -> first()
|
||||||
|
}
|
||||||
|
}
|
||||||
|
public override fun Polynomial<C>.asConstant(): C = asConstantOrNull() ?: error("Can not represent non-constant polynomial as a constant")
|
||||||
|
|
||||||
|
@Suppress("NOTHING_TO_INLINE")
|
||||||
|
public inline fun Polynomial<C>.substitute(argument: C): C = this.substitute(ring, argument)
|
||||||
|
@Suppress("NOTHING_TO_INLINE")
|
||||||
|
public inline fun Polynomial<C>.substitute(argument: Polynomial<C>): Polynomial<C> = this.substitute(ring, argument)
|
||||||
|
|
||||||
|
@Suppress("NOTHING_TO_INLINE")
|
||||||
|
public inline fun Polynomial<C>.asFunction(): (C) -> C = { this.substitute(ring, it) }
|
||||||
|
@Suppress("NOTHING_TO_INLINE")
|
||||||
|
public inline fun Polynomial<C>.asFunctionOnConstants(): (C) -> C = { this.substitute(ring, it) }
|
||||||
|
@Suppress("NOTHING_TO_INLINE")
|
||||||
|
public inline fun Polynomial<C>.asFunctionOnPolynomials(): (Polynomial<C>) -> Polynomial<C> = { this.substitute(ring, it) }
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Evaluates the polynomial for the given value [arg].
|
* Evaluates the polynomial for the given value [arg].
|
||||||
*/
|
*/
|
||||||
public operator fun Polynomial<T>.invoke(arg: T): T = value(ring, arg)
|
@Suppress("NOTHING_TO_INLINE")
|
||||||
|
public inline operator fun Polynomial<C>.invoke(argument: C): C = this.substitute(ring, argument)
|
||||||
public fun Polynomial<T>.asFunction(): (T) -> T = asFunction(ring)
|
@Suppress("NOTHING_TO_INLINE")
|
||||||
|
public inline operator fun Polynomial<C>.invoke(argument: Polynomial<C>): Polynomial<C> = this.substitute(ring, argument)
|
||||||
|
// endregion
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
public inline fun <T, C, R> C.polynomial(block: PolynomialSpace<T, C>.() -> R): R where C : Ring<T>, C : ScaleOperations<T> {
|
/**
|
||||||
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
|
* Space of polynomials constructed over ring.
|
||||||
return PolynomialSpace(this).block()
|
*
|
||||||
|
* @param C the type of constants. Polynomials have them as a coefficients in their terms.
|
||||||
|
* @param A type of underlying ring of constants. It's [Ring] of [C].
|
||||||
|
* @param ring underlying ring of constants of type [A].
|
||||||
|
*/
|
||||||
|
public class ScalablePolynomialSpace<C, A>(
|
||||||
|
ring: A,
|
||||||
|
) : PolynomialSpace<C, A>(ring), ScaleOperations<Polynomial<C>> where A : Ring<C>, A : ScaleOperations<C> {
|
||||||
|
|
||||||
|
override fun scale(a: Polynomial<C>, value: Double): Polynomial<C> =
|
||||||
|
ring { Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * value }) }
|
||||||
|
|
||||||
}
|
}
|
||||||
|
@ -0,0 +1,139 @@
|
|||||||
|
/*
|
||||||
|
* Copyright 2018-2021 KMath contributors.
|
||||||
|
* Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
|
||||||
|
*/
|
||||||
|
|
||||||
|
package space.kscience.kmath.functions
|
||||||
|
|
||||||
|
import space.kscience.kmath.misc.UnstableKMathAPI
|
||||||
|
import space.kscience.kmath.operations.*
|
||||||
|
import kotlin.contracts.InvocationKind
|
||||||
|
import kotlin.contracts.contract
|
||||||
|
import kotlin.jvm.JvmName
|
||||||
|
import kotlin.math.pow
|
||||||
|
|
||||||
|
|
||||||
|
// region Utilities
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Removes zeros on the end of the coefficient list of polynomial.
|
||||||
|
*/
|
||||||
|
//context(PolynomialSpace<C, A>)
|
||||||
|
//fun <C, A: Ring<C>> Polynomial<C>.removeZeros() : Polynomial<C> =
|
||||||
|
// if (degree > -1) Polynomial(coefficients.subList(0, degree + 1)) else zero
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Crates a [PolynomialSpace] over received ring.
|
||||||
|
*/
|
||||||
|
public fun <C, A : Ring<C>> A.polynomial(): PolynomialSpace<C, A> =
|
||||||
|
PolynomialSpace(this)
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Crates a [PolynomialSpace]'s scope over received ring.
|
||||||
|
*/
|
||||||
|
public inline fun <C, A : Ring<C>, R> A.polynomial(block: PolynomialSpace<C, A>.() -> R): R {
|
||||||
|
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
|
||||||
|
return PolynomialSpace(this).block()
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Crates a [ScalablePolynomialSpace] over received scalable ring.
|
||||||
|
*/
|
||||||
|
public fun <C, A> A.scalablePolynomial(): ScalablePolynomialSpace<C, A> where A : Ring<C>, A : ScaleOperations<C> =
|
||||||
|
ScalablePolynomialSpace(this)
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Crates a [ScalablePolynomialSpace]'s scope over received scalable ring.
|
||||||
|
*/
|
||||||
|
public inline fun <C, A, R> A.scalablePolynomial(block: ScalablePolynomialSpace<C, A>.() -> R): R where A : Ring<C>, A : ScaleOperations<C> {
|
||||||
|
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
|
||||||
|
return ScalablePolynomialSpace(this).block()
|
||||||
|
}
|
||||||
|
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Polynomial substitution and functional representation
|
||||||
|
|
||||||
|
// TODO: May be apply Horner's method too?
|
||||||
|
/**
|
||||||
|
* Evaluates the value of the given double polynomial for given double argument.
|
||||||
|
*/
|
||||||
|
public fun Polynomial<Double>.substitute(arg: Double): Double =
|
||||||
|
coefficients.reduceIndexedOrNull { index, acc, c ->
|
||||||
|
acc + c * arg.pow(index)
|
||||||
|
} ?: .0
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Evaluates the value of the given polynomial for given argument.
|
||||||
|
*
|
||||||
|
* It is an implementation of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
|
||||||
|
*/
|
||||||
|
public fun <C> Polynomial<C>.substitute(ring: Ring<C>, arg: C): C = ring {
|
||||||
|
if (coefficients.isEmpty()) return@ring zero
|
||||||
|
var result: C = coefficients.last()
|
||||||
|
for (j in coefficients.size - 2 downTo 0) {
|
||||||
|
result = (arg * result) + coefficients[j]
|
||||||
|
}
|
||||||
|
return result
|
||||||
|
}
|
||||||
|
|
||||||
|
// TODO: (Waiting for hero) Replace with optimisation: the [result] may be unboxed, and all operations may be performed
|
||||||
|
// as soon as possible on it
|
||||||
|
public fun <C> Polynomial<C>.substitute(ring: Ring<C>, arg: Polynomial<C>) : Polynomial<C> = ring.polynomial {
|
||||||
|
if (coefficients.isEmpty()) return zero
|
||||||
|
var result: Polynomial<C> = coefficients.last().asPolynomial()
|
||||||
|
for (j in coefficients.size - 2 downTo 0) {
|
||||||
|
result = (arg * result) + coefficients[j]
|
||||||
|
}
|
||||||
|
return result
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Represent the polynomial as a regular context-less function.
|
||||||
|
*/
|
||||||
|
public fun <C, A : Ring<C>> Polynomial<C>.asFunction(ring: A): (C) -> C = { substitute(ring, it) }
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Represent the polynomial as a regular context-less function.
|
||||||
|
*/
|
||||||
|
public fun <C, A : Ring<C>> Polynomial<C>.asPolynomialFunctionOver(ring: A): (Polynomial<C>) -> Polynomial<C> = { substitute(ring, it) }
|
||||||
|
|
||||||
|
// endregion
|
||||||
|
|
||||||
|
// region Algebraic derivative and antiderivative
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Returns algebraic derivative of received polynomial.
|
||||||
|
*/
|
||||||
|
@UnstableKMathAPI
|
||||||
|
public fun <C, A> Polynomial<C>.derivative(
|
||||||
|
algebra: A,
|
||||||
|
): Polynomial<C> where A : Ring<C>, A : NumericAlgebra<C> = algebra {
|
||||||
|
Polynomial(coefficients.drop(1).mapIndexed { index, t -> number(index) * t })
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Create a polynomial witch represents indefinite integral version of this polynomial
|
||||||
|
*/
|
||||||
|
@UnstableKMathAPI
|
||||||
|
public fun <C, A> Polynomial<C>.antiderivative(
|
||||||
|
algebra: A,
|
||||||
|
): Polynomial<C> where A : Field<C>, A : NumericAlgebra<C> = algebra {
|
||||||
|
val integratedCoefficients = buildList(coefficients.size + 1) {
|
||||||
|
add(zero)
|
||||||
|
coefficients.forEachIndexed{ index, t -> add(t / (number(index) + one)) }
|
||||||
|
}
|
||||||
|
Polynomial(integratedCoefficients)
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Compute a definite integral of a given polynomial in a [range]
|
||||||
|
*/
|
||||||
|
@UnstableKMathAPI
|
||||||
|
public fun <C : Comparable<C>> Polynomial<C>.integrate(
|
||||||
|
algebra: Field<C>,
|
||||||
|
range: ClosedRange<C>,
|
||||||
|
): C = algebra {
|
||||||
|
val integral = antiderivative(algebra)
|
||||||
|
integral.substitute(algebra, range.endInclusive) - integral.substitute(algebra, range.start)
|
||||||
|
}
|
@ -7,6 +7,7 @@ package space.kscience.kmath.integration
|
|||||||
|
|
||||||
import space.kscience.kmath.functions.PiecewisePolynomial
|
import space.kscience.kmath.functions.PiecewisePolynomial
|
||||||
import space.kscience.kmath.functions.integrate
|
import space.kscience.kmath.functions.integrate
|
||||||
|
import space.kscience.kmath.functions.antiderivative
|
||||||
import space.kscience.kmath.interpolation.PolynomialInterpolator
|
import space.kscience.kmath.interpolation.PolynomialInterpolator
|
||||||
import space.kscience.kmath.interpolation.SplineInterpolator
|
import space.kscience.kmath.interpolation.SplineInterpolator
|
||||||
import space.kscience.kmath.interpolation.interpolatePolynomials
|
import space.kscience.kmath.interpolation.interpolatePolynomials
|
||||||
@ -23,7 +24,7 @@ import space.kscience.kmath.structures.MutableBufferFactory
|
|||||||
@OptIn(PerformancePitfall::class)
|
@OptIn(PerformancePitfall::class)
|
||||||
@UnstableKMathAPI
|
@UnstableKMathAPI
|
||||||
public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(algebra: Field<T>): PiecewisePolynomial<T> =
|
public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(algebra: Field<T>): PiecewisePolynomial<T> =
|
||||||
PiecewisePolynomial(pieces.map { it.first to it.second.integrate(algebra) })
|
PiecewisePolynomial(pieces.map { it.first to it.second.antiderivative(algebra) })
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Compute definite integral of given [PiecewisePolynomial] piece by piece in a given [range]
|
* Compute definite integral of given [PiecewisePolynomial] piece by piece in a given [range]
|
||||||
|
@ -9,7 +9,7 @@ package space.kscience.kmath.interpolation
|
|||||||
|
|
||||||
import space.kscience.kmath.data.XYColumnarData
|
import space.kscience.kmath.data.XYColumnarData
|
||||||
import space.kscience.kmath.functions.PiecewisePolynomial
|
import space.kscience.kmath.functions.PiecewisePolynomial
|
||||||
import space.kscience.kmath.functions.value
|
import space.kscience.kmath.functions.substitute
|
||||||
import space.kscience.kmath.misc.UnstableKMathAPI
|
import space.kscience.kmath.misc.UnstableKMathAPI
|
||||||
import space.kscience.kmath.operations.Ring
|
import space.kscience.kmath.operations.Ring
|
||||||
import space.kscience.kmath.structures.Buffer
|
import space.kscience.kmath.structures.Buffer
|
||||||
@ -33,7 +33,7 @@ public interface PolynomialInterpolator<T : Comparable<T>> : Interpolator<T, T,
|
|||||||
public fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T>
|
public fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T>
|
||||||
|
|
||||||
override fun interpolate(points: XYColumnarData<T, T, T>): (T) -> T = { x ->
|
override fun interpolate(points: XYColumnarData<T, T, T>): (T) -> T = { x ->
|
||||||
interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue()
|
interpolatePolynomials(points).substitute(algebra, x) ?: getDefaultValue()
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -5,13 +5,22 @@
|
|||||||
|
|
||||||
package space.kscience.kmath.functions
|
package space.kscience.kmath.functions
|
||||||
|
|
||||||
|
import space.kscience.kmath.operations.algebra
|
||||||
import kotlin.test.Test
|
import kotlin.test.Test
|
||||||
import kotlin.test.assertEquals
|
import kotlin.test.assertEquals
|
||||||
|
|
||||||
class PolynomialTest {
|
class PolynomialTest {
|
||||||
|
@Test
|
||||||
|
fun simple_polynomial_test() {
|
||||||
|
Double.algebra.scalablePolynomial {
|
||||||
|
val x = Polynomial(listOf(0.0, 1.0))
|
||||||
|
val polynomial = x * x - 2 * x + 1
|
||||||
|
assertEquals(0.0, polynomial.substitute(1.0), 0.001)
|
||||||
|
}
|
||||||
|
}
|
||||||
@Test
|
@Test
|
||||||
fun testIntegration() {
|
fun testIntegration() {
|
||||||
val polynomial = Polynomial(1.0, -2.0, 1.0)
|
val polynomial = Polynomial(1.0, -2.0, 1.0)
|
||||||
assertEquals(0.0, polynomial.value(1.0), 0.001)
|
assertEquals(0.0, polynomial.substitute(1.0), 0.001)
|
||||||
}
|
}
|
||||||
}
|
}
|
@ -274,11 +274,6 @@ argparse@^2.0.1:
|
|||||||
resolved "https://registry.yarnpkg.com/argparse/-/argparse-2.0.1.tgz#246f50f3ca78a3240f6c997e8a9bd1eac49e4b38"
|
resolved "https://registry.yarnpkg.com/argparse/-/argparse-2.0.1.tgz#246f50f3ca78a3240f6c997e8a9bd1eac49e4b38"
|
||||||
integrity sha512-8+9WqebbFzpX9OR+Wa6O29asIogeRMzcGtAINdpMHHyAg10f05aSFVBbcEqGf/PXw1EjAZ+q2/bEBg3DvurK3Q==
|
integrity sha512-8+9WqebbFzpX9OR+Wa6O29asIogeRMzcGtAINdpMHHyAg10f05aSFVBbcEqGf/PXw1EjAZ+q2/bEBg3DvurK3Q==
|
||||||
|
|
||||||
astring@1.7.5:
|
|
||||||
version "1.7.5"
|
|
||||||
resolved "https://registry.yarnpkg.com/astring/-/astring-1.7.5.tgz#a7d47fceaf32b052d33a3d07c511efeec67447ca"
|
|
||||||
integrity sha512-lobf6RWXb8c4uZ7Mdq0U12efYmpD1UFnyOWVJPTa3ukqZrMopav+2hdNu0hgBF0JIBFK9QgrBDfwYvh3DFJDAA==
|
|
||||||
|
|
||||||
balanced-match@^1.0.0:
|
balanced-match@^1.0.0:
|
||||||
version "1.0.2"
|
version "1.0.2"
|
||||||
resolved "https://registry.yarnpkg.com/balanced-match/-/balanced-match-1.0.2.tgz#e83e3a7e3f300b34cb9d87f615fa0cbf357690ee"
|
resolved "https://registry.yarnpkg.com/balanced-match/-/balanced-match-1.0.2.tgz#e83e3a7e3f300b34cb9d87f615fa0cbf357690ee"
|
||||||
@ -294,24 +289,11 @@ base64id@2.0.0, base64id@~2.0.0:
|
|||||||
resolved "https://registry.yarnpkg.com/base64id/-/base64id-2.0.0.tgz#2770ac6bc47d312af97a8bf9a634342e0cd25cb6"
|
resolved "https://registry.yarnpkg.com/base64id/-/base64id-2.0.0.tgz#2770ac6bc47d312af97a8bf9a634342e0cd25cb6"
|
||||||
integrity sha512-lGe34o6EHj9y3Kts9R4ZYs/Gr+6N7MCaMlIFA3F1R2O5/m7K06AxfSeO5530PEERE6/WyEg3lsuyw4GHlPZHog==
|
integrity sha512-lGe34o6EHj9y3Kts9R4ZYs/Gr+6N7MCaMlIFA3F1R2O5/m7K06AxfSeO5530PEERE6/WyEg3lsuyw4GHlPZHog==
|
||||||
|
|
||||||
benchmark@*:
|
|
||||||
version "2.1.4"
|
|
||||||
resolved "https://registry.yarnpkg.com/benchmark/-/benchmark-2.1.4.tgz#09f3de31c916425d498cc2ee565a0ebf3c2a5629"
|
|
||||||
integrity sha1-CfPeMckWQl1JjMLuVloOvzwqVik=
|
|
||||||
dependencies:
|
|
||||||
lodash "^4.17.4"
|
|
||||||
platform "^1.3.3"
|
|
||||||
|
|
||||||
binary-extensions@^2.0.0:
|
binary-extensions@^2.0.0:
|
||||||
version "2.2.0"
|
version "2.2.0"
|
||||||
resolved "https://registry.yarnpkg.com/binary-extensions/-/binary-extensions-2.2.0.tgz#75f502eeaf9ffde42fc98829645be4ea76bd9e2d"
|
resolved "https://registry.yarnpkg.com/binary-extensions/-/binary-extensions-2.2.0.tgz#75f502eeaf9ffde42fc98829645be4ea76bd9e2d"
|
||||||
integrity sha512-jDctJ/IVQbZoJykoeHbhXpOlNBqGNcwXJKJog42E5HDPUwQTSdjCHdihjj0DlnheQ7blbT6dHOafNAiS8ooQKA==
|
integrity sha512-jDctJ/IVQbZoJykoeHbhXpOlNBqGNcwXJKJog42E5HDPUwQTSdjCHdihjj0DlnheQ7blbT6dHOafNAiS8ooQKA==
|
||||||
|
|
||||||
binaryen@101.0.0:
|
|
||||||
version "101.0.0"
|
|
||||||
resolved "https://registry.yarnpkg.com/binaryen/-/binaryen-101.0.0.tgz#42a9e4cc7a22e2c1d75a31d28005a9b518b2c555"
|
|
||||||
integrity sha512-FRmVxvrR8jtcf0qcukNAPZDM3dZ2sc9GmA/hKxBI7k3fFzREKh1cAs+ruQi+ITTKz7u/AuFMuVnbJwTh0ef/HQ==
|
|
||||||
|
|
||||||
body-parser@^1.19.0:
|
body-parser@^1.19.0:
|
||||||
version "1.19.1"
|
version "1.19.1"
|
||||||
resolved "https://registry.yarnpkg.com/body-parser/-/body-parser-1.19.1.tgz#1499abbaa9274af3ecc9f6f10396c995943e31d4"
|
resolved "https://registry.yarnpkg.com/body-parser/-/body-parser-1.19.1.tgz#1499abbaa9274af3ecc9f6f10396c995943e31d4"
|
||||||
@ -599,14 +581,6 @@ dom-serialize@^2.2.1:
|
|||||||
extend "^3.0.0"
|
extend "^3.0.0"
|
||||||
void-elements "^2.0.0"
|
void-elements "^2.0.0"
|
||||||
|
|
||||||
dukat@0.5.8-rc.4:
|
|
||||||
version "0.5.8-rc.4"
|
|
||||||
resolved "https://registry.yarnpkg.com/dukat/-/dukat-0.5.8-rc.4.tgz#90384dcb50b14c26f0e99dae92b2dea44f5fce21"
|
|
||||||
integrity sha512-ZnMt6DGBjlVgK2uQamXfd7uP/AxH7RqI0BL9GLrrJb2gKdDxvJChWy+M9AQEaL+7/6TmxzJxFOsRiInY9oGWTA==
|
|
||||||
dependencies:
|
|
||||||
google-protobuf "3.12.2"
|
|
||||||
typescript "3.9.5"
|
|
||||||
|
|
||||||
ee-first@1.1.1:
|
ee-first@1.1.1:
|
||||||
version "1.1.1"
|
version "1.1.1"
|
||||||
resolved "https://registry.yarnpkg.com/ee-first/-/ee-first-1.1.1.tgz#590c61156b0ae2f4f0255732a158b266bc56b21d"
|
resolved "https://registry.yarnpkg.com/ee-first/-/ee-first-1.1.1.tgz#590c61156b0ae2f4f0255732a158b266bc56b21d"
|
||||||
@ -881,11 +855,6 @@ glob@^7.1.3, glob@^7.1.7:
|
|||||||
once "^1.3.0"
|
once "^1.3.0"
|
||||||
path-is-absolute "^1.0.0"
|
path-is-absolute "^1.0.0"
|
||||||
|
|
||||||
google-protobuf@3.12.2:
|
|
||||||
version "3.12.2"
|
|
||||||
resolved "https://registry.yarnpkg.com/google-protobuf/-/google-protobuf-3.12.2.tgz#50ce9f9b6281235724eb243d6a83e969a2176e53"
|
|
||||||
integrity sha512-4CZhpuRr1d6HjlyrxoXoocoGFnRYgKULgMtikMddA9ztRyYR59Aondv2FioyxWVamRo0rF2XpYawkTCBEQOSkA==
|
|
||||||
|
|
||||||
graceful-fs@^4.1.2, graceful-fs@^4.1.6, graceful-fs@^4.2.0, graceful-fs@^4.2.4, graceful-fs@^4.2.6:
|
graceful-fs@^4.1.2, graceful-fs@^4.1.6, graceful-fs@^4.2.0, graceful-fs@^4.2.4, graceful-fs@^4.2.6:
|
||||||
version "4.2.9"
|
version "4.2.9"
|
||||||
resolved "https://registry.yarnpkg.com/graceful-fs/-/graceful-fs-4.2.9.tgz#041b05df45755e587a24942279b9d113146e1c96"
|
resolved "https://registry.yarnpkg.com/graceful-fs/-/graceful-fs-4.2.9.tgz#041b05df45755e587a24942279b9d113146e1c96"
|
||||||
@ -1065,11 +1034,6 @@ jest-worker@^27.4.1:
|
|||||||
merge-stream "^2.0.0"
|
merge-stream "^2.0.0"
|
||||||
supports-color "^8.0.0"
|
supports-color "^8.0.0"
|
||||||
|
|
||||||
js-base64@3.6.1:
|
|
||||||
version "3.6.1"
|
|
||||||
resolved "https://registry.yarnpkg.com/js-base64/-/js-base64-3.6.1.tgz#555aae398b74694b4037af1f8a5a6209d170efbe"
|
|
||||||
integrity sha512-Frdq2+tRRGLQUIQOgsIGSCd1VePCS2fsddTG5dTCqR0JHgltXWfsxnY0gIXPoMeRmdom6Oyq+UMOFg5suduOjQ==
|
|
||||||
|
|
||||||
js-yaml@4.1.0:
|
js-yaml@4.1.0:
|
||||||
version "4.1.0"
|
version "4.1.0"
|
||||||
resolved "https://registry.yarnpkg.com/js-yaml/-/js-yaml-4.1.0.tgz#c1fb65f8f5017901cdd2c951864ba18458a10602"
|
resolved "https://registry.yarnpkg.com/js-yaml/-/js-yaml-4.1.0.tgz#c1fb65f8f5017901cdd2c951864ba18458a10602"
|
||||||
@ -1179,7 +1143,7 @@ locate-path@^6.0.0:
|
|||||||
dependencies:
|
dependencies:
|
||||||
p-locate "^5.0.0"
|
p-locate "^5.0.0"
|
||||||
|
|
||||||
lodash@^4.17.15, lodash@^4.17.21, lodash@^4.17.4:
|
lodash@^4.17.15, lodash@^4.17.21:
|
||||||
version "4.17.21"
|
version "4.17.21"
|
||||||
resolved "https://registry.yarnpkg.com/lodash/-/lodash-4.17.21.tgz#679591c564c3bffaae8454cf0b3df370c3d6911c"
|
resolved "https://registry.yarnpkg.com/lodash/-/lodash-4.17.21.tgz#679591c564c3bffaae8454cf0b3df370c3d6911c"
|
||||||
integrity sha512-v2kDEe57lecTulaDIuNTPy3Ry4gLGJ6Z1O3vE1krgXZNrsQ+LFTGHVxVjcXPs17LhbZVGedAJv8XZ1tvj5FvSg==
|
integrity sha512-v2kDEe57lecTulaDIuNTPy3Ry4gLGJ6Z1O3vE1krgXZNrsQ+LFTGHVxVjcXPs17LhbZVGedAJv8XZ1tvj5FvSg==
|
||||||
@ -1437,11 +1401,6 @@ pkg-dir@^4.2.0:
|
|||||||
dependencies:
|
dependencies:
|
||||||
find-up "^4.0.0"
|
find-up "^4.0.0"
|
||||||
|
|
||||||
platform@^1.3.3:
|
|
||||||
version "1.3.6"
|
|
||||||
resolved "https://registry.yarnpkg.com/platform/-/platform-1.3.6.tgz#48b4ce983164b209c2d45a107adb31f473a6e7a7"
|
|
||||||
integrity sha512-fnWVljUchTro6RiCFvCXBbNhJc2NijN7oIQxbwsyL0buWJPG85v81ehlHI9fXrJsMNgTofEoWIQeClKpgxFLrg==
|
|
||||||
|
|
||||||
postcss-modules-extract-imports@^3.0.0:
|
postcss-modules-extract-imports@^3.0.0:
|
||||||
version "3.0.0"
|
version "3.0.0"
|
||||||
resolved "https://registry.yarnpkg.com/postcss-modules-extract-imports/-/postcss-modules-extract-imports-3.0.0.tgz#cda1f047c0ae80c97dbe28c3e76a43b88025741d"
|
resolved "https://registry.yarnpkg.com/postcss-modules-extract-imports/-/postcss-modules-extract-imports-3.0.0.tgz#cda1f047c0ae80c97dbe28c3e76a43b88025741d"
|
||||||
@ -1696,14 +1655,6 @@ source-map-loader@3.0.0:
|
|||||||
iconv-lite "^0.6.2"
|
iconv-lite "^0.6.2"
|
||||||
source-map-js "^0.6.2"
|
source-map-js "^0.6.2"
|
||||||
|
|
||||||
source-map-support@0.5.20:
|
|
||||||
version "0.5.20"
|
|
||||||
resolved "https://registry.yarnpkg.com/source-map-support/-/source-map-support-0.5.20.tgz#12166089f8f5e5e8c56926b377633392dd2cb6c9"
|
|
||||||
integrity sha512-n1lZZ8Ve4ksRqizaBQgxXDgKwttHDhyfQjA6YZZn8+AroHbsIz+JjwxQDxbp+7y5OYCI8t1Yk7etjD9CRd2hIw==
|
|
||||||
dependencies:
|
|
||||||
buffer-from "^1.0.0"
|
|
||||||
source-map "^0.6.0"
|
|
||||||
|
|
||||||
source-map-support@~0.5.20:
|
source-map-support@~0.5.20:
|
||||||
version "0.5.21"
|
version "0.5.21"
|
||||||
resolved "https://registry.yarnpkg.com/source-map-support/-/source-map-support-0.5.21.tgz#04fe7c7f9e1ed2d662233c28cb2b35b9f63f6e4f"
|
resolved "https://registry.yarnpkg.com/source-map-support/-/source-map-support-0.5.21.tgz#04fe7c7f9e1ed2d662233c28cb2b35b9f63f6e4f"
|
||||||
@ -1838,11 +1789,6 @@ type-is@~1.6.18:
|
|||||||
media-typer "0.3.0"
|
media-typer "0.3.0"
|
||||||
mime-types "~2.1.24"
|
mime-types "~2.1.24"
|
||||||
|
|
||||||
typescript@3.9.5:
|
|
||||||
version "3.9.5"
|
|
||||||
resolved "https://registry.yarnpkg.com/typescript/-/typescript-3.9.5.tgz#586f0dba300cde8be52dd1ac4f7e1009c1b13f36"
|
|
||||||
integrity sha512-hSAifV3k+i6lEoCJ2k6R2Z/rp/H3+8sdmcn5NrS3/3kE7+RyZXm9aqvxWqjEXHAd8b0pShatpcdMTvEdvAJltQ==
|
|
||||||
|
|
||||||
ua-parser-js@^0.7.28:
|
ua-parser-js@^0.7.28:
|
||||||
version "0.7.31"
|
version "0.7.31"
|
||||||
resolved "https://registry.yarnpkg.com/ua-parser-js/-/ua-parser-js-0.7.31.tgz#649a656b191dffab4f21d5e053e27ca17cbff5c6"
|
resolved "https://registry.yarnpkg.com/ua-parser-js/-/ua-parser-js-0.7.31.tgz#649a656b191dffab4f21d5e053e27ca17cbff5c6"
|
||||||
|
Loading…
Reference in New Issue
Block a user