kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialUtil.kt

@@ -50,24 +50,36 @@ public inline fun <C, A, R> A.scalablePolynomial(block: ScalablePolynomialSpace<
 }
 
 @Suppress("NOTHING_TO_INLINE")
-internal inline fun <C> iadd(
-    ring: Ring<C>,
-    augend: MutableList<C>,
-    addend: List<C>,
-    degree: Int
-) = ring {
-    for (deg in 0 .. degree) augend[deg] += addend[deg]
+internal inline fun <C> copyTo(
+    origin: List<C>,
+    originDegree: Int,
+    target: MutableList<C>,
+) {
+    for (deg in 0 .. originDegree) target[deg] = origin[deg]
 }
 
 @Suppress("NOTHING_TO_INLINE")
-internal inline fun <C> addTo(
+internal inline fun <C> multiplyAddingToUpdater(
     ring: Ring<C>,
-    augend: List<C>,
-    addend: List<C>,
-    degree: Int,
-    target: MutableList<C>
-) = ring {
-    for (deg in 0 .. degree) target[deg] = augend[deg] + addend[deg]
+    multiplicand: MutableList<C>,
+    multiplicandDegree: Int,
+    multiplier: List<C>,
+    multiplierDegree: Int,
+    updater: MutableList<C>,
+    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 <C> Polynomial<C>.substitute(ring: Ring<C>, arg: C): C = ring {
     return result
 }
 
-public fun <C> Polynomial<C>.substitute(ring: Ring<C>, arg: Polynomial<C>) : Polynomial<C> = ring.polynomial {
-    if (coefficients.isEmpty()) return zero
+public fun <C> Polynomial<C>.substitute(ring: Ring<C>, arg: Polynomial<C>) : Polynomial<C> = 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<C> = MutableList(thisDegree * argDegree + 1) { constantZero }
     val resultCoefsUpdate: MutableList<C> = MutableList(thisDegree * argDegree + 1) { constantZero }
     var resultDegree = 0
     for (deg in thisDegree downTo 0) {
         resultCoefsUpdate[0] = coefficients[deg]
-        multiplyAddingTo(
+        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<C>(resultCoefs)
 }

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 <C, A : Ring<C>, R> A.rationalFunction(block: RationalFunction
     return RationalFunctionSpace(this).block()
 }
 
+/**
+ * Evaluates the value of the given double polynomial for given double argument.
+ */
+public fun RationalFunction<Double>.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 <C> RationalFunction<C>.substitute(ring: Field<C>, 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 <C> Polynomial<C>.substituteRationalFunctionTakeNumerator(ring: Ring<C>, arg: RationalFunction<C>): Polynomial<C> = 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<Int>(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<List<C>>(thisDegreeLog2 + 1) {
+        add(arg.numerator.coefficients)
+        repeat(thisDegreeLog2) {
+            val next = MutableList<C>(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<List<C>>(thisDegreeLog2 + 1) {
+        add(arg.denominator.coefficients)
+        repeat(thisDegreeLog2) {
+            val next = MutableList<C>(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<MutableList<C>>(thisDegreeLog2 + 1) {
+        repeat(thisDegreeLog2 + 1) {
+            add(MutableList(hashes[it] * argDegree) { constantZero })
+        }
+    }
+    val edgedMultiplier = MutableList<C>(0) { TODO() }
+    val edgedMultiplierUpdater = MutableList<C>(0) { TODO() }
+
+    fun MutableList<C>.reset() {
+        for (i in indices) set(i, constantZero)
+    }
+
+    fun processLevel(level: Int, start: Int, end: Int) : List<C> {
+        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<C> {
+        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 <T: Field<T>> RationalFunction<T>.invoke(arg: T): T = numerator(arg) / denominator(arg)
 //
 //fun <T: Field<T>> RationalFunction<T>.reduced(): RationalFunction<T> =