Feature: Polynomials and rational functions #469

Merged
lounres merged 132 commits from feature/polynomials into dev 2022-07-28 18:04:06 +03:00
6 changed files with 45 additions and 72 deletions
Showing only changes of commit 5bc627f1d4 - Show all commits

View File

@ -117,16 +117,16 @@ public fun <T : Comparable<T>> PiecewisePolynomial(
* Return a value of polynomial function with given [ring] a given [arg] or null if argument is outside piecewise * Return a value of polynomial function with given [ring] a given [arg] or null if argument is outside piecewise
* definition. * definition.
*/ */
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.substitute(ring: C, arg: T): T? = public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.value(ring: C, arg: T): T? =
findPiece(arg)?.substitute(ring, arg) findPiece(arg)?.value(ring, arg)
/** /**
* Convert this polynomial to a function returning nullable value (null if argument is outside piecewise range). * Convert this polynomial to a function returning nullable value (null if argument is outside piecewise range).
*/ */
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { substitute(ring, it) } public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { value(ring, it) }
/** /**
* Convert this polynomial to a function using [defaultValue] for arguments outside the piecewise range. * Convert this polynomial to a function using [defaultValue] for arguments outside the piecewise range.
*/ */
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C, defaultValue: T): (T) -> T = public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C, defaultValue: T): (T) -> T =
{ substitute(ring, it) ?: defaultValue } { value(ring, it) ?: defaultValue }

View File

