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

Feature: Polynomials and rational functions #469

Merged
lounres merged 132 commits from feature/polynomials into dev 2022-07-28 18:04:06 +03:00
Conversation 19 Commits 132 Files Changed 51 +2007 -9
11 changed files with 2007 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit b50d8dcd23 - Show all commits

385
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/ListPolynomial.kt Normal file
View File

@ -0,0 +1,385 @@
/*
* 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.Ring
import space.kscience.kmath.operations.ScaleOperations
import space.kscience.kmath.operations.invoke
import kotlin.math.max
import kotlin.math.min
/**
* Represents univariate polynomial that stores its coefficients in a [List].
*
* @param coefficients constant is the leftmost coefficient.
*/
public data class ListPolynomial<C>(
/**
* List that contains coefficients of the polynomial. Every monomial `a x^d` is stored as a coefficient `a` placed
* into the list at index `d`. For example, coefficients of a polynomial `5 x^2 - 6` can be represented as
* ```
* listOf(
* -6, // -6 +
* 0, // 0 x +
* 5, // 5 x^2
* )
* ```
* and also as
* ```
* listOf(
* -6, // -6 +
* 0, // 0 x +
* 5, // 5 x^2
* 0, // 0 x^3
* 0, // 0 x^4
* )
* ```
* It is not prohibited to put extra zeros at end of the list (as for `0x^3` and `0x^4` in the example). But the
* longer the coefficients list the worse performance of arithmetical operations performed on it. Thus, it is
* recommended not to put (or even to remove) extra (or useless) coefficients at the end of the coefficients list.
*/
public val coefficients: List<C>
) : Polynomial<C> {
override fun toString(): String = "ListPolynomial$coefficients"
}
/**
* Arithmetic context for univariate polynomials with coefficients stored as a [List] constructed with the given [ring]
* of constants.
*
* @param C the type of constants. Polynomials have them a coefficients in their terms.
* @param A type of provided underlying ring of constants. It's [Ring] of [C].
* @param ring underlying ring of constants of type [A].
*/
public open class ListPolynomialSpace<C, A : Ring<C>>(
public override val ring: A,
) : PolynomialSpaceOverRing<C, ListPolynomial<C>, A> {
/**
* Returns sum of the polynomial and the integer represented as a polynomial.
*
* The operation is equivalent to adding [other] copies of unit polynomial to [this].
*/
public override operator fun ListPolynomial<C>.plus(other: Int): ListPolynomial<C> =
if (other == 0) this
else
ListPolynomial(
coefficients
.toMutableList()
.apply {
val result = getOrElse(0) { constantZero } + other
if(size == 0) add(result)
else this[0] = result
}
)
/**
* Returns difference between the polynomial and the integer represented as a polynomial.
*
* The operation is equivalent to subtraction [other] copies of unit polynomial from [this].
*/
public override operator fun ListPolynomial<C>.minus(other: Int): ListPolynomial<C> =
if (other == 0) this
else
ListPolynomial(
coefficients
.toMutableList()
.apply {
val result = getOrElse(0) { constantZero } - other
if(size == 0) add(result)
else this[0] = result
}
)
/**
* Returns product of the polynomial and the integer represented as a polynomial.
*
* The operation is equivalent to sum of [other] copies of [this].
*/
public override operator fun ListPolynomial<C>.times(other: Int): ListPolynomial<C> =
if (other == 0) zero
else ListPolynomial(
coefficients
.toMutableList()
.apply {
for (deg in indices) this[deg] = this[deg] * other
}
)
/**
* Returns sum of the integer represented as a polynomial and the polynomial.
*
* The operation is equivalent to adding [this] copies of unit polynomial to [other].
*/
public override operator fun Int.plus(other: ListPolynomial<C>): ListPolynomial<C> =
if (this == 0) other
else
ListPolynomial(
other.coefficients
.toMutableList()
.apply {
val result = this@plus + getOrElse(0) { constantZero }
if(size == 0) add(result)
else this[0] = result
}
)
/**
* Returns difference between the integer represented as a polynomial and the polynomial.
*
* The operation is equivalent to subtraction [this] copies of unit polynomial from [other].
*/
public override operator fun Int.minus(other: ListPolynomial<C>): ListPolynomial<C> =
if (this == 0) other
else
ListPolynomial(
other.coefficients
.toMutableList()
.apply {
forEachIndexed { index, c -> if (index != 0) this[index] = -c }
val result = this@minus - getOrElse(0) { constantZero }
if(size == 0) add(result)
else this[0] = result
}
)
/**
* Returns product of the integer represented as a polynomial and the polynomial.
*
* The operation is equivalent to sum of [this] copies of [other].
*/
public override operator fun Int.times(other: ListPolynomial<C>): ListPolynomial<C> =
if (this == 0) zero
else ListPolynomial(
other.coefficients
.toMutableList()
.apply {
for (deg in indices) this[deg] = this@times * this[deg]
}
)
/**
* Converts the integer [value] to polynomial.
*/
public override fun number(value: Int): ListPolynomial<C> = number(constantNumber(value))
/**
* Returns sum of the constant represented as a polynomial and the polynomial.
*/
public override operator fun C.plus(other: ListPolynomial<C>): ListPolynomial<C> =
with(other.coefficients) {
if (isEmpty()) ListPolynomial(listOf(this@plus))
else ListPolynomial(
toMutableList()
.apply {
val result = if (size == 0) this@plus else this@plus + get(0)
if(size == 0) add(result)
else this[0] = result
}
)
}
/**
* Returns difference between the constant represented as a polynomial and the polynomial.
*/
public override operator fun C.minus(other: ListPolynomial<C>): ListPolynomial<C> =
with(other.coefficients) {
if (isEmpty()) ListPolynomial(listOf(this@minus))
else ListPolynomial(
toMutableList()
.apply {
forEachIndexed { index, c -> if (index != 0) this[index] = -c }
val result = if (size == 0) this@minus else this@minus - get(0)
if(size == 0) add(result)
else this[0] = result
}
)
}
/**
* Returns product of the constant represented as a polynomial and the polynomial.
*/
public override operator fun C.times(other: ListPolynomial<C>): ListPolynomial<C> =
ListPolynomial(
other.coefficients
.toMutableList()
.apply {
for (deg in indices) this[deg] = this@times * this[deg]
}
)
/**
* Returns sum of the constant represented as a polynomial and the polynomial.
*/
public override operator fun ListPolynomial<C>.plus(other: C): ListPolynomial<C> =
with(coefficients) {
if (isEmpty()) ListPolynomial(listOf(other))
else ListPolynomial(
toMutableList()
.apply {
val result = if (size == 0) other else get(0) + other
if(size == 0) add(result)
else this[0] = result
}
)
}
/**
* Returns difference between the constant represented as a polynomial and the polynomial.
*/
public override operator fun ListPolynomial<C>.minus(other: C): ListPolynomial<C> =
with(coefficients) {
if (isEmpty()) ListPolynomial(listOf(-other))
else ListPolynomial(
toMutableList()
.apply {
val result = if (size == 0) other else get(0) - other
if(size == 0) add(result)
else this[0] = result
}
)
}
/**
* Returns product of the constant represented as a polynomial and the polynomial.
*/
public override operator fun ListPolynomial<C>.times(other: C): ListPolynomial<C> =
ListPolynomial(
coefficients
.toMutableList()
.apply {
for (deg in indices) this[deg] = this[deg] * other
}
)
/**
* Converts the constant [value] to polynomial.
*/
public override fun number(value: C): ListPolynomial<C> = ListPolynomial(listOf(value))
/**
* Returns negation of the polynomial.
*/
public override operator fun ListPolynomial<C>.unaryMinus(): ListPolynomial<C> =
ListPolynomial(coefficients.map { -it })
/**
* Returns sum of the polynomials.
*/
public override operator fun ListPolynomial<C>.plus(other: ListPolynomial<C>): ListPolynomial<C> {
val thisDegree = degree
val otherDegree = other.degree
return ListPolynomial(
List(max(thisDegree, otherDegree) + 1) {
when {
it > thisDegree -> other.coefficients[it]
it > otherDegree -> coefficients[it]
else -> coefficients[it] + other.coefficients[it]
}
}
)
}
/**
* Returns difference of the polynomials.
*/
public override operator fun ListPolynomial<C>.minus(other: ListPolynomial<C>): ListPolynomial<C> {
val thisDegree = degree
val otherDegree = other.degree
return ListPolynomial(
List(max(thisDegree, otherDegree) + 1) {
when {
it > thisDegree -> -other.coefficients[it]
it > otherDegree -> coefficients[it]
else -> coefficients[it] - other.coefficients[it]
}
}
)
}
/**
* Returns product of the polynomials.
*/
public override operator fun ListPolynomial<C>.times(other: ListPolynomial<C>): ListPolynomial<C> {
val thisDegree = degree
val otherDegree = other.degree
return ListPolynomial(
List(thisDegree + otherDegree + 1) { d ->
(max(0, d - otherDegree)..min(thisDegree, d))
.map { coefficients[it] * other.coefficients[d - it] }
.reduce { acc, rational -> acc + rational }
}
)
}
/**
* Instance of zero polynomial (zero of the polynomial ring).
*/
override val zero: ListPolynomial<C> = ListPolynomial(emptyList())
/**
* Instance of unit polynomial (unit of the polynomial ring).
*/
override val one: ListPolynomial<C> by lazy { ListPolynomial(listOf(constantOne)) }
/**
* Degree of the polynomial, [see also](https://en.wikipedia.org/wiki/Degree_of_a_polynomial). If the polynomial is
* zero, degree is -1.
*/
public override val ListPolynomial<C>.degree: Int get() = coefficients.lastIndex
// TODO: When context receivers will be ready move all of this substitutions and invocations to utilities with
// [ListPolynomialSpace] as a context receiver
/**
* Evaluates value of [this] polynomial on provided argument.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.substitute(argument: C): C = this.substitute(ring, argument)
/**
* Substitutes provided polynomial [argument] into [this] polynomial.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.substitute(argument: ListPolynomial<C>): ListPolynomial<C> = this.substitute(ring, argument)
/**
* Represent [this] polynomial as a regular context-less function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.asFunction(): (C) -> C = { this.substitute(ring, it) }
/**
* Represent [this] polynomial as a regular context-less function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.asFunctionOfConstant(): (C) -> C = { this.substitute(ring, it) }
/**
* Represent [this] polynomial as a regular context-less function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.asFunctionOfPolynomial(): (ListPolynomial<C>) -> ListPolynomial<C> = { this.substitute(ring, it) }
/**
* Evaluates value of [this] polynomial on provided [argument].
*/
@Suppress("NOTHING_TO_INLINE")
public inline operator fun ListPolynomial<C>.invoke(argument: C): C = this.substitute(ring, argument)
/**
* Evaluates value of [this] polynomial on provided [argument].
*/
@Suppress("NOTHING_TO_INLINE")
public inline operator fun ListPolynomial<C>.invoke(argument: ListPolynomial<C>): ListPolynomial<C> = this.substitute(ring, argument)
}
/**
* Space of polynomials constructed over ring.
*
* @param C the type of constants. Polynomials have them as a coefficients in their terms.
* @param A type of underlying ring of constants. It's [Ring] of [C].
* @param ring underlying ring of constants of type [A].
*/
public class ScalableListPolynomialSpace<C, A>(
ring: A,
) : ListPolynomialSpace<C, A>(ring), ScaleOperations<ListPolynomial<C>> where A : Ring<C>, A : ScaleOperations<C> {
override fun scale(a: ListPolynomial<C>, value: Double): ListPolynomial<C> =
ring { ListPolynomial(a.coefficients.map { scale(it, value) }) }
}

