Feature: Polynomials and rational functions #469
@ -220,7 +220,7 @@ public open class PolynomialSpace<C, A : Ring<C>>(
|
||||
public override operator fun C.minus(other: Polynomial<C>): Polynomial<C> =
|
||||
if (this.isZero()) other
|
||||
else with(other.coefficients) {
|
||||
if (isEmpty()) Polynomial(listOf(-this@minus))
|
||||
if (isEmpty()) Polynomial(listOf(this@minus))
|
||||
else Polynomial(
|
||||
toMutableList()
|
||||
.apply {
|
||||
|
@ -6,11 +6,11 @@
|
||||
package space.kscience.kmath.functions
|
||||
|
||||
import space.kscience.kmath.operations.algebra
|
||||
import space.kscience.kmath.test.misc.Rational
|
||||
import space.kscience.kmath.test.misc.RationalField
|
||||
import space.kscience.kmath.test.misc.*
|
||||
import kotlin.test.Test
|
||||
import kotlin.test.assertEquals
|
||||
|
||||
|
||||
class PolynomialTest {
|
||||
@Test
|
||||
fun test_Polynomial_Int_plus() {
|
||||
@ -93,6 +93,116 @@ class PolynomialTest {
|
||||
}
|
||||
}
|
||||
@Test
|
||||
fun test_Polynomial_Int_times() {
|
||||
IntModuloRing(35).polynomial {
|
||||
assertEquals(
|
||||
Polynomial(34, 2, 1, 20, 2),
|
||||
Polynomial(22, 26, 13, 15, 26) * 27,
|
||||
"test 1"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
Polynomial(7, 0, 49, 21, 14) * 15,
|
||||
"test 2"
|
||||
)
|
||||
}
|
||||
}
|
||||
@Test
|
||||
fun test_Int_Polynomial_plus() {
|
||||
RationalField.polynomial {
|
||||
assertEquals(
|
||||
Polynomial(Rational(-22, 9), Rational(-8, 9), Rational(-8, 7)),
|
||||
-3 + Polynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)),
|
||||
"test 1"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
|
||||
2 + Polynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
|
||||
"test 2"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
2 + Polynomial(Rational(-2)),
|
||||
"test 3"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
0 + Polynomial(),
|
||||
"test 4"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(-1), Rational(0), Rational(0), Rational(0)),
|
||||
1 + Polynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
|
||||
"test 5"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(-1)),
|
||||
1 + Polynomial(Rational(-2)),
|
||||
"test 6"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(2)),
|
||||
2 + Polynomial(),
|
||||
"test 7"
|
||||
)
|
||||
}
|
||||
}
|
||||
@Test
|
||||
fun test_Int_Polynomial_minus() {
|
||||
RationalField.polynomial {
|
||||
assertEquals(
|
||||
Polynomial(Rational(32, 9), Rational(-8, 9), Rational(-8, 7)),
|
||||
3 - Polynomial(Rational(-5, 9), Rational(8, 9), Rational(8, 7)),
|
||||
"test 1"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
|
||||
-2 - Polynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
|
||||
"test 2"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
-2 - Polynomial(Rational(-2)),
|
||||
"test 3"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
0 - Polynomial(),
|
||||
"test 4"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(1), Rational(0), Rational(0), Rational(0)),
|
||||
-1 - Polynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
|
||||
"test 5"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(1)),
|
||||
-1 - Polynomial(Rational(-2)),
|
||||
"test 6"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(-2)),
|
||||
-2 - Polynomial(),
|
||||
"test 7"
|
||||
)
|
||||
}
|
||||
}
|
||||
@Test
|
||||
fun test_Int_Polynomial_times() {
|
||||
IntModuloRing(35).polynomial {
|
||||
assertEquals(
|
||||
Polynomial(34, 2, 1, 20, 2),
|
||||
27 * Polynomial(22, 26, 13, 15, 26),
|
||||
"test 1"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
15 * Polynomial(7, 0, 49, 21, 14),
|
||||
"test 2"
|
||||
)
|
||||
}
|
||||
}
|
||||
@Test
|
||||
fun test_Polynomial_Constant_plus() {
|
||||
RationalField.polynomial {
|
||||
assertEquals(
|
||||
@ -173,6 +283,116 @@ class PolynomialTest {
|
||||
}
|
||||
}
|
||||
@Test
|
||||
fun test_Polynomial_Constant_times() {
|
||||
IntModuloRing(35).polynomial {
|
||||
assertEquals(
|
||||
Polynomial(34, 2, 1, 20, 2),
|
||||
Polynomial(22, 26, 13, 15, 26) * number(27),
|
||||
"test 1"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
Polynomial(7, 0, 49, 21, 14) * number(15),
|
||||
"test 2"
|
||||
)
|
||||
}
|
||||
}
|
||||
@Test
|
||||
fun test_Constant_Polynomial_plus() {
|
||||
RationalField.polynomial {
|
||||
assertEquals(
|
||||
Polynomial(Rational(-22, 9), Rational(-8, 9), Rational(-8, 7)),
|
||||
Rational(-3) + Polynomial(Rational(5, 9), Rational(-8, 9), Rational(-8, 7)),
|
||||
"test 1"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
|
||||
Rational(2) + Polynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
|
||||
"test 2"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
Rational(2) + Polynomial(Rational(-2)),
|
||||
"test 3"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
Rational(0) + Polynomial(),
|
||||
"test 4"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(-1), Rational(0), Rational(0), Rational(0)),
|
||||
Rational(1) + Polynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
|
||||
"test 5"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(-1)),
|
||||
Rational(1) + Polynomial(Rational(-2)),
|
||||
"test 6"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(2)),
|
||||
Rational(2) + Polynomial(),
|
||||
"test 7"
|
||||
)
|
||||
}
|
||||
}
|
||||
@Test
|
||||
fun test_Constant_Polynomial_minus() {
|
||||
RationalField.polynomial {
|
||||
assertEquals(
|
||||
Polynomial(Rational(32, 9), Rational(-8, 9), Rational(-8, 7)),
|
||||
Rational(3) - Polynomial(Rational(-5, 9), Rational(8, 9), Rational(8, 7)),
|
||||
"test 1"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(0), Rational(0), Rational(0), Rational(0)),
|
||||
Rational(-2) - Polynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
|
||||
"test 2"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
Rational(-2) - Polynomial(Rational(-2)),
|
||||
"test 3"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
Rational(0) - Polynomial(),
|
||||
"test 4"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(1), Rational(0), Rational(0), Rational(0)),
|
||||
Rational(-1) - Polynomial(Rational(-2), Rational(0), Rational(0), Rational(0)),
|
||||
"test 5"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(1)),
|
||||
Rational(-1) - Polynomial(Rational(-2)),
|
||||
"test 6"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(Rational(-2)),
|
||||
Rational(-2) - Polynomial(),
|
||||
"test 7"
|
||||
)
|
||||
}
|
||||
}
|
||||
@Test
|
||||
fun test_Constant_Polynomial_times() {
|
||||
IntModuloRing(35).polynomial {
|
||||
assertEquals(
|
||||
Polynomial(34, 2, 1, 20, 2),
|
||||
27 * Polynomial(22, 26, 13, 15, 26),
|
||||
"test 1"
|
||||
)
|
||||
assertEquals(
|
||||
Polynomial(),
|
||||
15 * Polynomial(7, 0, 49, 21, 14),
|
||||
"test 2"
|
||||
)
|
||||
}
|
||||
}
|
||||
@Test
|
||||
fun test_Polynomial_Polynomial_plus() {
|
||||
RationalField.polynomial {
|
||||
// (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
|
||||
|
@ -0,0 +1,142 @@
|
||||
/*
|
||||
* 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.Polynomial
|
||||
import space.kscience.kmath.functions.PolynomialSpace
|
||||
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(0, 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 PolynomialSpace<IntModulo, IntModuloRing>.number(arg: Int) = IntModulo(arg, ring.modulus, toCheckInput = false)
|
||||
|
||||
fun PolynomialSpace<IntModulo, IntModuloRing>.Polynomial(vararg coefs: Int): Polynomial<IntModulo> =
|
||||
Polynomial(coefs.map { IntModulo(it, ring.modulus) })
|
||||
fun IntModuloRing.Polynomial(vararg coefs: Int): Polynomial<IntModulo> =
|
||||
Polynomial(coefs.map { IntModulo(it, modulus) })
|
@ -114,10 +114,6 @@ class Rational {
|
||||
override fun toString(): String = if (denominator == 1L) "$numerator" else "$numerator/$denominator"
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Algebraic structure for rational numbers.
|
||||
*/
|
||||
@Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
|
||||
@OptIn(UnstableKMathAPI::class)
|
||||
object RationalField : Field<Rational>, NumbersAddOps<Rational> {
|
||||
|
@ -5,8 +5,27 @@
|
||||
|
||||
package space.kscience.kmath.test.misc
|
||||
|
||||
import kotlin.math.abs
|
||||
|
||||
// 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