@ -59,12 +59,12 @@ public value class Polynomial<C>(
* @param A type of provided underlying ring of constants. It's [Ring] of [C]. * @param A type of provided underlying ring of constants. It's [Ring] of [C].
* @param ring underlying ring of constants of type [A]. * @param ring underlying ring of constants of type [A].
*/ */
public open class PolynomialSpace<C, A : Ring<C>>( public open class PolynomialSpace<C, A>(
/** /**
* Underlying ring of constants. Its operations on constants are inherited by local operations on constants. * Underlying ring of constants. Its operations on constants are inherited by local operations on constants.
*/ */
public val ring: A, public val ring: A,
) : Ring<Polynomial<C>> { ) : Ring<Polynomial<C>>, ScaleOperations<Polynomial<C>> where A : Ring<C>, A : ScaleOperations<C> {
/** /**
* Instance of zero constant (zero of the underlying ring). * Instance of zero constant (zero of the underlying ring).
@ -267,13 +267,15 @@ public open class PolynomialSpace<C, A : Ring<C>>(
override fun add(left: Polynomial<C>, right: Polynomial<C>): Polynomial<C> = left + right override fun add(left: Polynomial<C>, right: Polynomial<C>): Polynomial<C> = left + right
override fun multiply(left: Polynomial<C>, right: Polynomial<C>): Polynomial<C> = left * right override fun multiply(left: Polynomial<C>, right: Polynomial<C>): Polynomial<C> = left * right
override fun scale(a: Polynomial<C>, value: Double): Polynomial<C> =
ring { Polynomial(a.coefficients.map { scale(it, value) }) }
// TODO: When context receivers will be ready move all of this substitutions and invocations to utilities with // TODO: When context receivers will be ready move all of this substitutions and invocations to utilities with
// [ListPolynomialSpace] as a context receiver // [ListPolynomialSpace] as a context receiver
/** /**
* Evaluates value of [this] polynomial on provided [argument]. * Evaluates value of [this] polynomial on provided [argument].
*/ */
public inline fun Polynomial<C>.substitute(argument: C): C = substitute(ring, argument) public inline fun Polynomial<C>.substitute(argument: C): C = value(ring, argument)
/** /**
* Represent [this] polynomial as a regular context-less function. * Represent [this] polynomial as a regular context-less function.
@ -283,19 +285,5 @@ public open class PolynomialSpace<C, A : Ring<C>>(
/** /**
* Evaluates value of [this] polynomial on provided [argument]. * Evaluates value of [this] polynomial on provided [argument].
*/ */
public inline operator fun Polynomial<C>.invoke(argument: C): C = substitute(ring, argument) public inline operator fun Polynomial<C>.invoke(argument: C): C = value(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 ScalablePolynomialSpace<C, A>(
ring: A,
) : PolynomialSpace<C, A>(ring), ScaleOperations<Polynomial<C>> where A : Ring<C>, A : ScaleOperations<C> {
override fun scale(a: Polynomial<C>, value: Double): Polynomial<C> =
ring { Polynomial(a.coefficients.map { scale(it, value) }) }
} }

View File

@ -16,36 +16,22 @@ import kotlin.math.pow
/** /**
* Creates a [PolynomialSpace] over a received ring. * Creates a [PolynomialSpace] over a received ring.
*/ */
public inline val <C, A : Ring<C>> A.polynomialSpace: PolynomialSpace<C, A> public inline val <C, A> A.polynomialSpace: PolynomialSpace<C, A> where A : Ring<C>, A : ScaleOperations<C>
get() = PolynomialSpace(this) get() = PolynomialSpace(this)
/** /**
* Creates a [PolynomialSpace]'s scope over a received ring. * Creates a [PolynomialSpace]'s scope over a received ring.
*/ // TODO: When context will be ready move [ListPolynomialSpace] and add [A] to context receivers of [block] */ // TODO: When context will be ready move [ListPolynomialSpace] and add [A] to context receivers of [block]
public inline fun <C, A : Ring<C>, R> A.polynomialSpace(block: PolynomialSpace<C, A>.() -> R): R { public inline fun <C, A, R> A.polynomialSpace(block: PolynomialSpace<C, A>.() -> R): R where A : Ring<C>, A : ScaleOperations<C> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) } contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return PolynomialSpace(this).block() return PolynomialSpace(this).block()
} }
/**
* Creates a [ScalablePolynomialSpace] over a received scalable ring.
*/
public inline val <C, A> A.scalablePolynomialSpace: ScalablePolynomialSpace<C, A> where A : Ring<C>, A : ScaleOperations<C>
get() = ScalablePolynomialSpace(this)
/**
* Creates a [ScalablePolynomialSpace]'s scope over a received scalable ring.
*/ // TODO: When context will be ready move [ListPolynomialSpace] and add [A] to context receivers of [block]
public inline fun <C, A, R> A.scalablePolynomialSpace(block: ScalablePolynomialSpace<C, A>.() -> R): R where A : Ring<C>, A : ScaleOperations<C> {
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
return ScalablePolynomialSpace(this).block()
}
/** /**
* Evaluates value of [this] Double polynomial on provided Double argument. * Evaluates value of [this] Double polynomial on provided Double argument.
*/ */
public fun Polynomial<Double>.substitute(arg: Double): Double = public fun Polynomial<Double>.value(arg: Double): Double =
coefficients.reduceIndexedOrNull { index, acc, c -> coefficients.reduceIndexedOrNull { index, acc, c ->
acc + c * arg.pow(index) acc + c * arg.pow(index)
} ?: .0 } ?: .0
@ -55,7 +41,7 @@ public fun Polynomial<Double>.substitute(arg: Double): Double =
* *
* It is an implementation of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method). * It is an implementation of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
*/ */
public fun <C> Polynomial<C>.substitute(ring: Ring<C>, arg: C): C = ring { public fun <C> Polynomial<C>.value(ring: Ring<C>, arg: C): C = ring {
if (coefficients.isEmpty()) return zero if (coefficients.isEmpty()) return zero
var result: C = coefficients.last() var result: C = coefficients.last()
for (j in coefficients.size - 2 downTo 0) { for (j in coefficients.size - 2 downTo 0) {
@ -67,13 +53,13 @@ public fun <C> Polynomial<C>.substitute(ring: Ring<C>, arg: C): C = ring {
/** /**
* Represent [this] polynomial as a regular context-less function. * Represent [this] polynomial as a regular context-less function.
*/ */
public fun <C, A : Ring<C>> Polynomial<C>.asFunctionOver(ring: A): (C) -> C = { substitute(ring, it) } public fun <C, A : Ring<C>> Polynomial<C>.asFunctionOver(ring: A): (C) -> C = { value(ring, it) }
/** /**
* Returns algebraic derivative of received polynomial. * Returns algebraic derivative of received polynomial.
*/ */
@UnstableKMathAPI @UnstableKMathAPI
public fun <C, A> Polynomial<C>.derivative( public fun <C, A> Polynomial<C>.differentiate(
ring: A, ring: A,
): Polynomial<C> where A : Ring<C>, A : NumericAlgebra<C> = ring { ): Polynomial<C> where A : Ring<C>, A : NumericAlgebra<C> = ring {
Polynomial( Polynomial(
@ -87,7 +73,7 @@ public fun <C, A> Polynomial<C>.derivative(
* Returns algebraic antiderivative of received polynomial. * Returns algebraic antiderivative of received polynomial.
*/ */
@UnstableKMathAPI @UnstableKMathAPI
public fun <C, A> Polynomial<C>.antiderivative( public fun <C, A> Polynomial<C>.integrate(
ring: A, ring: A,
): Polynomial<C> where A : Field<C>, A : NumericAlgebra<C> = ring { ): Polynomial<C> where A : Field<C>, A : NumericAlgebra<C> = ring {
Polynomial( Polynomial(
@ -106,6 +92,6 @@ public fun <C : Comparable<C>> Polynomial<C>.integrate(
ring: Field<C>, ring: Field<C>,
range: ClosedRange<C>, range: ClosedRange<C>,
): C { ): C {
val antiderivative = antiderivative(ring) val antiderivative = integrate(ring)
return ring { antiderivative.substitute(ring, range.endInclusive) - antiderivative.substitute(ring, range.start) } return ring { antiderivative.value(ring, range.endInclusive) - antiderivative.value(ring, range.start) }
} }

View File

@ -7,7 +7,6 @@ package space.kscience.kmath.integration
import space.kscience.kmath.functions.PiecewisePolynomial import space.kscience.kmath.functions.PiecewisePolynomial
import space.kscience.kmath.functions.integrate import space.kscience.kmath.functions.integrate
import space.kscience.kmath.functions.antiderivative
import space.kscience.kmath.interpolation.PolynomialInterpolator import space.kscience.kmath.interpolation.PolynomialInterpolator
import space.kscience.kmath.interpolation.SplineInterpolator import space.kscience.kmath.interpolation.SplineInterpolator
import space.kscience.kmath.interpolation.interpolatePolynomials import space.kscience.kmath.interpolation.interpolatePolynomials
@ -24,7 +23,7 @@ import space.kscience.kmath.structures.MutableBufferFactory
@OptIn(PerformancePitfall::class) @OptIn(PerformancePitfall::class)
@UnstableKMathAPI @UnstableKMathAPI
public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(algebra: Field<T>): PiecewisePolynomial<T> = public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(algebra: Field<T>): PiecewisePolynomial<T> =
PiecewisePolynomial(pieces.map { it.first to it.second.antiderivative(algebra) }) PiecewisePolynomial(pieces.map { it.first to it.second.integrate(algebra) })
/** /**
* Compute definite integral of given [PiecewisePolynomial] piece by piece in a given [range] * Compute definite integral of given [PiecewisePolynomial] piece by piece in a given [range]

View File

@ -10,7 +10,7 @@ package space.kscience.kmath.interpolation
import space.kscience.kmath.data.XYColumnarData import space.kscience.kmath.data.XYColumnarData
import space.kscience.kmath.functions.PiecewisePolynomial import space.kscience.kmath.functions.PiecewisePolynomial
import space.kscience.kmath.functions.asFunction import space.kscience.kmath.functions.asFunction
import space.kscience.kmath.functions.substitute import space.kscience.kmath.functions.value
import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.Ring import space.kscience.kmath.operations.Ring
import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.Buffer
@ -34,7 +34,7 @@ public interface PolynomialInterpolator<T : Comparable<T>> : Interpolator<T, T,
public fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T> public fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T>
override fun interpolate(points: XYColumnarData<T, T, T>): (T) -> T = { x -> override fun interpolate(points: XYColumnarData<T, T, T>): (T) -> T = { x ->
interpolatePolynomials(points).substitute(algebra, x) ?: getDefaultValue() interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue()
} }
} }

View File

@ -15,119 +15,119 @@ import kotlin.test.assertEquals
@OptIn(UnstableKMathAPI::class) @OptIn(UnstableKMathAPI::class)
class PolynomialUtilTest { class PolynomialUtilTest {
@Test @Test
fun test_Polynomial_substitute_Double() { fun test_Polynomial_value_Double() {
assertEquals( assertEquals(
0.0, 0.0,
Polynomial(1.0, -2.0, 1.0).substitute(1.0), Polynomial(1.0, -2.0, 1.0).value(1.0),
0.001, 0.001,
"test 1" "test 1"
) )
assertEquals( assertEquals(
0.0, 0.0,
Polynomial(1.0, -2.0, 1.0).substitute(1.0), Polynomial(1.0, -2.0, 1.0).value(1.0),
0.001, 0.001,
"test 1" "test 1"
) )
assertEquals( assertEquals(
1.1931904761904761, 1.1931904761904761,
Polynomial(0.625, 2.6666666666666665, 0.5714285714285714, 1.5).substitute(0.2), Polynomial(0.625, 2.6666666666666665, 0.5714285714285714, 1.5).value(0.2),
0.001, 0.001,
"test 2" "test 2"
) )
assertEquals( assertEquals(
0.5681904761904762, 0.5681904761904762,
Polynomial(0.0, 2.6666666666666665, 0.5714285714285714, 1.5).substitute(0.2), Polynomial(0.0, 2.6666666666666665, 0.5714285714285714, 1.5).value(0.2),
0.001, 0.001,
"test 3" "test 3"
) )
assertEquals( assertEquals(
1.1811904761904761, 1.1811904761904761,
Polynomial(0.625, 2.6666666666666665, 0.5714285714285714, 0.0).substitute(0.2), Polynomial(0.625, 2.6666666666666665, 0.5714285714285714, 0.0).value(0.2),
0.001, 0.001,
"test 4" "test 4"
) )
assertEquals( assertEquals(
1.1703333333333332, 1.1703333333333332,
Polynomial(0.625, 2.6666666666666665, 0.0, 1.5).substitute(0.2), Polynomial(0.625, 2.6666666666666665, 0.0, 1.5).value(0.2),
0.001, 0.001,
"test 5" "test 5"
) )
} }
@Test @Test
fun test_Polynomial_substitute_Constant() { fun test_Polynomial_value_Constant() {
assertEquals( assertEquals(
Rational(0), Rational(0),
Polynomial(Rational(1), Rational(-2), Rational(1)).substitute(RationalField, Rational(1)), Polynomial(Rational(1), Rational(-2), Rational(1)).value(RationalField, Rational(1)),
"test 1" "test 1"
) )
assertEquals( assertEquals(
Rational(25057, 21000), Rational(25057, 21000),
Polynomial(Rational(5, 8), Rational(8, 3), Rational(4, 7), Rational(3, 2)) Polynomial(Rational(5, 8), Rational(8, 3), Rational(4, 7), Rational(3, 2))
.substitute(RationalField, Rational(1, 5)), .value(RationalField, Rational(1, 5)),
"test 2" "test 2"
) )
assertEquals( assertEquals(
Rational(2983, 5250), Rational(2983, 5250),
Polynomial(Rational(0), Rational(8, 3), Rational(4, 7), Rational(3, 2)) Polynomial(Rational(0), Rational(8, 3), Rational(4, 7), Rational(3, 2))
.substitute(RationalField, Rational(1, 5)), .value(RationalField, Rational(1, 5)),
"test 3" "test 3"
) )
assertEquals( assertEquals(
Rational(4961, 4200), Rational(4961, 4200),
Polynomial(Rational(5, 8), Rational(8, 3), Rational(4, 7), Rational(0)) Polynomial(Rational(5, 8), Rational(8, 3), Rational(4, 7), Rational(0))
.substitute(RationalField, Rational(1, 5)), .value(RationalField, Rational(1, 5)),
"test 4" "test 4"
) )
assertEquals( assertEquals(
Rational(3511, 3000), Rational(3511, 3000),
Polynomial(Rational(5, 8), Rational(8, 3), Rational(0), Rational(3, 2)) Polynomial(Rational(5, 8), Rational(8, 3), Rational(0), Rational(3, 2))
.substitute(RationalField, Rational(1, 5)), .value(RationalField, Rational(1, 5)),
"test 5" "test 5"
) )
} }
@Test @Test
fun test_Polynomial_derivative() { fun test_Polynomial_differentiate() {
assertEquals( assertEquals(
Polynomial(Rational(-2), Rational(2)), Polynomial(Rational(-2), Rational(2)),
Polynomial(Rational(1), Rational(-2), Rational(1)).derivative(RationalField), Polynomial(Rational(1), Rational(-2), Rational(1)).differentiate(RationalField),
"test 1" "test 1"
) )
assertEquals( assertEquals(
Polynomial(Rational(-8, 3), Rational(8, 9), Rational(15, 7), Rational(-20, 9)), Polynomial(Rational(-8, 3), Rational(8, 9), Rational(15, 7), Rational(-20, 9)),
Polynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).derivative(RationalField), Polynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).differentiate(RationalField),
"test 2" "test 2"
) )
assertEquals( assertEquals(
Polynomial(Rational(0), Rational(8, 9), Rational(15, 7), Rational(-20, 9)), Polynomial(Rational(0), Rational(8, 9), Rational(15, 7), Rational(-20, 9)),
Polynomial(Rational(0), Rational(0), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).derivative(RationalField), Polynomial(Rational(0), Rational(0), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).differentiate(RationalField),
"test 3" "test 3"
) )
assertEquals( assertEquals(
Polynomial(Rational(-8, 3), Rational(8, 9), Rational(15, 7), Rational(0)), Polynomial(Rational(-8, 3), Rational(8, 9), Rational(15, 7), Rational(0)),
Polynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(0)).derivative(RationalField), Polynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(0)).differentiate(RationalField),
"test 4" "test 4"
) )
} }
@Test @Test
fun test_Polynomial_antiderivative() { fun test_Polynomial_integrate() {
assertEquals( assertEquals(
Polynomial(Rational(0), Rational(1), Rational(-1), Rational(1, 3)), Polynomial(Rational(0), Rational(1), Rational(-1), Rational(1, 3)),
Polynomial(Rational(1), Rational(-2), Rational(1)).antiderivative(RationalField), Polynomial(Rational(1), Rational(-2), Rational(1)).integrate(RationalField),
"test 1" "test 1"
) )
assertEquals( assertEquals(
Polynomial(Rational(0), Rational(1, 5), Rational(-4, 3), Rational(4, 27), Rational(5, 28), Rational(-1, 9)), Polynomial(Rational(0), Rational(1, 5), Rational(-4, 3), Rational(4, 27), Rational(5, 28), Rational(-1, 9)),
Polynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).antiderivative(RationalField), Polynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).integrate(RationalField),
"test 2" "test 2"
) )
assertEquals( assertEquals(
Polynomial(Rational(0), Rational(0), Rational(0), Rational(4, 27), Rational(5, 28), Rational(-1, 9)), Polynomial(Rational(0), Rational(0), Rational(0), Rational(4, 27), Rational(5, 28), Rational(-1, 9)),
Polynomial(Rational(0), Rational(0), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).antiderivative(RationalField), Polynomial(Rational(0), Rational(0), Rational(4, 9), Rational(5, 7), Rational(-5, 9)).integrate(RationalField),
"test 3" "test 3"
) )
assertEquals( assertEquals(
Polynomial(Rational(0), Rational(1, 5), Rational(-4, 3), Rational(4, 27), Rational(5, 28), Rational(0)), Polynomial(Rational(0), Rational(1, 5), Rational(-4, 3), Rational(4, 27), Rational(5, 28), Rational(0)),
Polynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(0)).antiderivative(RationalField), Polynomial(Rational(1, 5), Rational(-8, 3), Rational(4, 9), Rational(5, 7), Rational(0)).integrate(RationalField),
"test 4" "test 4"
) )
} }