137
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/ListRationalFunction.kt Normal file
View File

@ -0,0 +1,137 @@
/*
* 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.Ring
/**
* Represents rational function that stores its numerator and denominator as [ListPolynomial]s.
*/
public data class ListRationalFunction<C>(
public override val numerator: ListPolynomial<C>,
public override val denominator: ListPolynomial<C>
) : RationalFunction<C, ListPolynomial<C>> {
override fun toString(): String = "ListRationalFunction${numerator.coefficients}/${denominator.coefficients}"
}
/**
* Arithmetic context for univariate rational functions with numerator and denominator represented as [ListPolynomial]s.
*
* @param C the type of constants. Polynomials have them a coefficients in their terms.
* @param A type of provided underlying ring of constants. It's [Ring] of [C].
* @param ring underlying ring of constants of type [A].
*/
public class ListRationalFunctionSpace<C, A : Ring<C>> (
public val ring: A,
) :
RationalFunctionalSpaceOverPolynomialSpace<
C,
ListPolynomial<C>,
ListRationalFunction<C>,
ListPolynomialSpace<C, A>,
>,
PolynomialSpaceOfFractions<
C,
ListPolynomial<C>,
ListRationalFunction<C>,
>() {
/**
* Underlying polynomial ring. Its polynomial operations are inherited by local polynomial operations.
*/
override val polynomialRing : ListPolynomialSpace<C, A> = ListPolynomialSpace(ring)
/**
* Constructor of [ListRationalFunction] from numerator and denominator [ListPolynomial].
*/
override fun constructRationalFunction(numerator: ListPolynomial<C>, denominator: ListPolynomial<C>): ListRationalFunction<C> =
ListRationalFunction(numerator, denominator)
// TODO: When context receivers will be ready move all of this substitutions and invocations to utilities with
// [ListPolynomialSpace] as a context receiver
/**
* Evaluates value of [this] polynomial on provided argument.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.substitute(argument: C): C = this.substitute(ring, argument)
/**
* Substitutes provided polynomial [argument] into [this] polynomial.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.substitute(argument: ListPolynomial<C>): ListPolynomial<C> = this.substitute(ring, argument)
/**
* Substitutes provided rational function [argument] into [this] polynomial.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.substitute(argument: ListRationalFunction<C>): ListRationalFunction<C> = this.substitute(ring, argument)
/**
* Substitutes provided polynomial [argument] into [this] rational function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListRationalFunction<C>.substitute(argument: ListPolynomial<C>): ListRationalFunction<C> = this.substitute(ring, argument)
/**
* Substitutes provided rational function [argument] into [this] rational function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListRationalFunction<C>.substitute(argument: ListRationalFunction<C>): ListRationalFunction<C> = this.substitute(ring, argument)
/**
* Represent [this] polynomial as a regular context-less function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.asFunction(): (C) -> C = { this.substitute(ring, it) }
/**
* Represent [this] polynomial as a regular context-less function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.asFunctionOfConstant(): (C) -> C = { this.substitute(ring, it) }
/**
* Represent [this] polynomial as a regular context-less function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.asFunctionOfPolynomial(): (ListPolynomial<C>) -> ListPolynomial<C> = { this.substitute(ring, it) }
/**
* Represent [this] polynomial as a regular context-less function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListPolynomial<C>.asFunctionOfRationalFunction(): (ListRationalFunction<C>) -> ListRationalFunction<C> = { this.substitute(ring, it) }
/**
* Represent [this] rational function as a regular context-less function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListRationalFunction<C>.asFunctionOfPolynomial(): (ListPolynomial<C>) -> ListRationalFunction<C> = { this.substitute(ring, it) }
/**
* Represent [this] rational function as a regular context-less function.
*/
@Suppress("NOTHING_TO_INLINE")
public inline fun ListRationalFunction<C>.asFunctionOfRationalFunction(): (ListRationalFunction<C>) -> ListRationalFunction<C> = { this.substitute(ring, it) }
/**
* Evaluates value of [this] polynomial on provided argument.
*/
@Suppress("NOTHING_TO_INLINE")
public inline operator fun ListPolynomial<C>.invoke(argument: C): C = this.substitute(ring, argument)
/**
* Evaluates value of [this] polynomial on provided argument.
*/
@Suppress("NOTHING_TO_INLINE")
public inline operator fun ListPolynomial<C>.invoke(argument: ListPolynomial<C>): ListPolynomial<C> = this.substitute(ring, argument)
/**
* Evaluates value of [this] polynomial on provided argument.
*/
@Suppress("NOTHING_TO_INLINE")
public inline operator fun ListPolynomial<C>.invoke(argument: ListRationalFunction<C>): ListRationalFunction<C> = this.substitute(ring, argument)
/**
* Evaluates value of [this] rational function on provided argument.
*/
@Suppress("NOTHING_TO_INLINE")
public inline operator fun ListRationalFunction<C>.invoke(argument: ListPolynomial<C>): ListRationalFunction<C> = this.substitute(ring, argument)
/**
* Evaluates value of [this] rational function on provided argument.
*/
@Suppress("NOTHING_TO_INLINE")
public inline operator fun ListRationalFunction<C>.invoke(argument: ListRationalFunction<C>): ListRationalFunction<C> = this.substitute(ring, argument)
}

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

