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

Feature: Polynomials and rational functions #469

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

32
kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/ListPolynomialTest.kt
View File

@ -12,7 +12,7 @@ import kotlin.test.*
class ListPolynomialTest {
@Test
fun test_Polynomial_Int_plus() {
RationalField.listPolynomial {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(-22, 9), Rational(-8, 9), Rational(-8, 7)),
ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)) + -3,
@ -52,7 +52,7 @@ class ListPolynomialTest {
}
@Test
fun test_Polynomial_Int_minus() {
RationalField.listPolynomial {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(32, 9), Rational(-8, 9), Rational(-8, 7)),
ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)) - -3,
@ -92,7 +92,7 @@ class ListPolynomialTest {
}
@Test
fun test_Polynomial_Int_times() {
IntModuloRing(35).listPolynomial {
IntModuloRing(35).listPolynomialSpace {
assertEquals(
ListPolynomial(34, 2, 1, 20, 2),
ListPolynomial(22, 26, 13, 15, 26) * 27,
@ -107,7 +107,7 @@ class ListPolynomialTest {
}
@Test
fun test_Int_Polynomial_plus() {
RationalField.listPolynomial {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(-22, 9), Rational(-8, 9), Rational(-8, 7)),
-3 + ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)),
@ -147,7 +147,7 @@ class ListPolynomialTest {
}
@Test
fun test_Int_Polynomial_minus() {
RationalField.listPolynomial {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(32, 9), Rational(-8, 9), Rational(-8, 7)),
3 - ListPolynomial(Rational(-5, 9), Rational(8, 9), Rational(8, 7)),
@ -187,7 +187,7 @@ class ListPolynomialTest {
}
@Test
fun test_Int_Polynomial_times() {
IntModuloRing(35).listPolynomial {
IntModuloRing(35).listPolynomialSpace {
assertEquals(
ListPolynomial(34, 2, 1, 20, 2),
27 * ListPolynomial(22, 26, 13, 15, 26),
@ -202,7 +202,7 @@ class ListPolynomialTest {
}
@Test
fun test_Polynomial_Constant_plus() {
RationalField.listPolynomial {
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),
@ -242,7 +242,7 @@ class ListPolynomialTest {
}
@Test
fun test_Polynomial_Constant_minus() {
RationalField.listPolynomial {
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),
@ -282,7 +282,7 @@ class ListPolynomialTest {
}
@Test
fun test_Polynomial_Constant_times() {
IntModuloRing(35).listPolynomial {
IntModuloRing(35).listPolynomialSpace {
assertEquals(
ListPolynomial(34, 2, 1, 20, 2),
ListPolynomial(22, 26, 13, 15, 26) * 27.asConstant(),
@ -297,7 +297,7 @@ class ListPolynomialTest {
}
@Test
fun test_Constant_Polynomial_plus() {
RationalField.listPolynomial {
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)),
@ -337,7 +337,7 @@ class ListPolynomialTest {
}
@Test
fun test_Constant_Polynomial_minus() {
RationalField.listPolynomial {
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)),
@ -377,7 +377,7 @@ class ListPolynomialTest {
}
@Test
fun test_Constant_Polynomial_times() {
IntModuloRing(35).listPolynomial {
IntModuloRing(35).listPolynomialSpace {
assertEquals(
ListPolynomial(34, 2, 1, 20, 2),
27 * ListPolynomial(22, 26, 13, 15, 26),
@ -392,7 +392,7 @@ class ListPolynomialTest {
}
@Test
fun test_Polynomial_unaryMinus() {
RationalField.listPolynomial {
RationalField.listPolynomialSpace {
assertEquals(
ListPolynomial(Rational(-5, 9), Rational(8, 9), Rational(8, 7)),
-ListPolynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)),
@ -407,7 +407,7 @@ class ListPolynomialTest {
}
@Test
fun test_Polynomial_Polynomial_plus() {
RationalField.listPolynomial {
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)),
@ -440,7 +440,7 @@ class ListPolynomialTest {
}
@Test
fun test_Polynomial_Polynomial_minus() {
RationalField.listPolynomial {
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)),
@ -473,7 +473,7 @@ class ListPolynomialTest {
}
@Test
fun test_Polynomial_Polynomial_times() {
IntModuloRing(35).listPolynomial {
IntModuloRing(35).listPolynomialSpace {
// (1 + x + x^2) * (1 - x + x^2) ?= 1 + x^2 + x^4
assertEquals(
ListPolynomial(1, 0, 1, 0, 1),

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

@ -0,0 +1,140 @@
/*
* 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>.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) })

20
kmath-functions/src/commonTest/kotlin/space/kscience/kmath/test/misc/misc.kt
View File

@ -8,4 +8,22 @@ package space.kscience.kmath.test.misc
import kotlin.math.abs
tailrec fun gcd(a: Long, b: Long): Long = if (a == 0L) abs(b) else gcd(b % a, a)
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)
}