From ed2f14b68e7f799fcbe370bf744256a278cfcf74 Mon Sep 17 00:00:00 2001 From: Gleb Minaev <43728100+lounres@users.noreply.github.com> Date: Fri, 18 Mar 2022 01:47:03 +0300 Subject: [PATCH] Optimised existing substitution function. Prototyped substitution for RFs. --- .../kmath/functions/polynomialUtil.kt | 63 ++++--- .../kmath/functions/rationalFunctionUtil.kt | 174 ++++++++++++++++++ 2 files changed, 210 insertions(+), 27 deletions(-) diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt index c624380be..1a7245b9f 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt @@ -50,24 +50,36 @@ public inline fun A.scalablePolynomial(block: ScalablePolynomialSpace< } @Suppress("NOTHING_TO_INLINE") -internal inline fun iadd( - ring: Ring, - augend: MutableList, - addend: List, - degree: Int -) = ring { - for (deg in 0 .. degree) augend[deg] += addend[deg] +internal inline fun copyTo( + origin: List, + originDegree: Int, + target: MutableList, +) { + for (deg in 0 .. originDegree) target[deg] = origin[deg] } @Suppress("NOTHING_TO_INLINE") -internal inline fun addTo( +internal inline fun multiplyAddingToUpdater( ring: Ring, - augend: List, - addend: List, - degree: Int, - target: MutableList -) = ring { - for (deg in 0 .. degree) target[deg] = augend[deg] + addend[deg] + 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 + } } @Suppress("NOTHING_TO_INLINE") @@ -107,32 +119,29 @@ public fun Polynomial.substitute(ring: Ring, arg: C): C = ring { return result } -public fun Polynomial.substitute(ring: Ring, arg: Polynomial) : Polynomial = ring.polynomial { - if (coefficients.isEmpty()) return zero +public fun Polynomial.substitute(ring: Ring, arg: Polynomial) : Polynomial = ring { + if (coefficients.isEmpty()) return Polynomial(emptyList()) - val thisDegree = degree - if (thisDegree == -1) return zero - val argDegree = arg.degree + val thisDegree = coefficients.indexOfLast { it != zero } + if (thisDegree == -1) return Polynomial(emptyList()) + val argDegree = arg.coefficients.indexOfLast { it != zero } if (argDegree == -1) return coefficients[0].asPolynomial() - val constantZero = constantZero + val constantZero = zero val resultCoefs: MutableList = MutableList(thisDegree * argDegree + 1) { constantZero } val resultCoefsUpdate: MutableList = MutableList(thisDegree * argDegree + 1) { constantZero } var resultDegree = 0 for (deg in thisDegree downTo 0) { resultCoefsUpdate[0] = coefficients[deg] - multiplyAddingTo( - ring=ring, + multiplyAddingToUpdater( + ring = ring, multiplicand = resultCoefs, multiplicandDegree = resultDegree, multiplier = arg.coefficients, multiplierDegree = argDegree, - target = resultCoefsUpdate + updater = resultCoefsUpdate, + zero = constantZero ) resultDegree += argDegree - for (updateDeg in 0 .. resultDegree) { - resultCoefs[updateDeg] = resultCoefsUpdate[updateDeg] - resultCoefsUpdate[updateDeg] = constantZero - } } return Polynomial(resultCoefs) } diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/rationalFunctionUtil.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/rationalFunctionUtil.kt index 456a1fa2b..8d2073834 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/rationalFunctionUtil.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/rationalFunctionUtil.kt @@ -5,9 +5,12 @@ package space.kscience.kmath.functions +import space.kscience.kmath.operations.Field import space.kscience.kmath.operations.Ring +import space.kscience.kmath.operations.invoke import kotlin.contracts.InvocationKind import kotlin.contracts.contract +import kotlin.math.max /** @@ -24,6 +27,177 @@ public inline fun , R> A.rationalFunction(block: RationalFunction return RationalFunctionSpace(this).block() } +/** + * Evaluates the value of the given double polynomial for given double argument. + */ +public fun RationalFunction.substitute(arg: Double): Double = + numerator.substitute(arg) / denominator.substitute(arg) + +/** + * Evaluates the value of the given polynomial for given argument. + * + * It is an implementation of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method). + */ +public fun RationalFunction.substitute(ring: Field, arg: C): C = ring { + numerator.substitute(ring, arg) / denominator.substitute(ring, arg) +} + +/** + * 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 [Polynomial.substitute] and [RationalFunction.substitute] for performance optimisation. + */ // TODO: Дописать +internal fun Polynomial.substituteRationalFunctionTakeNumerator(ring: Ring, arg: RationalFunction): Polynomial = ring { + if (coefficients.isEmpty()) return Polynomial(emptyList()) + + val thisDegree = coefficients.indexOfLast { it != zero } + if (thisDegree == -1) return Polynomial(emptyList()) + val thisDegreeLog2 = 31 - thisDegree.countLeadingZeroBits() + val numeratorDegree = arg.numerator.coefficients.indexOfLast { it != zero } + val denominatorDegree = arg.denominator.coefficients.indexOfLast { it != zero } + 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 Polynomial( + processLevelEdged( + level = thisDegreeLog2, + start = 0, + end = thisDegree + 1 + ) + ) +} + //operator fun > RationalFunction.invoke(arg: T): T = numerator(arg) / denominator(arg) // //fun > RationalFunction.reduced(): RationalFunction =