@ -435,10 +435,12 @@ public interface MultivariatePolynomialSpace<C, V, P: Polynomial<C>>: Polynomial
/**
* Represents the [variable] as a monic monomial.
*/
@JvmName("numberVariable")
public fun number(variable: V): P = +variable
/**
* Represents the variable as a monic monomial.
*/
@JvmName("asPolynomialVariable")
public fun V.asPolynomial(): P = number(this)
/**

44
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/RationalFunction.kt
View File

@ -834,7 +834,12 @@ public abstract class PolynomialSpaceOfFractions<
numerator * other,
denominator
)
/**
* Returns quotient of the rational function and the integer represented as a rational function.
*
* The operation is equivalent to creating a new rational function by preserving numerator of [this] and
* multiplication denominator of [this] to [other].
*/
public override operator fun R.div(other: Int): R =
constructRationalFunction(
numerator,
@ -871,7 +876,12 @@ public abstract class PolynomialSpaceOfFractions<
this * other.numerator,
other.denominator
)
/**
* Returns quotient of the integer represented as a rational function and the rational function.
*
* The operation is equivalent to creating a new rational function which numerator is [this] times denominator of
* [other] and which denominator is [other]'s numerator.
*/
public override operator fun Int.div(other: R): R =
constructRationalFunction(
this * other.denominator,
@ -912,7 +922,9 @@ public abstract class PolynomialSpaceOfFractions<
this * other.numerator,
other.denominator
)
/**
* Returns quotient of the constant represented as a polynomial and the rational function.
*/
public override operator fun C.div(other: R): R =
constructRationalFunction(
this * other.denominator,
@ -943,7 +955,9 @@ public abstract class PolynomialSpaceOfFractions<
numerator * other,
denominator
)
/**
* Returns quotient of the rational function and the constant represented as a rational function.
*/
public override operator fun R.div(other: C): R =
constructRationalFunction(
numerator,
@ -979,7 +993,9 @@ public abstract class PolynomialSpaceOfFractions<
this * other.numerator,
other.denominator
)
/**
* Returns quotient of the polynomial represented as a polynomial and the rational function.
*/
public override operator fun P.div(other: R): R =
constructRationalFunction(
this * other.denominator,
@ -1010,7 +1026,9 @@ public abstract class PolynomialSpaceOfFractions<
numerator * other,
denominator
)
/**
* Returns quotient of the rational function and the polynomial represented as a rational function.
*/
public override operator fun R.div(other: P): R =
constructRationalFunction(
numerator,
@ -1050,7 +1068,9 @@ public abstract class PolynomialSpaceOfFractions<
numerator * other.numerator,
denominator * other.denominator
)
/**
* Returns quotient of the rational functions.
*/
public override operator fun R.div(other: R): R =
constructRationalFunction(
numerator * other.denominator,
@ -1060,12 +1080,12 @@ public abstract class PolynomialSpaceOfFractions<
/**
* Instance of zero rational function (zero of the rational functions ring).
*/
public override val zero: R get() = constructRationalFunction(polynomialZero)
public override val zero: R by lazy { constructRationalFunction(polynomialZero) }
/**
* Instance of unit polynomial (unit of the rational functions ring).
*/
public override val one: R get() = constructRationalFunction(polynomialOne)
public override val one: R by lazy { constructRationalFunction(polynomialOne) }
}
/**
@ -1177,19 +1197,23 @@ public interface MultivariateRationalFunctionalSpace<
/**
* Represents the [variable] as a monic monomial.
*/
@JvmName("polynomialNumberVariable")
public fun polynomialNumber(variable: V): P = +variable
/**
* Represents the variable as a monic monomial.
*/
@JvmName("asPolynomialVariable")
public fun V.asPolynomial(): P = polynomialNumber(this)
/**
* Represents the [variable] as a rational function.
*/
@JvmName("numberVariable")
public fun number(variable: V): R = number(polynomialNumber(variable))
/**
* Represents the variable as a rational function.
*/
@JvmName("asRationalFunctionVariable")
public fun V.asRationalFunction(): R = number(this)
/**
@ -1403,10 +1427,12 @@ public interface MultivariateRationalFunctionalSpaceOverMultivariatePolynomialSp
/**
* Represents the [variable] as a monic monomial.
*/
@JvmName("polynomialNumberVariable")
public override fun polynomialNumber(variable: V): P = polynomialRing { number(variable) }
/**
* Represents the variable as a monic monomial.
*/
@JvmName("asPolynomialVariable")
public override fun V.asPolynomial(): P = polynomialRing { this@asPolynomial.asPolynomial() }
/**

268
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/listUtil.kt Normal file
View File

@ -0,0 +1,268 @@
/*
* 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.math.max
import kotlin.math.pow
/**
* Creates a [ListPolynomialSpace] over a received ring.
*/
public fun <C, A : Ring<C>> A.listPolynomialSpace(): ListPolynomialSpace<C, A> =
ListPolynomialSpace(this)
/**
* Creates a [ListPolynomialSpace]'s scope over a received ring.
*/ // TODO: When context will be ready move [ListPolynomialSpace] and add [A] to context receivers of [block]
public inline fun <C, A : Ring<C>, R> A.listPolynomialSpace(block: ListPolynomialSpace<C, A>.() -> R): R {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return ListPolynomialSpace(this).block()
}
/**
* Creates a [ScalableListPolynomialSpace] over a received scalable ring.
*/
public fun <C, A> A.scalableListPolynomialSpace(): ScalableListPolynomialSpace<C, A> where A : Ring<C>, A : ScaleOperations<C> =
ScalableListPolynomialSpace(this)
/**
* Creates a [ScalableListPolynomialSpace]'s scope over a received scalable ring.
*/ // TODO: When context will be ready move [ListPolynomialSpace] and add [A] to context receivers of [block]
public inline fun <C, A, R> A.scalableListPolynomialSpace(block: ScalableListPolynomialSpace<C, A>.() -> R): R where A : Ring<C>, A : ScaleOperations<C> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return ScalableListPolynomialSpace(this).block()
}
/**
* Creates a [ListRationalFunctionSpace] over a received ring.
*/
public fun <C, A : Ring<C>> A.listRationalFunctionSpace(): ListRationalFunctionSpace<C, A> =
ListRationalFunctionSpace(this)
/**
* Creates a [ListRationalFunctionSpace]'s scope over a received ring.
*/ // TODO: When context will be ready move [ListRationalFunctionSpace] and add [A] to context receivers of [block]
public inline fun <C, A : Ring<C>, R> A.listRationalFunctionSpace(block: ListRationalFunctionSpace<C, A>.() -> R): R {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return ListRationalFunctionSpace(this).block()
}
/**
* Evaluates value of [this] Double polynomial on provided Double argument.
*/
public fun ListPolynomial<Double>.substitute(arg: Double): Double =
coefficients.reduceIndexedOrNull { index, acc, c ->
acc + c * arg.pow(index)
} ?: .0
/**
* Evaluates value of [this] polynomial on provided argument.
*
* It is an implementation of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
*/
public fun <C> ListPolynomial<C>.substitute(ring: Ring<C>, arg: C): C = ring {
if (coefficients.isEmpty()) return zero
var result: C = coefficients.last()
for (j in coefficients.size - 2 downTo 0) {
result = (arg * result) + coefficients[j]
}
return result
}
/**
* Substitutes provided polynomial [arg] into [this] polynomial.
*
* It is an implementation of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
*/ // TODO: To optimize boxing
public fun <C> ListPolynomial<C>.substitute(ring: Ring<C>, arg: ListPolynomial<C>) : ListPolynomial<C> =
ring.listPolynomialSpace {
if (coefficients.isEmpty()) return zero
var result: ListPolynomial<C> = coefficients.last().asPolynomial()
for (j in coefficients.size - 2 downTo 0) {
result = (arg * result) + coefficients[j]
}
return result
}
/**
* Substitutes provided rational function [arg] into [this] polynomial.
*
* It is an implementation of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
*/ // TODO: To optimize boxing
public fun <C> ListPolynomial<C>.substitute(ring: Ring<C>, arg: ListRationalFunction<C>) : ListRationalFunction<C> =
ring.listRationalFunctionSpace {
if (coefficients.isEmpty()) return zero
var result: ListRationalFunction<C> = coefficients.last().asRationalFunction()
for (j in coefficients.size - 2 downTo 0) {
result = (arg * result) + coefficients[j]
}
return result
}
/**
* Evaluates value of [this] Double rational function in provided Double argument.
*/
public fun ListRationalFunction<Double>.substitute(arg: Double): Double =
numerator.substitute(arg) / denominator.substitute(arg)
/**
* Evaluates value of [this] polynomial for provided argument.
*
* It is an implementation of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
*/
public fun <C> ListRationalFunction<C>.substitute(ring: Field<C>, arg: C): C = ring {
numerator.substitute(ring, arg) / denominator.substitute(ring, arg)
}
/**
* Substitutes provided polynomial [arg] into [this] rational function.
*/ // TODO: To optimize boxing
public fun <C> ListRationalFunction<C>.substitute(ring: Ring<C>, arg: ListPolynomial<C>) : ListRationalFunction<C> =
ring.listRationalFunctionSpace {
numerator.substitute(ring, arg) / denominator.substitute(ring, arg)
}
/**
* Substitutes provided rational function [arg] into [this] rational function.
*/ // TODO: To optimize boxing
public fun <C> ListRationalFunction<C>.substitute(ring: Ring<C>, arg: ListRationalFunction<C>) : ListRationalFunction<C> =
ring.listRationalFunctionSpace {
numerator.substitute(ring, arg) / denominator.substitute(ring, arg)
}
/**
* Represent [this] polynomial as a regular context-less function.
*/
public fun <C, A : Ring<C>> ListPolynomial<C>.asFunctionOver(ring: A): (C) -> C = { substitute(ring, it) }
/**
* Represent [this] polynomial as a regular context-less function.
*/
public fun <C, A : Ring<C>> ListPolynomial<C>.asPolynomialFunctionOver(ring: A): (ListPolynomial<C>) -> ListPolynomial<C> = { substitute(ring, it) }
/**
* Represent [this] polynomial as a regular context-less function.
*/
public fun <C, A : Ring<C>> ListPolynomial<C>.asFunctionOfRationalFunctionOver(ring: A): (ListPolynomial<C>) -> ListPolynomial<C> = { substitute(ring, it) }
/**
* Represent [this] rational function as a regular context-less function.
*/
public fun <C, A : Field<C>> ListRationalFunction<C>.asFunctionOver(ring: A): (C) -> C = { substitute(ring, it) }
/**
* Represent [this] rational function as a regular context-less function.
*/
public fun <C, A : Ring<C>> ListRationalFunction<C>.asPolynomialFunctionOver(ring: A): (ListPolynomial<C>) -> ListRationalFunction<C> = { substitute(ring, it) }
/**
* Represent [this] rational function as a regular context-less function.
*/
public fun <C, A : Ring<C>> ListRationalFunction<C>.asFunctionOfRationalFunctionOver(ring: A): (ListPolynomial<C>) -> ListRationalFunction<C> = { substitute(ring, it) }
/**
* Returns algebraic derivative of received polynomial.
*/
@UnstableKMathAPI
public fun <C, A> ListPolynomial<C>.derivative(
ring: A,
): ListPolynomial<C> where A : Ring<C>, A : NumericAlgebra<C> = ring {
ListPolynomial(
buildList(max(0, coefficients.size - 1)) {
for (deg in 1 .. coefficients.lastIndex) add(number(deg) * coefficients[deg])
}
)
}
/**
* Returns algebraic derivative of received polynomial of specified [order]. The [order] should be non-negative integer.
*/
@UnstableKMathAPI
public fun <C, A> ListPolynomial<C>.nthDerivative(
ring: A,
order: Int,
): ListPolynomial<C> where A : Ring<C>, A : NumericAlgebra<C> = ring {
require(order >= 0) { "Order of derivative must be non-negative" }
ListPolynomial(
buildList(max(0, coefficients.size - order)) {
for (deg in order.. coefficients.lastIndex)
add((deg - order + 1 .. deg).fold(coefficients[deg]) { acc, d -> acc * number(d) })
}
)
}
/**
* Returns algebraic antiderivative of received polynomial.
*/
@UnstableKMathAPI
public fun <C, A> ListPolynomial<C>.antiderivative(
ring: A,
): ListPolynomial<C> where A : Field<C>, A : NumericAlgebra<C> = ring {
ListPolynomial(
buildList(coefficients.size + 1) {
add(zero)
coefficients.mapIndexedTo(this) { index, t -> t / number(index + 1) }
}
)
}
/**
* Returns algebraic antiderivative of received polynomial of specified [order]. The [order] should be non-negative integer.
*/
@UnstableKMathAPI
public fun <C, A> ListPolynomial<C>.nthAntiderivative(
ring: A,
order: Int,
): ListPolynomial<C> where A : Field<C>, A : NumericAlgebra<C> = ring {
require(order >= 0) { "Order of antiderivative must be non-negative" }
ListPolynomial(
buildList(coefficients.size + order) {
repeat(order) { add(zero) }
coefficients.mapIndexedTo(this) { index, c -> (1..order).fold(c) { acc, i -> acc / number(index + i) } }
}
)
}
/**
* Computes a definite integral of [this] polynomial in the specified [range].
*/
@UnstableKMathAPI
public fun <C : Comparable<C>> ListPolynomial<C>.integrate(
ring: Field<C>,
range: ClosedRange<C>,
): C = ring {
val antiderivative = antiderivative(ring)
antiderivative.substitute(ring, range.endInclusive) - antiderivative.substitute(ring, range.start)
}
/**
* Returns algebraic derivative of received rational function.
*/
@UnstableKMathAPI
public fun <C, A> ListRationalFunction<C>.derivative(
ring: A,
): ListRationalFunction<C> where A : Ring<C>, A : NumericAlgebra<C> = ring.listRationalFunctionSpace {
ListRationalFunction(
numerator.derivative(ring) * denominator - numerator * denominator.derivative(ring),
denominator * denominator
)
}
/**
* Returns algebraic derivative of received rational function of specified [order]. The [order] should be non-negative integer.
*/
@UnstableKMathAPI
public tailrec fun <C, A> ListRationalFunction<C>.nthDerivative(
ring: A,
order: Int,
): ListRationalFunction<C> where A : Ring<C>, A : NumericAlgebra<C> =
if (order == 0) this else derivative(ring).nthDerivative(ring, order - 1)

