Feature: Polynomials and rational functions #469
@ -1,92 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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.
|
|
||||||
*/
|
|
||||||
|
|
||||||
@file:OptIn(UnstableKMathAPI::class)
|
|
||||||
|
|
||||||
package space.kscience.kmath.interpolation
|
|
||||||
|
|
||||||
import space.kscience.kmath.data.XYColumnarData
|
|
||||||
import space.kscience.kmath.functions.PiecewisePolynomial
|
|
||||||
import space.kscience.kmath.functions.asFunction
|
|
||||||
import space.kscience.kmath.functions.substitute
|
|
||||||
import space.kscience.kmath.misc.UnstableKMathAPI
|
|
||||||
import space.kscience.kmath.operations.Ring
|
|
||||||
import space.kscience.kmath.structures.Buffer
|
|
||||||
import space.kscience.kmath.structures.asBuffer
|
|
||||||
|
|
||||||
/**
|
|
||||||
* And interpolator for data with x column type [X], y column type [Y].
|
|
||||||
*/
|
|
||||||
public fun interface Interpolator<T, in X : T, Y : T> {
|
|
||||||
public fun interpolate(points: XYColumnarData<T, X, Y>): (X) -> Y
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* And interpolator returning [PiecewisePolynomial] function
|
|
||||||
*/
|
|
||||||
public interface PolynomialInterpolator<T : Comparable<T>> : Interpolator<T, T, T> {
|
|
||||||
public val algebra: Ring<T>
|
|
||||||
|
|
||||||
public fun getDefaultValue(): T = error("Out of bounds")
|
|
||||||
|
|
||||||
public fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T>
|
|
||||||
|
|
||||||
override fun interpolate(points: XYColumnarData<T, T, T>): (T) -> T = { x ->
|
|
||||||
interpolatePolynomials(points).substitute(algebra, x) ?: getDefaultValue()
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
|
|
||||||
x: Buffer<T>,
|
|
||||||
y: Buffer<T>,
|
|
||||||
): PiecewisePolynomial<T> {
|
|
||||||
val pointSet = XYColumnarData.of(x, y)
|
|
||||||
return interpolatePolynomials(pointSet)
|
|
||||||
}
|
|
||||||
|
|
||||||
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
|
|
||||||
data: Map<T, T>,
|
|
||||||
): PiecewisePolynomial<T> {
|
|
||||||
val pointSet = XYColumnarData.of(data.keys.toList().asBuffer(), data.values.toList().asBuffer())
|
|
||||||
return interpolatePolynomials(pointSet)
|
|
||||||
}
|
|
||||||
|
|
||||||
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
|
|
||||||
data: List<Pair<T, T>>,
|
|
||||||
): PiecewisePolynomial<T> {
|
|
||||||
val pointSet = XYColumnarData.of(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer())
|
|
||||||
return interpolatePolynomials(pointSet)
|
|
||||||
}
|
|
||||||
|
|
||||||
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
|
|
||||||
x: Buffer<T>,
|
|
||||||
y: Buffer<T>,
|
|
||||||
): (T) -> T? = interpolatePolynomials(x, y).asFunction(algebra)
|
|
||||||
|
|
||||||
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
|
|
||||||
data: Map<T, T>,
|
|
||||||
): (T) -> T? = interpolatePolynomials(data).asFunction(algebra)
|
|
||||||
|
|
||||||
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
|
|
||||||
data: List<Pair<T, T>>,
|
|
||||||
): (T) -> T? = interpolatePolynomials(data).asFunction(algebra)
|
|
||||||
|
|
||||||
|
|
||||||
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
|
|
||||||
x: Buffer<T>,
|
|
||||||
y: Buffer<T>,
|
|
||||||
defaultValue: T,
|
|
||||||
): (T) -> T = interpolatePolynomials(x, y).asFunction(algebra, defaultValue)
|
|
||||||
|
|
||||||
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
|
|
||||||
data: Map<T, T>,
|
|
||||||
defaultValue: T,
|
|
||||||
): (T) -> T = interpolatePolynomials(data).asFunction(algebra, defaultValue)
|
|
||||||
|
|
||||||
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
|
|
||||||
data: List<Pair<T, T>>,
|
|
||||||
defaultValue: T,
|
|
||||||
): (T) -> T = interpolatePolynomials(data).asFunction(algebra, defaultValue)
|
|
||||||
@ -1,43 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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.interpolation
|
|
||||||
|
|
||||||
import space.kscience.kmath.data.XYColumnarData
|
|
||||||
import space.kscience.kmath.functions.PiecewisePolynomial
|
|
||||||
import space.kscience.kmath.functions.ListPolynomial
|
|
||||||
import space.kscience.kmath.misc.UnstableKMathAPI
|
|
||||||
import space.kscience.kmath.operations.Field
|
|
||||||
import space.kscience.kmath.operations.invoke
|
|
||||||
|
|
||||||
@OptIn(UnstableKMathAPI::class)
|
|
||||||
internal fun <T : Comparable<T>> insureSorted(points: XYColumnarData<*, T, *>) {
|
|
||||||
for (i in 0 until points.size - 1)
|
|
||||||
require(points.x[i + 1] > points.x[i]) { "Input data is not sorted at index $i" }
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Reference JVM implementation: https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/analysis/interpolation/LinearInterpolator.java
|
|
||||||
*/
|
|
||||||
public class LinearInterpolator<T : Comparable<T>>(override val algebra: Field<T>) : PolynomialInterpolator<T> {
|
|
||||||
|
|
||||||
@OptIn(UnstableKMathAPI::class)
|
|
||||||
override fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T> = algebra {
|
|
||||||
require(points.size > 0) { "Point array should not be empty" }
|
|
||||||
insureSorted(points)
|
|
||||||
|
|
||||||
PiecewisePolynomial(points.x[0]) {
|
|
||||||
for (i in 0 until points.size - 1) {
|
|
||||||
val slope = (points.y[i + 1] - points.y[i]) / (points.x[i + 1] - points.x[i])
|
|
||||||
val const = points.y[i] - slope * points.x[i]
|
|
||||||
val polynomial = ListPolynomial(const, slope)
|
|
||||||
putRight(points.x[i + 1], polynomial)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
public val <T : Comparable<T>> Field<T>.linearInterpolator: LinearInterpolator<T>
|
|
||||||
get() = LinearInterpolator(this)
|
|
||||||
@ -1,83 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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.interpolation
|
|
||||||
|
|
||||||
import space.kscience.kmath.data.XYColumnarData
|
|
||||||
import space.kscience.kmath.functions.PiecewisePolynomial
|
|
||||||
import space.kscience.kmath.functions.ListPolynomial
|
|
||||||
import space.kscience.kmath.misc.UnstableKMathAPI
|
|
||||||
import space.kscience.kmath.operations.DoubleField
|
|
||||||
import space.kscience.kmath.operations.Field
|
|
||||||
import space.kscience.kmath.operations.invoke
|
|
||||||
import space.kscience.kmath.structures.DoubleBuffer
|
|
||||||
import space.kscience.kmath.structures.MutableBufferFactory
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Generic spline interpolator. Not recommended for performance critical places, use platform-specific and type
|
|
||||||
* specific ones.
|
|
||||||
*
|
|
||||||
* Based on
|
|
||||||
* https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/SplineInterpolator.java
|
|
||||||
*/
|
|
||||||
public class SplineInterpolator<T : Comparable<T>>(
|
|
||||||
override val algebra: Field<T>,
|
|
||||||
public val bufferFactory: MutableBufferFactory<T>,
|
|
||||||
) : PolynomialInterpolator<T> {
|
|
||||||
//TODO possibly optimize zeroed buffers
|
|
||||||
|
|
||||||
@OptIn(UnstableKMathAPI::class)
|
|
||||||
override fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T> = algebra {
|
|
||||||
require(points.size >= 3) { "Can't use spline interpolator with less than 3 points" }
|
|
||||||
insureSorted(points)
|
|
||||||
// Number of intervals. The number of data points is n + 1.
|
|
||||||
val n = points.size - 1
|
|
||||||
// Differences between knot points
|
|
||||||
val h = bufferFactory(n) { i -> points.x[i + 1] - points.x[i] }
|
|
||||||
val mu = bufferFactory(n) { zero }
|
|
||||||
val z = bufferFactory(n + 1) { zero }
|
|
||||||
|
|
||||||
for (i in 1 until n) {
|
|
||||||
val g = 2.0 * (points.x[i + 1] - points.x[i - 1]) - h[i - 1] * mu[i - 1]
|
|
||||||
mu[i] = h[i] / g
|
|
||||||
z[i] =
|
|
||||||
((points.y[i + 1] * h[i - 1] - points.y[i] * (points.x[i + 1] - points.x[i - 1]) + points.y[i - 1] * h[i]) * 3.0 /
|
|
||||||
(h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g
|
|
||||||
}
|
|
||||||
|
|
||||||
// cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)
|
|
||||||
|
|
||||||
PiecewisePolynomial(points.x[points.size - 1]) {
|
|
||||||
var cOld = zero
|
|
||||||
|
|
||||||
for (j in n - 1 downTo 0) {
|
|
||||||
val c = z[j] - mu[j] * cOld
|
|
||||||
val a = points.y[j]
|
|
||||||
val b = (points.y[j + 1] - points.y[j]) / h[j] - h[j] * (cOld + 2.0 * c) / 3.0
|
|
||||||
val d = (cOld - c) / (3.0 * h[j])
|
|
||||||
val x0 = points.x[j]
|
|
||||||
val x02 = x0 * x0
|
|
||||||
val x03 = x02 * x0
|
|
||||||
//Shift coefficients to represent absolute polynomial instead of one with an offset
|
|
||||||
val polynomial = ListPolynomial(
|
|
||||||
a - b * x0 + c * x02 - d * x03,
|
|
||||||
b - 2 * c * x0 + 3 * d * x02,
|
|
||||||
c - 3 * d * x0,
|
|
||||||
d
|
|
||||||
)
|
|
||||||
cOld = c
|
|
||||||
putLeft(x0, polynomial)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
public fun <T : Comparable<T>> Field<T>.splineInterpolator(
|
|
||||||
bufferFactory: MutableBufferFactory<T>,
|
|
||||||
): SplineInterpolator<T> = SplineInterpolator(this, bufferFactory)
|
|
||||||
|
|
||||||
public val DoubleField.splineInterpolator: SplineInterpolator<Double>
|
|
||||||
get() = SplineInterpolator(this, ::DoubleBuffer)
|
|
||||||
@ -1,48 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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.integration
|
|
||||||
|
|
||||||
import space.kscience.kmath.functions.ListPolynomial
|
|
||||||
import space.kscience.kmath.functions.integrate
|
|
||||||
import space.kscience.kmath.misc.UnstableKMathAPI
|
|
||||||
import space.kscience.kmath.operations.DoubleField
|
|
||||||
import kotlin.math.PI
|
|
||||||
import kotlin.math.sin
|
|
||||||
import kotlin.test.Test
|
|
||||||
import kotlin.test.assertEquals
|
|
||||||
|
|
||||||
@OptIn(UnstableKMathAPI::class)
|
|
||||||
class SplineIntegralTest {
|
|
||||||
|
|
||||||
@Test
|
|
||||||
fun integratePolynomial(){
|
|
||||||
val polynomial = ListPolynomial(1.0, 2.0, 3.0)
|
|
||||||
val integral = polynomial.integrate(DoubleField,1.0..2.0)
|
|
||||||
assertEquals(11.0, integral, 0.001)
|
|
||||||
}
|
|
||||||
|
|
||||||
@Test
|
|
||||||
fun gaussSin() {
|
|
||||||
val res = DoubleField.splineIntegrator.integrate(0.0..2 * PI, IntegrandMaxCalls(5)) { x ->
|
|
||||||
sin(x)
|
|
||||||
}
|
|
||||||
assertEquals(0.0, res.value, 1e-2)
|
|
||||||
}
|
|
||||||
|
|
||||||
@Test
|
|
||||||
fun gaussUniform() {
|
|
||||||
val res = DoubleField.splineIntegrator.integrate(35.0..100.0, IntegrandMaxCalls(20)) { x ->
|
|
||||||
if(x in 30.0..50.0){
|
|
||||||
1.0
|
|
||||||
} else {
|
|
||||||
0.0
|
|
||||||
}
|
|
||||||
}
|
|
||||||
assertEquals(15.0, res.value, 0.5)
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
}
|
|
||||||
@ -1,29 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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.interpolation
|
|
||||||
|
|
||||||
import space.kscience.kmath.operations.DoubleField
|
|
||||||
import kotlin.test.Test
|
|
||||||
import kotlin.test.assertEquals
|
|
||||||
|
|
||||||
internal class LinearInterpolatorTest {
|
|
||||||
@Test
|
|
||||||
fun testInterpolation() {
|
|
||||||
val data = listOf(
|
|
||||||
0.0 to 0.0,
|
|
||||||
1.0 to 1.0,
|
|
||||||
2.0 to 3.0,
|
|
||||||
3.0 to 4.0
|
|
||||||
)
|
|
||||||
|
|
||||||
//val polynomial: PiecewisePolynomial<Double> = DoubleField.linearInterpolator.interpolatePolynomials(data)
|
|
||||||
val function = DoubleField.linearInterpolator.interpolate(data)
|
|
||||||
assertEquals(null, function(-1.0))
|
|
||||||
assertEquals(0.5, function(0.5))
|
|
||||||
assertEquals(2.0, function(1.5))
|
|
||||||
assertEquals(3.0, function(2.0))
|
|
||||||
}
|
|
||||||
}
|
|
||||||
@ -1,31 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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.interpolation
|
|
||||||
|
|
||||||
import space.kscience.kmath.operations.DoubleField
|
|
||||||
import kotlin.math.PI
|
|
||||||
import kotlin.math.sin
|
|
||||||
import kotlin.test.Test
|
|
||||||
import kotlin.test.assertEquals
|
|
||||||
|
|
||||||
internal class SplineInterpolatorTest {
|
|
||||||
@Test
|
|
||||||
fun testInterpolation() {
|
|
||||||
val data = (0..10).map {
|
|
||||||
val x = it.toDouble() / 5 * PI
|
|
||||||
x to sin(x)
|
|
||||||
}
|
|
||||||
|
|
||||||
//val polynomial: PiecewisePolynomial<Double> = DoubleField.splineInterpolator.interpolatePolynomials(data)
|
|
||||||
|
|
||||||
val function = DoubleField.splineInterpolator.interpolate(data, Double.NaN)
|
|
||||||
|
|
||||||
assertEquals(Double.NaN, function(-1.0))
|
|
||||||
assertEquals(sin(0.5), function(0.5), 0.1)
|
|
||||||
assertEquals(sin(1.5), function(1.5), 0.1)
|
|
||||||
assertEquals(sin(2.0), function(2.0), 0.1)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
@ -1,142 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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.misc.UnstableKMathAPI
|
|
||||||
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")
|
|
||||||
@OptIn(UnstableKMathAPI::class)
|
|
||||||
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>.number(arg: Int) = IntModulo(arg, ring.modulus, toCheckInput = false)
|
|
||||||
|
|
||||||
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) })
|
|
||||||
@ -1,135 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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.misc.UnstableKMathAPI
|
|
||||||
import space.kscience.kmath.operations.Field
|
|
||||||
import space.kscience.kmath.operations.NumbersAddOps
|
|
||||||
|
|
||||||
class Rational {
|
|
||||||
companion object {
|
|
||||||
val ZERO: Rational = Rational(0L)
|
|
||||||
val ONE: Rational = Rational(1L)
|
|
||||||
}
|
|
||||||
|
|
||||||
val numerator: Long
|
|
||||||
val denominator: Long
|
|
||||||
|
|
||||||
internal constructor(numerator: Long, denominator: Long, toCheckInput: Boolean = true) {
|
|
||||||
if (toCheckInput) {
|
|
||||||
if (denominator == 0L) throw ArithmeticException("/ by zero")
|
|
||||||
|
|
||||||
val greatestCommonDivider = gcd(numerator, denominator).let { if (denominator < 0L) -it else it }
|
|
||||||
|
|
||||||
this.numerator = numerator / greatestCommonDivider
|
|
||||||
this.denominator = denominator / greatestCommonDivider
|
|
||||||
} else {
|
|
||||||
this.numerator = numerator
|
|
||||||
this.denominator = denominator
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
constructor(numerator: Int, denominator: Int) : this(numerator.toLong(), denominator.toLong(), true)
|
|
||||||
constructor(numerator: Int, denominator: Long) : this(numerator.toLong(), denominator, true)
|
|
||||||
constructor(numerator: Long, denominator: Int) : this(numerator, denominator.toLong(), true)
|
|
||||||
constructor(numerator: Long, denominator: Long) : this(numerator, denominator, true)
|
|
||||||
constructor(numerator: Int) : this(numerator.toLong(), 1L, false)
|
|
||||||
constructor(numerator: Long) : this(numerator, 1L, false)
|
|
||||||
|
|
||||||
operator fun unaryPlus(): Rational = this
|
|
||||||
operator fun unaryMinus(): Rational = Rational(-this.numerator, this.denominator)
|
|
||||||
operator fun plus(other: Rational): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator * other.denominator + denominator * other.numerator,
|
|
||||||
denominator * other.denominator
|
|
||||||
)
|
|
||||||
operator fun plus(other: Int): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator + denominator * other.toLong(),
|
|
||||||
denominator
|
|
||||||
)
|
|
||||||
operator fun plus(other: Long): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator + denominator * other,
|
|
||||||
denominator
|
|
||||||
)
|
|
||||||
operator fun minus(other: Rational): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator * other.denominator - denominator * other.numerator,
|
|
||||||
denominator * other.denominator
|
|
||||||
)
|
|
||||||
operator fun minus(other: Int): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator - denominator * other.toLong(),
|
|
||||||
denominator
|
|
||||||
)
|
|
||||||
operator fun minus(other: Long): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator - denominator * other,
|
|
||||||
denominator
|
|
||||||
)
|
|
||||||
operator fun times(other: Rational): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator * other.numerator,
|
|
||||||
denominator * other.denominator
|
|
||||||
)
|
|
||||||
operator fun times(other: Int): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator * other.toLong(),
|
|
||||||
denominator
|
|
||||||
)
|
|
||||||
operator fun times(other: Long): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator * other,
|
|
||||||
denominator
|
|
||||||
)
|
|
||||||
operator fun div(other: Rational): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator * other.denominator,
|
|
||||||
denominator * other.numerator
|
|
||||||
)
|
|
||||||
operator fun div(other: Int): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator,
|
|
||||||
denominator * other.toLong()
|
|
||||||
)
|
|
||||||
operator fun div(other: Long): Rational =
|
|
||||||
Rational(
|
|
||||||
numerator,
|
|
||||||
denominator * other
|
|
||||||
)
|
|
||||||
override fun equals(other: Any?): Boolean =
|
|
||||||
when (other) {
|
|
||||||
is Rational -> numerator == other.numerator && denominator == other.denominator
|
|
||||||
is Int -> numerator == other && denominator == 1L
|
|
||||||
is Long -> numerator == other && denominator == 1L
|
|
||||||
else -> false
|
|
||||||
}
|
|
||||||
|
|
||||||
override fun hashCode(): Int = 31 * numerator.hashCode() + denominator.hashCode()
|
|
||||||
|
|
||||||
override fun toString(): String = if (denominator == 1L) "$numerator" else "$numerator/$denominator"
|
|
||||||
}
|
|
||||||
|
|
||||||
@Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
|
|
||||||
@OptIn(UnstableKMathAPI::class)
|
|
||||||
object RationalField : Field<Rational>, NumbersAddOps<Rational> {
|
|
||||||
override inline val zero: Rational get() = Rational.ZERO
|
|
||||||
override inline val one: Rational get() = Rational.ONE
|
|
||||||
|
|
||||||
override inline fun number(value: Number): Rational = Rational(value.toLong())
|
|
||||||
|
|
||||||
override inline fun add(left: Rational, right: Rational): Rational = left + right
|
|
||||||
override inline fun multiply(left: Rational, right: Rational): Rational = left * right
|
|
||||||
override inline fun divide(left: Rational, right: Rational): Rational = left / right
|
|
||||||
override inline fun scale(a: Rational, value: Double): Rational = a * number(value)
|
|
||||||
|
|
||||||
override inline fun Rational.unaryMinus(): Rational = -this
|
|
||||||
override inline fun Rational.plus(arg: Rational): Rational = this + arg
|
|
||||||
override inline fun Rational.minus(arg: Rational): Rational = this - arg
|
|
||||||
override inline fun Rational.times(arg: Rational): Rational = this * arg
|
|
||||||
override inline fun Rational.div(arg: Rational): Rational = this / arg
|
|
||||||
}
|
|
||||||
@ -1,107 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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.misc.UnstableKMathAPI
|
|
||||||
import space.kscience.kmath.operations.*
|
|
||||||
|
|
||||||
class RationalWithMemorization private constructor(
|
|
||||||
val value: Rational,
|
|
||||||
override val memory : OperationsMemory
|
|
||||||
): WithMemorization {
|
|
||||||
companion object {
|
|
||||||
/**
|
|
||||||
* Constant containing the zero (the additive identity) of the [Rational] field.
|
|
||||||
*/
|
|
||||||
public val ZERO: RationalWithMemorization = RationalWithMemorization(Rational.ZERO, object : Endpoint {})
|
|
||||||
/**
|
|
||||||
* Constant containing the one (the multiplicative identity) of the [Rational] field.
|
|
||||||
*/
|
|
||||||
public val ONE: RationalWithMemorization = RationalWithMemorization(Rational.ONE, object : Endpoint {})
|
|
||||||
}
|
|
||||||
constructor(numerator: Int, denominator: Int) : this(Rational(numerator, denominator), object : Endpoint {})
|
|
||||||
constructor(numerator: Int, denominator: Long) : this(Rational(numerator, denominator), object : Endpoint {})
|
|
||||||
constructor(numerator: Long, denominator: Int) : this(Rational(numerator, denominator), object : Endpoint {})
|
|
||||||
constructor(numerator: Long, denominator: Long) : this(Rational(numerator, denominator), object : Endpoint {})
|
|
||||||
constructor(numerator: Int) : this(Rational(numerator), object : Endpoint {})
|
|
||||||
constructor(numerator: Long) : this(Rational(numerator), object : Endpoint {})
|
|
||||||
|
|
||||||
operator fun unaryPlus(): RationalWithMemorization = this
|
|
||||||
operator fun unaryMinus(): RationalWithMemorization = RationalWithMemorization(
|
|
||||||
-value,
|
|
||||||
object : Negation {
|
|
||||||
override val negated: OperationsMemory = memory
|
|
||||||
}
|
|
||||||
)
|
|
||||||
operator fun plus(other: RationalWithMemorization): RationalWithMemorization = RationalWithMemorization(
|
|
||||||
value + other.value,
|
|
||||||
object : Sum {
|
|
||||||
override val augend: OperationsMemory = memory
|
|
||||||
override val addend: OperationsMemory = other.memory
|
|
||||||
}
|
|
||||||
)
|
|
||||||
operator fun minus(other: RationalWithMemorization): RationalWithMemorization = RationalWithMemorization(
|
|
||||||
value - other.value,
|
|
||||||
object : Difference {
|
|
||||||
override val minuend: OperationsMemory = memory
|
|
||||||
override val subtrahend: OperationsMemory = other.memory
|
|
||||||
}
|
|
||||||
)
|
|
||||||
operator fun times(other: RationalWithMemorization): RationalWithMemorization = RationalWithMemorization(
|
|
||||||
value * other.value,
|
|
||||||
object : Product {
|
|
||||||
override val multiplicand: OperationsMemory = memory
|
|
||||||
override val multiplier: OperationsMemory = other.memory
|
|
||||||
}
|
|
||||||
)
|
|
||||||
operator fun div(other: RationalWithMemorization): RationalWithMemorization = RationalWithMemorization(
|
|
||||||
value / other.value,
|
|
||||||
object : Quotient {
|
|
||||||
override val dividend: OperationsMemory = memory
|
|
||||||
override val divisor: OperationsMemory = other.memory
|
|
||||||
}
|
|
||||||
)
|
|
||||||
|
|
||||||
override fun equals(other: Any?): Boolean =
|
|
||||||
other is RationalWithMemorization && value == other.value
|
|
||||||
|
|
||||||
override fun hashCode(): Int = value.hashCode()
|
|
||||||
|
|
||||||
override fun toString(): String = value.toString()
|
|
||||||
}
|
|
||||||
|
|
||||||
@Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
|
|
||||||
object RationalWithMemorizationRing : Ring<RationalWithMemorization> {
|
|
||||||
override inline val zero: RationalWithMemorization get() = RationalWithMemorization.ZERO
|
|
||||||
override inline val one: RationalWithMemorization get() = RationalWithMemorization.ONE
|
|
||||||
|
|
||||||
override inline fun add(left: RationalWithMemorization, right: RationalWithMemorization): RationalWithMemorization = left + right
|
|
||||||
override inline fun multiply(left: RationalWithMemorization, right: RationalWithMemorization): RationalWithMemorization = left * right
|
|
||||||
|
|
||||||
override inline fun RationalWithMemorization.unaryMinus(): RationalWithMemorization = -this
|
|
||||||
override inline fun RationalWithMemorization.plus(arg: RationalWithMemorization): RationalWithMemorization = this + arg
|
|
||||||
override inline fun RationalWithMemorization.minus(arg: RationalWithMemorization): RationalWithMemorization = this - arg
|
|
||||||
override inline fun RationalWithMemorization.times(arg: RationalWithMemorization): RationalWithMemorization = this * arg
|
|
||||||
}
|
|
||||||
|
|
||||||
@Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
|
|
||||||
object RationalWithMemorizationField : Field<RationalWithMemorization> {
|
|
||||||
override inline val zero: RationalWithMemorization get() = RationalWithMemorization.ZERO
|
|
||||||
override inline val one: RationalWithMemorization get() = RationalWithMemorization.ONE
|
|
||||||
|
|
||||||
override inline fun number(value: Number): RationalWithMemorization = RationalWithMemorization(value.toLong())
|
|
||||||
|
|
||||||
override inline fun add(left: RationalWithMemorization, right: RationalWithMemorization): RationalWithMemorization = left + right
|
|
||||||
override inline fun multiply(left: RationalWithMemorization, right: RationalWithMemorization): RationalWithMemorization = left * right
|
|
||||||
override inline fun divide(left: RationalWithMemorization, right: RationalWithMemorization): RationalWithMemorization = left / right
|
|
||||||
override inline fun scale(a: RationalWithMemorization, value: Double): RationalWithMemorization = a * number(value)
|
|
||||||
|
|
||||||
override inline fun RationalWithMemorization.unaryMinus(): RationalWithMemorization = -this
|
|
||||||
override inline fun RationalWithMemorization.plus(arg: RationalWithMemorization): RationalWithMemorization = this + arg
|
|
||||||
override inline fun RationalWithMemorization.minus(arg: RationalWithMemorization): RationalWithMemorization = this - arg
|
|
||||||
override inline fun RationalWithMemorization.times(arg: RationalWithMemorization): RationalWithMemorization = this * arg
|
|
||||||
override inline fun RationalWithMemorization.div(arg: RationalWithMemorization): RationalWithMemorization = this / arg
|
|
||||||
}
|
|
||||||
@ -1,51 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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
|
|
||||||
|
|
||||||
sealed interface OperationsMemory
|
|
||||||
|
|
||||||
interface Endpoint: OperationsMemory
|
|
||||||
|
|
||||||
interface Negation: OperationsMemory {
|
|
||||||
val negated: OperationsMemory
|
|
||||||
}
|
|
||||||
|
|
||||||
interface Sum: OperationsMemory {
|
|
||||||
val augend: OperationsMemory
|
|
||||||
val addend: OperationsMemory
|
|
||||||
}
|
|
||||||
|
|
||||||
interface Difference: OperationsMemory {
|
|
||||||
val minuend: OperationsMemory
|
|
||||||
val subtrahend: OperationsMemory
|
|
||||||
}
|
|
||||||
|
|
||||||
interface Product: OperationsMemory {
|
|
||||||
val multiplicand: OperationsMemory
|
|
||||||
val multiplier: OperationsMemory
|
|
||||||
}
|
|
||||||
|
|
||||||
interface Quotient: OperationsMemory {
|
|
||||||
val dividend: OperationsMemory
|
|
||||||
val divisor: OperationsMemory
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
fun equalMemories(one: OperationsMemory, other: OperationsMemory) : Boolean =
|
|
||||||
when(one) {
|
|
||||||
is Negation -> other is Negation && equalMemories(one.negated, other.negated)
|
|
||||||
is Sum -> other is Sum && equalMemories(one.augend, other.augend) && equalMemories(one.addend, other.addend)
|
|
||||||
is Difference -> other is Difference && equalMemories(one.minuend, other.minuend) && equalMemories(one.subtrahend, other.subtrahend)
|
|
||||||
is Product -> other is Product && equalMemories(one.multiplicand, other.multiplicand) && equalMemories(one.multiplier, other.multiplier)
|
|
||||||
is Quotient -> other is Quotient && equalMemories(one.dividend, other.dividend) && equalMemories(one.divisor, other.divisor)
|
|
||||||
is Endpoint -> one === other
|
|
||||||
}
|
|
||||||
|
|
||||||
interface WithMemorization {
|
|
||||||
val memory: OperationsMemory
|
|
||||||
}
|
|
||||||
|
|
||||||
fun equalMemories(one: WithMemorization, other: WithMemorization) : Boolean = equalMemories(one.memory, other.memory)
|
|
||||||
@ -1,31 +0,0 @@
|
|||||||
/*
|
|
||||||
* 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
|
|
||||||
|
|
||||||
// TODO: Move to corresponding module "kmath-number-theory"
|
|
||||||
|
|
||||||
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)
|
|
||||||
}
|
|
||||||
Loading…
Reference in New Issue
Block a user