diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/ListPolynomial.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/ListPolynomial.kt new file mode 100644 index 000000000..fce179fc8 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/ListPolynomial.kt @@ -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( + /** + * 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 +) : Polynomial { + 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>( + public override val ring: A, +) : PolynomialSpaceOverRing, 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.plus(other: Int): ListPolynomial = + 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.minus(other: Int): ListPolynomial = + 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.times(other: Int): ListPolynomial = + 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): ListPolynomial = + 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): ListPolynomial = + 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): ListPolynomial = + 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 = number(constantNumber(value)) + + /** + * Returns sum of the constant represented as a polynomial and the polynomial. + */ + public override operator fun C.plus(other: ListPolynomial): ListPolynomial = + 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): ListPolynomial = + 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): ListPolynomial = + 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.plus(other: C): ListPolynomial = + 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.minus(other: C): ListPolynomial = + 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.times(other: C): ListPolynomial = + 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 = ListPolynomial(listOf(value)) + + /** + * Returns negation of the polynomial. + */ + public override operator fun ListPolynomial.unaryMinus(): ListPolynomial = + ListPolynomial(coefficients.map { -it }) + /** + * Returns sum of the polynomials. + */ + public override operator fun ListPolynomial.plus(other: ListPolynomial): ListPolynomial { + 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.minus(other: ListPolynomial): ListPolynomial { + 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.times(other: ListPolynomial): ListPolynomial { + 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 = ListPolynomial(emptyList()) + /** + * Instance of unit polynomial (unit of the polynomial ring). + */ + override val one: ListPolynomial 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.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.substitute(argument: C): C = this.substitute(ring, argument) + /** + * Substitutes provided polynomial [argument] into [this] polynomial. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListPolynomial.substitute(argument: ListPolynomial): ListPolynomial = this.substitute(ring, argument) + + /** + * Represent [this] polynomial as a regular context-less function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListPolynomial.asFunction(): (C) -> C = { this.substitute(ring, it) } + /** + * Represent [this] polynomial as a regular context-less function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListPolynomial.asFunctionOfConstant(): (C) -> C = { this.substitute(ring, it) } + /** + * Represent [this] polynomial as a regular context-less function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListPolynomial.asFunctionOfPolynomial(): (ListPolynomial) -> ListPolynomial = { this.substitute(ring, it) } + + /** + * Evaluates value of [this] polynomial on provided [argument]. + */ + @Suppress("NOTHING_TO_INLINE") + public inline operator fun ListPolynomial.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.invoke(argument: ListPolynomial): ListPolynomial = 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( + ring: A, +) : ListPolynomialSpace(ring), ScaleOperations> where A : Ring, A : ScaleOperations { + override fun scale(a: ListPolynomial, value: Double): ListPolynomial = + ring { ListPolynomial(a.coefficients.map { scale(it, value) }) } +} diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/ListRationalFunction.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/ListRationalFunction.kt new file mode 100644 index 000000000..f3e352bcd --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/ListRationalFunction.kt @@ -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( + public override val numerator: ListPolynomial, + public override val denominator: ListPolynomial +) : RationalFunction> { + 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> ( + public val ring: A, +) : + RationalFunctionalSpaceOverPolynomialSpace< + C, + ListPolynomial, + ListRationalFunction, + ListPolynomialSpace, + >, + PolynomialSpaceOfFractions< + C, + ListPolynomial, + ListRationalFunction, + >() { + + /** + * Underlying polynomial ring. Its polynomial operations are inherited by local polynomial operations. + */ + override val polynomialRing : ListPolynomialSpace = ListPolynomialSpace(ring) + /** + * Constructor of [ListRationalFunction] from numerator and denominator [ListPolynomial]. + */ + override fun constructRationalFunction(numerator: ListPolynomial, denominator: ListPolynomial): ListRationalFunction = + 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.substitute(argument: C): C = this.substitute(ring, argument) + /** + * Substitutes provided polynomial [argument] into [this] polynomial. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListPolynomial.substitute(argument: ListPolynomial): ListPolynomial = this.substitute(ring, argument) + /** + * Substitutes provided rational function [argument] into [this] polynomial. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListPolynomial.substitute(argument: ListRationalFunction): ListRationalFunction = this.substitute(ring, argument) + /** + * Substitutes provided polynomial [argument] into [this] rational function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListRationalFunction.substitute(argument: ListPolynomial): ListRationalFunction = this.substitute(ring, argument) + /** + * Substitutes provided rational function [argument] into [this] rational function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListRationalFunction.substitute(argument: ListRationalFunction): ListRationalFunction = this.substitute(ring, argument) + + /** + * Represent [this] polynomial as a regular context-less function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListPolynomial.asFunction(): (C) -> C = { this.substitute(ring, it) } + /** + * Represent [this] polynomial as a regular context-less function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListPolynomial.asFunctionOfConstant(): (C) -> C = { this.substitute(ring, it) } + /** + * Represent [this] polynomial as a regular context-less function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListPolynomial.asFunctionOfPolynomial(): (ListPolynomial) -> ListPolynomial = { this.substitute(ring, it) } + /** + * Represent [this] polynomial as a regular context-less function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListPolynomial.asFunctionOfRationalFunction(): (ListRationalFunction) -> ListRationalFunction = { this.substitute(ring, it) } + /** + * Represent [this] rational function as a regular context-less function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListRationalFunction.asFunctionOfPolynomial(): (ListPolynomial) -> ListRationalFunction = { this.substitute(ring, it) } + /** + * Represent [this] rational function as a regular context-less function. + */ + @Suppress("NOTHING_TO_INLINE") + public inline fun ListRationalFunction.asFunctionOfRationalFunction(): (ListRationalFunction) -> ListRationalFunction = { this.substitute(ring, it) } + + /** + * Evaluates value of [this] polynomial on provided argument. + */ + @Suppress("NOTHING_TO_INLINE") + public inline operator fun ListPolynomial.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.invoke(argument: ListPolynomial): ListPolynomial = this.substitute(ring, argument) + /** + * Evaluates value of [this] polynomial on provided argument. + */ + @Suppress("NOTHING_TO_INLINE") + public inline operator fun ListPolynomial.invoke(argument: ListRationalFunction): ListRationalFunction = this.substitute(ring, argument) + /** + * Evaluates value of [this] rational function on provided argument. + */ + @Suppress("NOTHING_TO_INLINE") + public inline operator fun ListRationalFunction.invoke(argument: ListPolynomial): ListRationalFunction = this.substitute(ring, argument) + /** + * Evaluates value of [this] rational function on provided argument. + */ + @Suppress("NOTHING_TO_INLINE") + public inline operator fun ListRationalFunction.invoke(argument: ListRationalFunction): ListRationalFunction = this.substitute(ring, argument) +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt index f2eba10d5..12490d133 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt @@ -435,10 +435,12 @@ public interface MultivariatePolynomialSpace>: 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) /** diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/RationalFunction.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/RationalFunction.kt index edc9dfa5c..01911f980 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/RationalFunction.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/RationalFunction.kt @@ -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() } /** diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/listUtil.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/listUtil.kt new file mode 100644 index 000000000..127dd8c7a --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/listUtil.kt @@ -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 > A.listPolynomialSpace(): ListPolynomialSpace = + 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 , R> A.listPolynomialSpace(block: ListPolynomialSpace.() -> R): R { + contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) } + return ListPolynomialSpace(this).block() +} + +/** + * Creates a [ScalableListPolynomialSpace] over a received scalable ring. + */ +public fun A.scalableListPolynomialSpace(): ScalableListPolynomialSpace where A : Ring, A : ScaleOperations = + 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 A.scalableListPolynomialSpace(block: ScalableListPolynomialSpace.() -> R): R where A : Ring, A : ScaleOperations { + contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) } + return ScalableListPolynomialSpace(this).block() +} + +/** + * Creates a [ListRationalFunctionSpace] over a received ring. + */ +public fun > A.listRationalFunctionSpace(): ListRationalFunctionSpace = + 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 , R> A.listRationalFunctionSpace(block: ListRationalFunctionSpace.() -> 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.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 ListPolynomial.substitute(ring: Ring, 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 ListPolynomial.substitute(ring: Ring, arg: ListPolynomial) : ListPolynomial = + ring.listPolynomialSpace { + if (coefficients.isEmpty()) return zero + var result: ListPolynomial = 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 ListPolynomial.substitute(ring: Ring, arg: ListRationalFunction) : ListRationalFunction = + ring.listRationalFunctionSpace { + if (coefficients.isEmpty()) return zero + var result: ListRationalFunction = 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.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 ListRationalFunction.substitute(ring: Field, 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 ListRationalFunction.substitute(ring: Ring, arg: ListPolynomial) : ListRationalFunction = + 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 ListRationalFunction.substitute(ring: Ring, arg: ListRationalFunction) : ListRationalFunction = + ring.listRationalFunctionSpace { + numerator.substitute(ring, arg) / denominator.substitute(ring, arg) + } + +/** + * Represent [this] polynomial as a regular context-less function. + */ +public fun > ListPolynomial.asFunctionOver(ring: A): (C) -> C = { substitute(ring, it) } + +/** + * Represent [this] polynomial as a regular context-less function. + */ +public fun > ListPolynomial.asPolynomialFunctionOver(ring: A): (ListPolynomial) -> ListPolynomial = { substitute(ring, it) } + +/** + * Represent [this] polynomial as a regular context-less function. + */ +public fun > ListPolynomial.asFunctionOfRationalFunctionOver(ring: A): (ListPolynomial) -> ListPolynomial = { substitute(ring, it) } + +/** + * Represent [this] rational function as a regular context-less function. + */ +public fun > ListRationalFunction.asFunctionOver(ring: A): (C) -> C = { substitute(ring, it) } + +/** + * Represent [this] rational function as a regular context-less function. + */ +public fun > ListRationalFunction.asPolynomialFunctionOver(ring: A): (ListPolynomial) -> ListRationalFunction = { substitute(ring, it) } + +/** + * Represent [this] rational function as a regular context-less function. + */ +public fun > ListRationalFunction.asFunctionOfRationalFunctionOver(ring: A): (ListPolynomial) -> ListRationalFunction = { substitute(ring, it) } + +/** + * Returns algebraic derivative of received polynomial. + */ +@UnstableKMathAPI +public fun ListPolynomial.derivative( + ring: A, +): ListPolynomial where A : Ring, A : NumericAlgebra = 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 ListPolynomial.nthDerivative( + ring: A, + order: Int, +): ListPolynomial where A : Ring, A : NumericAlgebra = 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 ListPolynomial.antiderivative( + ring: A, +): ListPolynomial where A : Field, A : NumericAlgebra = 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 ListPolynomial.nthAntiderivative( + ring: A, + order: Int, +): ListPolynomial where A : Field, A : NumericAlgebra = 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 > ListPolynomial.integrate( + ring: Field, + range: ClosedRange, +): 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 ListRationalFunction.derivative( + ring: A, +): ListRationalFunction where A : Ring, A : NumericAlgebra = 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 ListRationalFunction.nthDerivative( + ring: A, + order: Int, +): ListRationalFunction where A : Ring, A : NumericAlgebra = + if (order == 0) this else derivative(ring).nthDerivative(ring, order - 1) \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/listUtilOptimized.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/listUtilOptimized.kt new file mode 100644 index 000000000..6eb3a1dc7 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/listUtilOptimized.kt @@ -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 copyTo( + origin: List, + originDegree: Int, + target: MutableList, +) { + for (deg in 0 .. originDegree) target[deg] = origin[deg] +} + +@UnstablePolynomialBoxingOptimization +@Suppress("NOTHING_TO_INLINE") +internal inline fun multiplyAddingToUpdater( + ring: Ring, + multiplicand: MutableList, + multiplicandDegree: Int, + multiplier: List, + multiplierDegree: Int, + updater: MutableList, + 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 multiplyAddingTo( + ring: Ring, + multiplicand: List, + multiplicandDegree: Int, + multiplier: List, + multiplierDegree: Int, + target: MutableList +) = 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 ListPolynomial.substitute2(ring: Ring, arg: ListPolynomial) : ListPolynomial = 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 = MutableList(thisDegree * argDegree + 1) { constantZero } + resultCoefs[0] = coefficients[thisDegree] + val resultCoefsUpdate: MutableList = 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(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 ListPolynomial.substituteRationalFunctionTakeNumerator(ring: Ring, arg: ListRationalFunction): ListPolynomial = 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(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>(thisDegreeLog2 + 1) { + add(arg.numerator.coefficients) + repeat(thisDegreeLog2) { + val next = MutableList(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>(thisDegreeLog2 + 1) { + add(arg.denominator.coefficients) + repeat(thisDegreeLog2) { + val next = MutableList(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>(thisDegreeLog2 + 1) { + repeat(thisDegreeLog2 + 1) { + add(MutableList(hashes[it] * argDegree) { constantZero }) + } + } + val edgedMultiplier = MutableList(0) { TODO() } + val edgedMultiplierUpdater = MutableList(0) { TODO() } + + fun MutableList.reset() { + for (i in indices) set(i, constantZero) + } + + fun processLevel(level: Int, start: Int, end: Int) : List { + 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 { + 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 + ) + ) +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/misc.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/misc.kt new file mode 100644 index 000000000..8b6fac39e --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/misc.kt @@ -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 \ No newline at end of file diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/ListPolynomialTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/ListPolynomialTest.kt new file mode 100644 index 000000000..c9950fac5 --- /dev/null +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/ListPolynomialTest.kt @@ -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() + 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() + 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() - 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() - 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(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(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(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(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" + ) + } + } +} \ No newline at end of file diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/ListPolynomialUtilTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/ListPolynomialUtilTest.kt new file mode 100644 index 000000000..69c1611f3 --- /dev/null +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/ListPolynomialUtilTest.kt @@ -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("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("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" + ) + } +} \ No newline at end of file diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/test/misc/IntModulo.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/test/misc/IntModulo.kt new file mode 100644 index 000000000..89764db46 --- /dev/null +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/test/misc/IntModulo.kt @@ -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 { + + 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.ListPolynomial(vararg coefs: Int): ListPolynomial = + ListPolynomial(coefs.map { IntModulo(it, ring.modulus) }) +fun IntModuloRing.ListPolynomial(vararg coefs: Int): ListPolynomial = + ListPolynomial(coefs.map { IntModulo(it, modulus) }) \ No newline at end of file diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/test/misc/misc.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/test/misc/misc.kt new file mode 100644 index 000000000..ed41b9245 --- /dev/null +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/test/misc/misc.kt @@ -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(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 = + 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 = + 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) + } \ No newline at end of file