250
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/listUtilOptimized.kt Normal file
View File

@ -0,0 +1,250 @@
/*
* 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.Ring
import space.kscience.kmath.operations.invoke
import kotlin.math.max
import kotlin.math.min
// TODO: Optimized copies of substitution and invocation
@UnstablePolynomialBoxingOptimization
@Suppress("NOTHING_TO_INLINE")
internal inline fun <C> copyTo(
origin: List<C>,
originDegree: Int,
target: MutableList<C>,
) {
for (deg in 0 .. originDegree) target[deg] = origin[deg]
}
@UnstablePolynomialBoxingOptimization
@Suppress("NOTHING_TO_INLINE")
internal inline fun <C> multiplyAddingToUpdater(
ring: Ring<C>,
multiplicand: MutableList<C>,
multiplicandDegree: Int,
multiplier: List<C>,
multiplierDegree: Int,
updater: MutableList<C>,
zero: C,
) {
multiplyAddingTo(
ring = ring,
multiplicand = multiplicand,
multiplicandDegree = multiplicandDegree,
multiplier = multiplier,
multiplierDegree = multiplierDegree,
target = updater
)
for (updateDeg in 0 .. multiplicandDegree + multiplierDegree) {
multiplicand[updateDeg] = updater[updateDeg]
updater[updateDeg] = zero
}
}
@UnstablePolynomialBoxingOptimization
@Suppress("NOTHING_TO_INLINE")
internal inline fun <C> multiplyAddingTo(
ring: Ring<C>,
multiplicand: List<C>,
multiplicandDegree: Int,
multiplier: List<C>,
multiplierDegree: Int,
target: MutableList<C>
) = ring {
for (d in 0 .. multiplicandDegree + multiplierDegree)
for (k in max(0, d - multiplierDegree)..min(multiplicandDegree, d))
target[d] += multiplicand[k] * multiplier[d - k]
}
@UnstablePolynomialBoxingOptimization
public fun <C> ListPolynomial<C>.substitute2(ring: Ring<C>, arg: ListPolynomial<C>) : ListPolynomial<C> = ring {
if (coefficients.isEmpty()) return ListPolynomial(emptyList())
val thisDegree = coefficients.lastIndex
if (thisDegree == -1) return ListPolynomial(emptyList())
val argDegree = arg.coefficients.lastIndex
if (argDegree == -1) return coefficients[0].asListPolynomial()
val constantZero = zero
val resultCoefs: MutableList<C> = MutableList(thisDegree * argDegree + 1) { constantZero }
resultCoefs[0] = coefficients[thisDegree]
val resultCoefsUpdate: MutableList<C> = MutableList(thisDegree * argDegree + 1) { constantZero }
var resultDegree = 0
for (deg in thisDegree - 1 downTo 0) {
resultCoefsUpdate[0] = coefficients[deg]
multiplyAddingToUpdater(
ring = ring,
multiplicand = resultCoefs,
multiplicandDegree = resultDegree,
multiplier = arg.coefficients,
multiplierDegree = argDegree,
updater = resultCoefsUpdate,
zero = constantZero
)
resultDegree += argDegree
}
return ListPolynomial<C>(resultCoefs)
}
/**
* Returns numerator (polynomial) of rational function gotten by substitution rational function [arg] to the polynomial instance.
* More concrete, if [arg] is a fraction `f(x)/g(x)` and the receiving instance is `p(x)`, then
* ```
* p(f/g) * g^deg(p)
* ```
* is returned.
*
* Used in [ListPolynomial.substitute] and [ListRationalFunction.substitute] for performance optimisation.
*/ // TODO: Дописать
@UnstablePolynomialBoxingOptimization
internal fun <C> ListPolynomial<C>.substituteRationalFunctionTakeNumerator(ring: Ring<C>, arg: ListRationalFunction<C>): ListPolynomial<C> = ring {
if (coefficients.isEmpty()) return ListPolynomial(emptyList())
val thisDegree = coefficients.lastIndex
if (thisDegree == -1) return ListPolynomial(emptyList())
val thisDegreeLog2 = 31 - thisDegree.countLeadingZeroBits()
val numeratorDegree = arg.numerator.coefficients.lastIndex
val denominatorDegree = arg.denominator.coefficients.lastIndex
val argDegree = max(numeratorDegree, denominatorDegree)
val constantZero = zero
val powersOf2 = buildList<Int>(thisDegreeLog2 + 1) {
var result = 1
for (exp in 0 .. thisDegreeLog2) {
add(result)
result = result shl 1
}
}
val hashes = powersOf2.runningReduce { acc, i -> acc + i }
val numeratorPowers = buildList<List<C>>(thisDegreeLog2 + 1) {
add(arg.numerator.coefficients)
repeat(thisDegreeLog2) {
val next = MutableList<C>(powersOf2[it + 1] * numeratorDegree + 1) { constantZero }
add(next)
val last = last()
multiplyAddingTo(
ring = ring,
multiplicand = last,
multiplicandDegree = powersOf2[it] * numeratorDegree + 1,
multiplier = last,
multiplierDegree = powersOf2[it] * numeratorDegree + 1,
target = next,
)
}
}
val denominatorPowers = buildList<List<C>>(thisDegreeLog2 + 1) {
add(arg.denominator.coefficients)
repeat(thisDegreeLog2) {
val next = MutableList<C>(powersOf2[it + 1] * denominatorDegree + 1) { constantZero }
add(next)
val last = last()
multiplyAddingTo(
ring = ring,
multiplicand = last,
multiplicandDegree = powersOf2[it] * denominatorDegree + 1,
multiplier = last,
multiplierDegree = powersOf2[it] * denominatorDegree + 1,
target = next,
)
}
}
val levelResultCoefsPool = buildList<MutableList<C>>(thisDegreeLog2 + 1) {
repeat(thisDegreeLog2 + 1) {
add(MutableList(hashes[it] * argDegree) { constantZero })
}
}
val edgedMultiplier = MutableList<C>(0) { TODO() }
val edgedMultiplierUpdater = MutableList<C>(0) { TODO() }
fun MutableList<C>.reset() {
for (i in indices) set(i, constantZero)
}
fun processLevel(level: Int, start: Int, end: Int) : List<C> {
val levelResultCoefs = levelResultCoefsPool[level + 1]
if (level == -1) {
levelResultCoefs[0] = coefficients[start]
} else {
levelResultCoefs.reset()
multiplyAddingTo(
ring = ring,
multiplicand = processLevel(level = level - 1, start = start, end = (start + end) / 2),
multiplicandDegree = hashes[level] * argDegree,
multiplier = denominatorPowers[level],
multiplierDegree = powersOf2[level] * denominatorDegree,
target = levelResultCoefs
)
multiplyAddingTo(
ring = ring,
multiplicand = processLevel(level = level - 1, start = (start + end) / 2, end = end),
multiplicandDegree = hashes[level] * argDegree,
multiplier = numeratorPowers[level],
multiplierDegree = powersOf2[level] * numeratorDegree,
target = levelResultCoefs
)
}
return levelResultCoefs
}
fun processLevelEdged(level: Int, start: Int, end: Int) : List<C> {
val levelResultCoefs = levelResultCoefsPool[level + 1]
if (level == -1) {
levelResultCoefs[0] = coefficients[start]
} else {
val levelsPowerOf2 = powersOf2[level]
if (end - start >= levelsPowerOf2) {
multiplyAddingTo(
ring = ring,
multiplicand = processLevelEdged(level = level - 1, start = start + levelsPowerOf2, end = end),
multiplicandDegree = hashes[level] * argDegree, // TODO: Ввести переменную
multiplier = numeratorPowers[level],
multiplierDegree = powersOf2[level] * numeratorDegree,
target = levelResultCoefs
)
multiplyAddingTo(
ring = ring,
multiplicand = processLevel(level = level - 1, start = start, end = start + levelsPowerOf2),
multiplicandDegree = hashes[level] * argDegree,
multiplier = edgedMultiplier,
multiplierDegree = max((hashes[level] and thisDegree) - powersOf2[level] + 1, 0) * denominatorDegree, // TODO: Ввести переменную
target = levelResultCoefs
)
if (level != thisDegreeLog2) {
multiplyAddingToUpdater(
ring = ring,
multiplicand = edgedMultiplier,
multiplicandDegree = max((hashes[level] and thisDegree) - powersOf2[level] + 1, 0) * denominatorDegree, // TODO: Ввести переменную
multiplier = denominatorPowers[level],
multiplierDegree = powersOf2[level] * denominatorDegree,
updater = edgedMultiplierUpdater,
zero = constantZero
)
}
} else {
copyTo(
origin = processLevelEdged(level = level - 1, start = start + levelsPowerOf2, end = end),
originDegree = hashes[level] * argDegree, // TODO: Ввести переменную
target = levelResultCoefs
)
}
}
return levelResultCoefs
}
return ListPolynomial(
processLevelEdged(
level = thisDegreeLog2,
start = 0,
end = thisDegree + 1
)
)
}

13
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/misc.kt Normal file
View File

@ -0,0 +1,13 @@
/*
* 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
@RequiresOptIn(
message = "It's copy of operation with optimized boxing. It's currently unstable.",
level = RequiresOptIn.Level.ERROR
)
internal annotation class UnstablePolynomialBoxingOptimization

491
kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/ListPolynomialTest.kt Normal file
View File

@ -0,0 +1,491 @@
/*
* 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.test.misc.*
import kotlin.test.*
class ListPolynomialTest {
@Test
fun test_Polynomial_Int_plus() {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(-22, 9), Rational(-8, 9), Rational(-8, 7)),
ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)) + -3,
"test 1"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)) + 2,
"test 2"
)
assertEquals(
ListPolynomial(Rational(0)),
ListPolynomial(Rational(-2)) + 2,
"test 3"
)
assertEquals(
ListPolynomial(),
ListPolynomial<Rational>() + 0,
"test 4"
)
assertEquals(
ListPolynomial(Rational(-1), Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)) + 1,
"test 5"
)
assertEquals(
ListPolynomial(Rational(-1)),
ListPolynomial(Rational(-2)) + 1,
"test 6"
)
assertEquals(
ListPolynomial(Rational(2)),
ListPolynomial<Rational>() + 2,
"test 7"
)
}
}
@Test
fun test_Polynomial_Int_minus() {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(32, 9), Rational(-8, 9), Rational(-8, 7)),
ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)) - -3,
"test 1"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(2), Rational(0), Rational(0), Rational(0)) - 2,
"test 2"
)
assertEquals(
ListPolynomial(Rational(0)),
ListPolynomial(Rational(2)) - 2,
"test 3"
)
assertEquals(
ListPolynomial(),
ListPolynomial<Rational>() - 0,
"test 4"
)
assertEquals(
ListPolynomial(Rational(1), Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(2), Rational(0), Rational(0), Rational(0)) - 1,
"test 5"
)
assertEquals(
ListPolynomial(Rational(1)),
ListPolynomial(Rational(2)) - 1,
"test 6"
)
assertEquals(
ListPolynomial(Rational(-2)),
ListPolynomial<Rational>() - 2,
"test 7"
)
}
}
@Test
fun test_Polynomial_Int_times() {
IntModuloRing(35).listPolynomialSpace {
assertEquals(
ListPolynomial(34, 2, 1, 20, 2),
ListPolynomial(22, 26, 13, 15, 26) * 27,
"test 1"
)
assertEquals(
ListPolynomial(0, 0, 0, 0, 0),
ListPolynomial(7, 0, 49, 21, 14) * 15,
"test 2"
)
}
}
@Test
fun test_Int_Polynomial_plus() {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(-22, 9), Rational(-8, 9), Rational(-8, 7)),
-3 + ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)),
"test 1"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
2 + ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
"test 2"
)
assertEquals(
ListPolynomial(Rational(0)),
2 + ListPolynomial(Rational(-2)),
"test 3"
)
assertEquals(
ListPolynomial(),
0 + ListPolynomial(),
"test 4"
)
assertEquals(
ListPolynomial(Rational(-1), Rational(0), Rational(0), Rational(0)),
1 + ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
"test 5"
)
assertEquals(
ListPolynomial(Rational(-1)),
1 + ListPolynomial(Rational(-2)),
"test 6"
)
assertEquals(
ListPolynomial(Rational(2)),
2 + ListPolynomial(),
"test 7"
)
}
}
@Test
fun test_Int_Polynomial_minus() {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(32, 9), Rational(-8, 9), Rational(-8, 7)),
3 - ListPolynomial(Rational(-5, 9), Rational(8, 9), Rational(8, 7)),
"test 1"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
-2 - ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
"test 2"
)
assertEquals(
ListPolynomial(Rational(0)),
-2 - ListPolynomial(Rational(-2)),
"test 3"
)
assertEquals(
ListPolynomial(),
0 - ListPolynomial(),
"test 4"
)
assertEquals(
ListPolynomial(Rational(1), Rational(0), Rational(0), Rational(0)),
-1 - ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
"test 5"
)
assertEquals(
ListPolynomial(Rational(1)),
-1 - ListPolynomial(Rational(-2)),
"test 6"
)
assertEquals(
ListPolynomial(Rational(-2)),
-2 - ListPolynomial(),
"test 7"
)
}
}
@Test
fun test_Int_Polynomial_times() {
IntModuloRing(35).listPolynomialSpace {
assertEquals(
ListPolynomial(34, 2, 1, 20, 2),
27 * ListPolynomial(22, 26, 13, 15, 26),
"test 1"
)
assertEquals(
ListPolynomial(0, 0, 0, 0, 0),
15 * ListPolynomial(7, 0, 49, 21, 14),
"test 2"
)
}
}
@Test
fun test_Polynomial_Constant_plus() {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(-22, 9), Rational(-8, 9), Rational(-8, 7)),
ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)) + Rational(-3),
"test 1"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)) + Rational(2),
"test 2"
)
assertEquals(
ListPolynomial(Rational(0)),
ListPolynomial(Rational(-2)) + Rational(2),
"test 3"
)
assertEquals(
ListPolynomial(Rational(0)),
ListPolynomial<Rational>() + Rational(0),
"test 4"
)
assertEquals(
ListPolynomial(Rational(-1), Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)) + Rational(1),
"test 5"
)
assertEquals(
ListPolynomial(Rational(-1)),
ListPolynomial(Rational(-2)) + Rational(1),
"test 6"
)
assertEquals(
ListPolynomial(Rational(2)),
ListPolynomial<Rational>() + Rational(2),
"test 7"
)
}
}
@Test
fun test_Polynomial_Constant_minus() {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(32, 9), Rational(-8, 9), Rational(-8, 7)),
ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)) - Rational(-3),
"test 1"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(2), Rational(0), Rational(0), Rational(0)) - Rational(2),
"test 2"
)
assertEquals(
ListPolynomial(Rational(0)),
ListPolynomial(Rational(2)) - Rational(2),
"test 3"
)
assertEquals(
ListPolynomial(Rational(0)),
ListPolynomial<Rational>() - Rational(0),
"test 4"
)
assertEquals(
ListPolynomial(Rational(1), Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(2), Rational(0), Rational(0), Rational(0)) - Rational(1),
"test 5"
)
assertEquals(
ListPolynomial(Rational(1)),
ListPolynomial(Rational(2)) - Rational(1),
"test 6"
)
assertEquals(
ListPolynomial(Rational(-2)),
ListPolynomial<Rational>() - Rational(2),
"test 7"
)
}
}
@Test
fun test_Polynomial_Constant_times() {
IntModuloRing(35).listPolynomialSpace {
assertEquals(
ListPolynomial(34, 2, 1, 20, 2),
ListPolynomial(22, 26, 13, 15, 26) * 27.asConstant(),
"test 1"
)
assertEquals(
ListPolynomial(0, 0, 0, 0, 0),
ListPolynomial(7, 0, 49, 21, 14) * 15.asConstant(),
"test 2"
)
}
}
@Test
fun test_Constant_Polynomial_plus() {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(-22, 9), Rational(-8, 9), Rational(-8, 7)),
Rational(-3) + ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)),
"test 1"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
Rational(2) + ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
"test 2"
)
assertEquals(
ListPolynomial(Rational(0)),
Rational(2) + ListPolynomial(Rational(-2)),
"test 3"
)
assertEquals(
ListPolynomial(Rational(0)),
Rational(0) + ListPolynomial(),
"test 4"
)
assertEquals(
ListPolynomial(Rational(-1), Rational(0), Rational(0), Rational(0)),
Rational(1) + ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
"test 5"
)
assertEquals(
ListPolynomial(Rational(-1)),
Rational(1) + ListPolynomial(Rational(-2)),
"test 6"
)
assertEquals(
ListPolynomial(Rational(2)),
Rational(2) + ListPolynomial(),
"test 7"
)
}
}
@Test
fun test_Constant_Polynomial_minus() {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(32, 9), Rational(-8, 9), Rational(-8, 7)),
Rational(3) - ListPolynomial(Rational(-5, 9), Rational(8, 9), Rational(8, 7)),
"test 1"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
Rational(-2) - ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
"test 2"
)
assertEquals(
ListPolynomial(Rational(0)),
Rational(-2) - ListPolynomial(Rational(-2)),
"test 3"
)
assertEquals(
ListPolynomial(Rational(0)),
Rational(0) - ListPolynomial(),
"test 4"
)
assertEquals(
ListPolynomial(Rational(1), Rational(0), Rational(0), Rational(0)),
Rational(-1) - ListPolynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
"test 5"
)
assertEquals(
ListPolynomial(Rational(1)),
Rational(-1) - ListPolynomial(Rational(-2)),
"test 6"
)
assertEquals(
ListPolynomial(Rational(-2)),
Rational(-2) - ListPolynomial(),
"test 7"
)
}
}
@Test
fun test_Constant_Polynomial_times() {
IntModuloRing(35).listPolynomialSpace {
assertEquals(
ListPolynomial(34, 2, 1, 20, 2),
27 * ListPolynomial(22, 26, 13, 15, 26),
"test 1"
)
assertEquals(
ListPolynomial(0, 0, 0, 0, 0),
15 * ListPolynomial(7, 0, 49, 21, 14),
"test 2"
)
}
}
@Test
fun test_Polynomial_unaryMinus() {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(-5, 9), Rational(8, 9), Rational(8, 7)),
-ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)),
"test 1"
)
assertEquals(
ListPolynomial(Rational(-5, 9), Rational(8, 9), Rational(8, 7), Rational(0), Rational(0)),
-ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7), Rational(0), Rational(0)),
"test 2"
)
}
}
@Test
fun test_Polynomial_Polynomial_plus() {
RationalField.listPolynomialSpace {
// (5/9 - 8/9 x - 8/7 x^2) + (-5/7 + 5/1 x + 5/8 x^2) ?= -10/63 + 37/9 x - 29/56 x^2
assertEquals(
ListPolynomial(Rational(-10, 63), Rational(37, 9), Rational(-29, 56)),
ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)) +
ListPolynomial(Rational(-5, 7), Rational(5, 1), Rational(5, 8)),
"test 1"
)
// (-2/9 - 8/3 x) + (0 + 9/4 x + 2/4 x^2) ?= -2/9 - 5/12 x + 2/4 x^2
assertEquals(
ListPolynomial(Rational(-2, 9), Rational(-5, 12), Rational(2, 4)),
ListPolynomial(Rational(-2, 9), Rational(-8, 3)) +
ListPolynomial(Rational(0), Rational(9, 4), Rational(2, 4)),
"test 2"
)
// (-4/7 - 2/6 x + 0 x^2 + 0 x^3) + (-6/3 - 7/2 x + 2/3 x^2) ?= -18/7 - 23/6 x + 2/3 x^2
assertEquals(
ListPolynomial(Rational(-18, 7), Rational(-23, 6), Rational(2, 3), Rational(0)),
ListPolynomial(Rational(-4, 7), Rational(-2, 6), Rational(0), Rational(0)) +
ListPolynomial(Rational(-6, 3), Rational(-7, 2), Rational(2, 3)),
"test 3"
)
// (-2/4 - 6/9 x - 4/9 x^2) + (2/4 + 6/9 x + 4/9 x^2) ?= 0
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(-2, 4), Rational(-6, 9), Rational(-4, 9)) +
ListPolynomial(Rational(2, 4), Rational(6, 9), Rational(4, 9)),
"test 4"
)
}
}
@Test
fun test_Polynomial_Polynomial_minus() {
RationalField.listPolynomialSpace {
// (5/9 - 8/9 x - 8/7 x^2) - (-5/7 + 5/1 x + 5/8 x^2) ?= 80/63 - 53/9 x - 99/56 x^2
assertEquals(
ListPolynomial(Rational(80, 63), Rational(-53, 9), Rational(-99, 56)),
ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)) -
ListPolynomial(Rational(-5, 7), Rational(5, 1), Rational(5, 8)),
"test 1"
)
// (-2/9 - 8/3 x) - (0 + 9/4 x + 2/4 x^2) ?= -2/9 - 59/12 x - 2/4 x^2
assertEquals(
ListPolynomial(Rational(-2, 9), Rational(-59, 12), Rational(-2, 4)),
ListPolynomial(Rational(-2, 9), Rational(-8, 3)) -
ListPolynomial(Rational(0), Rational(9, 4), Rational(2, 4)),
"test 2"
)
// (-4/7 - 2/6 x + 0 x^2 + 0 x^3) - (-6/3 - 7/2 x + 2/3 x^2) ?= 10/7 + 19/6 x - 2/3 x^2
assertEquals(
ListPolynomial(Rational(10, 7), Rational(19, 6), Rational(-2, 3), Rational(0)),
ListPolynomial(Rational(-4, 7), Rational(-2, 6), Rational(0), Rational(0)) -
ListPolynomial(Rational(-6, 3), Rational(-7, 2), Rational(2, 3)),
"test 3"
)
// (-2/4 - 6/9 x - 4/9 x^2) - (-2/4 - 6/9 x - 4/9 x^2) ?= 0
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(-2, 4), Rational(-6, 9), Rational(-4, 9)) -
ListPolynomial(Rational(-2, 4), Rational(-6, 9), Rational(-4, 9)),
"test 4"
)
}
}
@Test
fun test_Polynomial_Polynomial_times() {
IntModuloRing(35).listPolynomialSpace {
// (1 + x + x^2) * (1 - x + x^2) ?= 1 + x^2 + x^4
assertEquals(
ListPolynomial(1, 0, 1, 0, 1),
ListPolynomial(1, -1, 1) * ListPolynomial(1, 1, 1),
"test 1"
)
// Spoiler: 5 * 7 = 0
assertEquals(
ListPolynomial(0, 0, 0, 0, 0),
ListPolynomial(5, -25, 10) * ListPolynomial(21, 14, -7),
"test 2"
)
}
}
}

259
kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/ListPolynomialUtilTest.kt Normal file
View File

@ -0,0 +1,259 @@
/*
* 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.test.misc.Rational
import space.kscience.kmath.test.misc.RationalField
import kotlin.test.Test
import kotlin.test.assertEquals
import kotlin.test.assertFailsWith
@OptIn(UnstableKMathAPI::class)
class ListPolynomialUtilTest {
@Test
fun test_substitute_Double() {
assertEquals(
0.0,
ListPolynomial(1.0, -2.0, 1.0).substitute(1.0),
0.001,
"test 1"
)
assertEquals(
1.1931904761904761,
ListPolynomial(0.625, 2.6666666666666665, 0.5714285714285714, 1.5).substitute(0.2),
0.001,
"test 2"
)
assertEquals(
0.5681904761904762,
ListPolynomial(0.0, 2.6666666666666665, 0.5714285714285714, 1.5).substitute(0.2),
0.001,
"test 3"
)
assertEquals(
1.1811904761904761,
ListPolynomial(0.625, 2.6666666666666665, 0.5714285714285714, 0.0).substitute(0.2),
0.001,
"test 4"
)
assertEquals(
1.1703333333333332,
ListPolynomial(0.625, 2.6666666666666665, 0.0, 1.5).substitute(0.2),
0.001,
"test 5"
)
}
@Test
fun test_substitute_Constant() {
assertEquals(
Rational(0),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).substitute(RationalField, Rational(1)),
"test 1"
)
assertEquals(
Rational(25057, 21000),
ListPolynomial(Rational(5,8), Rational(8, 3), Rational(4, 7), Rational(3, 2))
.substitute(RationalField, Rational(1, 5)),
"test 2"
)
assertEquals(
Rational(2983, 5250),
ListPolynomial(Rational(0), Rational(8, 3), Rational(4, 7), Rational(3, 2))
.substitute(RationalField, Rational(1, 5)),
"test 3"
)
assertEquals(
Rational(4961, 4200),
ListPolynomial(Rational(5,8), Rational(8, 3), Rational(4, 7), Rational(0))
.substitute(RationalField, Rational(1, 5)),
"test 4"
)
assertEquals(
Rational(3511, 3000),
ListPolynomial(Rational(5,8), Rational(8, 3), Rational(0), Rational(3, 2))
.substitute(RationalField, Rational(1, 5)),
"test 5"
)
}
@Test
fun test_substitute_Polynomial() {
assertEquals(
ListPolynomial(Rational(0)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).substitute(RationalField, ListPolynomial(Rational(1))),
"test 1"
)
assertEquals(
ListPolynomial(Rational(709, 378), Rational(155, 252), Rational(19, 525), Rational(2, 875)),
ListPolynomial(Rational(1, 7), Rational(9, 4), Rational(1, 3), Rational(2, 7))
.substitute(RationalField, ListPolynomial(Rational(6, 9), Rational(1, 5))),
"test 2"
)
assertEquals(
ListPolynomial(Rational(655, 378), Rational(155, 252), Rational(19, 525), Rational(2, 875)),
ListPolynomial(Rational(0), Rational(9, 4), Rational(1, 3), Rational(2, 7))
.substitute(RationalField, ListPolynomial(Rational(6, 9), Rational(1, 5))),
"test 3"
)
assertEquals(
ListPolynomial(Rational(677, 378), Rational(97, 180), Rational(1, 75), Rational(0)),
ListPolynomial(Rational(1, 7), Rational(9, 4), Rational(1, 3), Rational(0))
.substitute(RationalField, ListPolynomial(Rational(6, 9), Rational(1, 5))),
"test 4"
)
assertEquals(
ListPolynomial(Rational(653, 378), Rational(221, 420), Rational(4, 175), Rational(2, 875)),
ListPolynomial(Rational(1, 7), Rational(9, 4), Rational(0), Rational(2, 7))
.substitute(RationalField, ListPolynomial(Rational(6, 9), Rational(1, 5))),
"test 5"
)
assertEquals(
ListPolynomial(Rational(89, 54), Rational(0), Rational(0), Rational(0)),
ListPolynomial(Rational(0), Rational(9, 4), Rational(1, 3), Rational(0))
.substitute(RationalField, ListPolynomial(Rational(6, 9), Rational(0))),
"test 6"
)
}
@Test
fun test_derivative() {
assertEquals(
ListPolynomial(Rational(-2), Rational(2)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).derivative(RationalField),
"test 1"
)
assertEquals(
ListPolynomial(Rational(-8, 3), Rational(8, 9), Rational(15, 7), Rational(-20, 9)),
ListPolynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).derivative(RationalField),
"test 2"
)
assertEquals(
ListPolynomial(Rational(0), Rational(8, 9), Rational(15, 7), Rational(-20, 9)),
ListPolynomial(Rational(0), Rational(0), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).derivative(RationalField),
"test 3"
)
assertEquals(
ListPolynomial(Rational(-8, 3), Rational(8, 9), Rational(15, 7), Rational(0)),
ListPolynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(0)).derivative(RationalField),
"test 4"
)
}
@Test
fun test_nthDerivative() {
assertEquals(
ListPolynomial(Rational(-2), Rational(2)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthDerivative(RationalField, 1),
"test 1"
)
assertFailsWith<IllegalArgumentException>("test2") {
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthDerivative(RationalField, -1)
}
assertEquals(
ListPolynomial(Rational(1), Rational(-2), Rational(1)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthDerivative(RationalField, 0),
"test 3"
)
assertEquals(
ListPolynomial(Rational(2)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthDerivative(RationalField, 2),
"test 4"
)
assertEquals(
ListPolynomial(),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthDerivative(RationalField, 3),
"test 5"
)
assertEquals(
ListPolynomial(),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthDerivative(RationalField, 4),
"test 6"
)
assertEquals(
ListPolynomial(Rational(8, 9), Rational(30, 7), Rational(-20, 3)),
ListPolynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).nthDerivative(RationalField, 2),
"test 7"
)
assertEquals(
ListPolynomial(Rational(8, 9), Rational(30, 7), Rational(-20, 3)),
ListPolynomial(Rational(0), Rational(0), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).nthDerivative(RationalField, 2),
"test 8"
)
assertEquals(
ListPolynomial(Rational(8, 9), Rational(30, 7), Rational(0)),
ListPolynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(0)).nthDerivative(RationalField, 2),
"test 9"
)
}
@Test
fun test_antiderivative() {
assertEquals(
ListPolynomial(Rational(0), Rational(1), Rational(-1), Rational(1, 3)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).antiderivative(RationalField),
"test 1"
)
assertEquals(
ListPolynomial(Rational(0), Rational(1, 5), Rational(-4, 3), Rational(4, 27), Rational(5, 28), Rational(-1, 9)),
ListPolynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).antiderivative(RationalField),
"test 2"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(4, 27), Rational(5, 28), Rational(-1, 9)),
ListPolynomial(Rational(0), Rational(0), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).antiderivative(RationalField),
"test 3"
)
assertEquals(
ListPolynomial(Rational(0), Rational(1, 5), Rational(-4, 3), Rational(4, 27), Rational(5, 28), Rational(0)),
ListPolynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(0)).antiderivative(RationalField),
"test 4"
)
}
@Test
fun test_nthAntiderivative() {
assertEquals(
ListPolynomial(Rational(0), Rational(1), Rational(-1), Rational(1, 3)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthAntiderivative(RationalField, 1),
"test 1"
)
assertFailsWith<IllegalArgumentException>("test2") {
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthAntiderivative(RationalField, -1)
}
assertEquals(
ListPolynomial(Rational(1), Rational(-2), Rational(1)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthAntiderivative(RationalField, 0),
"test 3"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(1, 2), Rational(-1, 3), Rational(1, 12)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthAntiderivative(RationalField, 2),
"test 4"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(1, 6), Rational(-1, 12), Rational(1, 60)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthAntiderivative(RationalField, 3),
"test 5"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(0), Rational(1, 24), Rational(-1, 60), Rational(1, 360)),
ListPolynomial(Rational(1), Rational(-2), Rational(1)).nthAntiderivative(RationalField, 4),
"test 6"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(1, 10), Rational(-4, 9), Rational(1, 27), Rational(1, 28), Rational(-1, 54)),
ListPolynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).nthAntiderivative(RationalField, 2),
"test 7"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(0), Rational(0), Rational(1, 27), Rational(1, 28), Rational(-1, 54)),
ListPolynomial(Rational(0), Rational(0), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).nthAntiderivative(RationalField, 2),
"test 8"
)
assertEquals(
ListPolynomial(Rational(0), Rational(0), Rational(1, 10), Rational(-4, 9), Rational(1, 27), Rational(1, 28), Rational(0)),
ListPolynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(0)).nthAntiderivative(RationalField, 2),
"test 9"
)
}
}

138
kmath-functions/src/commonTest/kotlin/space/kscience/kmath/test/misc/IntModulo.kt Normal file
View File

@ -0,0 +1,138 @@
/*
* 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.test.misc
import space.kscience.kmath.functions.ListPolynomial
import space.kscience.kmath.functions.ListPolynomialSpace
import space.kscience.kmath.operations.Ring
class IntModulo {
val residue: Int
val modulus: Int
@PublishedApi
internal constructor(residue: Int, modulus: Int, toCheckInput: Boolean = true) {
if (toCheckInput) {
require(modulus != 0) { "modulus can not be zero" }
this.modulus = if (modulus < 0) -modulus else modulus
this.residue = residue.mod(modulus)
} else {
this.residue = residue
this.modulus = modulus
}
}
constructor(residue: Int, modulus: Int) : this(residue, modulus, true)
operator fun unaryPlus(): IntModulo = this
operator fun unaryMinus(): IntModulo =
IntModulo(
if (residue == 0) 0 else modulus - residue,
modulus,
toCheckInput = false
)
operator fun plus(other: IntModulo): IntModulo {
require(modulus == other.modulus) { "can not add two residue different modulo" }
return IntModulo(
(residue + other.residue) % modulus,
modulus,
toCheckInput = false
)
}
operator fun plus(other: Int): IntModulo =
IntModulo(
(residue + other) % modulus,
modulus,
toCheckInput = false
)
operator fun minus(other: IntModulo): IntModulo {
require(modulus == other.modulus) { "can not subtract two residue different modulo" }
return IntModulo(
(residue - other.residue) % modulus,
modulus,
toCheckInput = false
)
}
operator fun minus(other: Int): IntModulo =
IntModulo(
(residue - other) % modulus,
modulus,
toCheckInput = false
)
operator fun times(other: IntModulo): IntModulo {
require(modulus == other.modulus) { "can not multiply two residue different modulo" }
return IntModulo(
(residue * other.residue) % modulus,
modulus,
toCheckInput = false
)
}
operator fun times(other: Int): IntModulo =
IntModulo(
(residue * other) % modulus,
modulus,
toCheckInput = false
)
operator fun div(other: IntModulo): IntModulo {
require(modulus == other.modulus) { "can not divide two residue different modulo" }
val (reciprocalCandidate, gcdOfOtherResidueAndModulus) = bezoutIdentityWithGCD(other.residue, modulus)
require(gcdOfOtherResidueAndModulus == 1) { "can not divide to residue that has non-trivial GCD with modulo" }
return IntModulo(
(residue * reciprocalCandidate) % modulus,
modulus,
toCheckInput = false
)
}
operator fun div(other: Int): IntModulo {
val (reciprocalCandidate, gcdOfOtherResidueAndModulus) = bezoutIdentityWithGCD(other, modulus)
require(gcdOfOtherResidueAndModulus == 1) { "can not divide to residue that has non-trivial GCD with modulo" }
return IntModulo(
(residue * reciprocalCandidate) % modulus,
modulus,
toCheckInput = false
)
}
override fun equals(other: Any?): Boolean =
when (other) {
is IntModulo -> residue == other.residue && modulus == other.modulus
else -> false
}
override fun hashCode(): Int = residue.hashCode()
override fun toString(): String = "$residue mod $modulus"
}
@Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
class IntModuloRing : Ring<IntModulo> {
val modulus: Int
constructor(modulus: Int) {
require(modulus != 0) { "modulus can not be zero" }
this.modulus = if (modulus < 0) -modulus else modulus
}
override inline val zero: IntModulo get() = IntModulo(0, modulus, toCheckInput = false)
override inline val one: IntModulo get() = IntModulo(1, modulus, toCheckInput = false)
fun number(arg: Int) = IntModulo(arg, modulus, toCheckInput = false)
override inline fun add(left: IntModulo, right: IntModulo): IntModulo = left + right
override inline fun multiply(left: IntModulo, right: IntModulo): IntModulo = left * right
override inline fun IntModulo.unaryMinus(): IntModulo = -this
override inline fun IntModulo.plus(arg: IntModulo): IntModulo = this + arg
override inline fun IntModulo.minus(arg: IntModulo): IntModulo = this - arg
override inline fun IntModulo.times(arg: IntModulo): IntModulo = this * arg
inline fun IntModulo.div(arg: IntModulo): IntModulo = this / arg
}
fun ListPolynomialSpace<IntModulo, IntModuloRing>.ListPolynomial(vararg coefs: Int): ListPolynomial<IntModulo> =
ListPolynomial(coefs.map { IntModulo(it, ring.modulus) })
fun IntModuloRing.ListPolynomial(vararg coefs: Int): ListPolynomial<IntModulo> =
ListPolynomial(coefs.map { IntModulo(it, modulus) })

29
kmath-functions/src/commonTest/kotlin/space/kscience/kmath/test/misc/misc.kt Normal file
View File

@ -0,0 +1,29 @@
/*
* 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.test.misc
import kotlin.math.abs
data class BezoutIdentityWithGCD<T>(val first: T, val second: T, val gcd: T)
tailrec fun gcd(a: Long, b: Long): Long = if (a == 0L) abs(b) else gcd(b % a, a)
fun bezoutIdentityWithGCD(a: Int, b: Int): BezoutIdentityWithGCD<Int> =
when {
a < 0 && b < 0 -> with(bezoutIdentityWithGCDInternalLogic(-a, -b, 1, 0, 0, 1)) { BezoutIdentityWithGCD(-first, -second, gcd) }
a < 0 -> with(bezoutIdentityWithGCDInternalLogic(-a, b, 1, 0, 0, 1)) { BezoutIdentityWithGCD(-first, second, gcd) }
b < 0 -> with(bezoutIdentityWithGCDInternalLogic(a, -b, 1, 0, 0, 1)) { BezoutIdentityWithGCD(first, -second, gcd) }
else -> bezoutIdentityWithGCDInternalLogic(a, b, 1, 0, 0, 1)
}
internal tailrec fun bezoutIdentityWithGCDInternalLogic(a: Int, b: Int, m1: Int, m2: Int, m3: Int, m4: Int): BezoutIdentityWithGCD<Int> =
if (b == 0) BezoutIdentityWithGCD(m1, m3, a)
else {
val quotient = a / b
val reminder = a % b
bezoutIdentityWithGCDInternalLogic(b, reminder, m2, m1 - quotient * m2, m4, m3 - quotient * m4)
}