diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt
new file mode 100644
index 000000000..612b00535
--- /dev/null
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt
@@ -0,0 +1,132 @@
+/*
+ * 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.functions
+
+import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.operations.Ring
+
+/**
+ * Represents piecewise-defined function.
+ *
+ * @param T the piece key type.
+ * @param R the sub-function type.
+ */
+public fun interface Piecewise<in T, out R> {
+    /**
+     * Returns the appropriate sub-function for given piece key.
+     */
+    public fun findPiece(arg: T): R?
+}
+
+/**
+ * Represents piecewise-defined function where all the sub-functions are polynomials.
+ *
+ * @property pieces An ordered list of range-polynomial pairs. The list does not in general guarantee that there are no
+ * "holes" in it.
+ */
+public interface PiecewisePolynomial<T : Comparable<T>> : Piecewise<T, ListPolynomial<T>> {
+    public val pieces: Collection<Pair<ClosedRange<T>, ListPolynomial<T>>>
+
+    override fun findPiece(arg: T): ListPolynomial<T>?
+}
+
+/**
+ * A generic piecewise without constraints on how pieces are placed
+ */
+@PerformancePitfall("findPiece method of resulting piecewise is slow")
+public fun <T : Comparable<T>> PiecewisePolynomial(
+    pieces: Collection<Pair<ClosedRange<T>, ListPolynomial<T>>>,
+): PiecewisePolynomial<T> = object : PiecewisePolynomial<T> {
+    override val pieces: Collection<Pair<ClosedRange<T>, ListPolynomial<T>>> = pieces
+
+    override fun findPiece(arg: T): ListPolynomial<T>? = pieces.firstOrNull { arg in it.first }?.second
+}
+
+/**
+ * An optimized piecewise that uses not separate pieces, but a range separated by delimiters.
+ * The pieces search is logarithmic.
+ */
+private class OrderedPiecewisePolynomial<T : Comparable<T>>(
+    override val pieces: List<Pair<ClosedRange<T>, ListPolynomial<T>>>,
+) : PiecewisePolynomial<T> {
+
+    override fun findPiece(arg: T): ListPolynomial<T>? {
+        val index = pieces.binarySearch { (range, _) ->
+            when {
+                arg >= range.endInclusive -> -1
+                arg < range.start -> +1
+                else -> 0
+            }
+        }
+        return if (index < 0) null else pieces[index].second
+    }
+
+}
+
+/**
+ * A [Piecewise]  builder where all the pieces are ordered by the [Comparable] type instances.
+ *
+ * @param T the comparable piece key type.
+ * @param delimiter the initial piecewise separator
+ */
+public class PiecewiseBuilder<T : Comparable<T>>(delimiter: T) {
+    private val delimiters: MutableList<T> = arrayListOf(delimiter)
+    private val pieces: MutableList<ListPolynomial<T>> = arrayListOf()
+
+    /**
+     * Dynamically adds a piece to the right side (beyond maximum argument value of previous piece)
+     *
+     * @param right new rightmost position. If is less than current rightmost position, an error is thrown.
+     * @param piece the sub-function.
+     */
+    public fun putRight(right: T, piece: ListPolynomial<T>) {
+        require(right > delimiters.last()) { "New delimiter should be to the right of old one" }
+        delimiters += right
+        pieces += piece
+    }
+
+    /**
+     * Dynamically adds a piece to the left side (beyond maximum argument value of previous piece)
+     *
+     * @param left the new leftmost position. If is less than current rightmost position, an error is thrown.
+     * @param piece the sub-function.
+     */
+    public fun putLeft(left: T, piece: ListPolynomial<T>) {
+        require(left < delimiters.first()) { "New delimiter should be to the left of old one" }
+        delimiters.add(0, left)
+        pieces.add(0, piece)
+    }
+
+    public fun build(): PiecewisePolynomial<T> = OrderedPiecewisePolynomial(delimiters.zipWithNext { l, r ->
+        l..r
+    }.zip(pieces))
+}
+
+/**
+ * A builder for [PiecewisePolynomial]
+ */
+public fun <T : Comparable<T>> PiecewisePolynomial(
+    startingPoint: T,
+    builder: PiecewiseBuilder<T>.() -> Unit,
+): PiecewisePolynomial<T> = PiecewiseBuilder(startingPoint).apply(builder).build()
+
+/**
+ * Return a value of polynomial function with given [ring] a given [arg] or null if argument is outside piecewise
+ * definition.
+ */
+public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.substitute(ring: C, arg: T): T? =
+    findPiece(arg)?.substitute(ring, arg)
+
+/**
+ * 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) }
+
+/**
+ * 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 =
+    { substitute(ring, it) ?: defaultValue }
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt
new file mode 100644
index 000000000..80006c2de
--- /dev/null
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt
@@ -0,0 +1,112 @@
+/*
+ * 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.
+ */
+
+/*
+ * 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.integration
+
+import space.kscience.kmath.functions.PiecewisePolynomial
+import space.kscience.kmath.functions.integrate
+import space.kscience.kmath.functions.antiderivative
+import space.kscience.kmath.interpolation.PolynomialInterpolator
+import space.kscience.kmath.interpolation.SplineInterpolator
+import space.kscience.kmath.interpolation.interpolatePolynomials
+import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.MutableBufferFactory
+
+/**
+ * Compute analytical indefinite integral of this [PiecewisePolynomial], keeping all intervals intact
+ */
+@OptIn(PerformancePitfall::class)
+@UnstableKMathAPI
+public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(algebra: Field<T>): PiecewisePolynomial<T> =
+    PiecewisePolynomial(pieces.map { it.first to it.second.antiderivative(algebra) })
+
+/**
+ * Compute definite integral of given [PiecewisePolynomial] piece by piece in a given [range]
+ * Requires [UnivariateIntegrationNodes] or [IntegrationRange] and [IntegrandMaxCalls]
+ *
+ * TODO use context receiver for algebra
+ */
+@UnstableKMathAPI
+public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(
+    algebra: Field<T>, range: ClosedRange<T>,
+): T = algebra.sum(
+    pieces.map { (region, poly) ->
+        val intersectedRange = maxOf(range.start, region.start)..minOf(range.endInclusive, region.endInclusive)
+        //Check if polynomial range is not used
+        if (intersectedRange.start == intersectedRange.endInclusive) algebra.zero
+        else poly.integrate(algebra, intersectedRange)
+    }
+)
+
+/**
+ * A generic spline-interpolation-based analytic integration
+ * * [IntegrationRange]&mdash;the univariate range of integration. By default, uses `0..1` interval.
+ * * [IntegrandMaxCalls]&mdash;the maximum number of function calls during integration. For non-iterative rules, always uses
+ * the maximum number of points. By default, uses 10 points.
+ */
+@UnstableKMathAPI
+public class SplineIntegrator<T : Comparable<T>>(
+    public val algebra: Field<T>,
+    public val bufferFactory: MutableBufferFactory<T>,
+) : UnivariateIntegrator<T> {
+    override fun process(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = algebra {
+        val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
+
+        val interpolator: PolynomialInterpolator<T> = SplineInterpolator(algebra, bufferFactory)
+
+        val nodes: Buffer<Double> = integrand.getFeature<UnivariateIntegrationNodes>()?.nodes ?: run {
+            val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
+            val step = (range.endInclusive - range.start) / (numPoints - 1)
+            DoubleBuffer(numPoints) { i -> range.start + i * step }
+        }
+
+        val values = nodes.map(bufferFactory) { integrand.function(it) }
+        val polynomials = interpolator.interpolatePolynomials(
+            nodes.map(bufferFactory) { number(it) },
+            values
+        )
+        val res = polynomials.integrate(algebra, number(range.start)..number(range.endInclusive))
+        integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + nodes.size)
+    }
+}
+
+/**
+ * A simplified double-based spline-interpolation-based analytic integration
+ * * [IntegrationRange]&mdash;the univariate range of integration. By default, uses `0.0..1.0` interval.
+ * * [IntegrandMaxCalls]&mdash;the maximum number of function calls during integration. For non-iterative rules, always
+ * uses the maximum number of points. By default, uses 10 points.
+ */
+@UnstableKMathAPI
+public object DoubleSplineIntegrator : UnivariateIntegrator<Double> {
+    override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
+        val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
+        val interpolator: PolynomialInterpolator<Double> = SplineInterpolator(DoubleField, ::DoubleBuffer)
+
+        val nodes: Buffer<Double> = integrand.getFeature<UnivariateIntegrationNodes>()?.nodes ?: run {
+            val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
+            val step = (range.endInclusive - range.start) / (numPoints - 1)
+            DoubleBuffer(numPoints) { i -> range.start + i * step }
+        }
+
+        val values = nodes.map { integrand.function(it) }
+        val polynomials = interpolator.interpolatePolynomials(nodes, values)
+        val res = polynomials.integrate(DoubleField, range)
+        return integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + nodes.size)
+    }
+}
+
+@Suppress("unused")
+@UnstableKMathAPI
+public inline val DoubleField.splineIntegrator: UnivariateIntegrator<Double>
+    get() = DoubleSplineIntegrator
\ No newline at end of file
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt
new file mode 100644
index 000000000..62819be0c
--- /dev/null
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt
@@ -0,0 +1,92 @@
+/*
+ * 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.
+ */
+
+@file:OptIn(UnstableKMathAPI::class)
+
+package space.kscience.kmath.interpolation
+
+import space.kscience.kmath.data.XYColumnarData
+import space.kscience.kmath.functions.PiecewisePolynomial
+import space.kscience.kmath.functions.asFunction
+import space.kscience.kmath.functions.substitute
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.Ring
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.asBuffer
+
+/**
+ * And interpolator for data with x column type [X], y column type [Y].
+ */
+public fun interface Interpolator<T, in X : T, Y : T> {
+    public fun interpolate(points: XYColumnarData<T, X, Y>): (X) -> Y
+}
+
+/**
+ * And interpolator returning [PiecewisePolynomial] function
+ */
+public interface PolynomialInterpolator<T : Comparable<T>> : Interpolator<T, T, T> {
+    public val algebra: Ring<T>
+
+    public fun getDefaultValue(): T = error("Out of bounds")
+
+    public fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T>
+
+    override fun interpolate(points: XYColumnarData<T, T, T>): (T) -> T = { x ->
+        interpolatePolynomials(points).substitute(algebra, x) ?: getDefaultValue()
+    }
+}
+
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
+    x: Buffer<T>,
+    y: Buffer<T>,
+): PiecewisePolynomial<T> {
+    val pointSet = XYColumnarData.of(x, y)
+    return interpolatePolynomials(pointSet)
+}
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
+    data: Map<T, T>,
+): PiecewisePolynomial<T> {
+    val pointSet = XYColumnarData.of(data.keys.toList().asBuffer(), data.values.toList().asBuffer())
+    return interpolatePolynomials(pointSet)
+}
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
+    data: List<Pair<T, T>>,
+): PiecewisePolynomial<T> {
+    val pointSet = XYColumnarData.of(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer())
+    return interpolatePolynomials(pointSet)
+}
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    x: Buffer<T>,
+    y: Buffer<T>,
+): (T) -> T? = interpolatePolynomials(x, y).asFunction(algebra)
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    data: Map<T, T>,
+): (T) -> T? = interpolatePolynomials(data).asFunction(algebra)
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    data: List<Pair<T, T>>,
+): (T) -> T? = interpolatePolynomials(data).asFunction(algebra)
+
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    x: Buffer<T>,
+    y: Buffer<T>,
+    defaultValue: T,
+): (T) -> T = interpolatePolynomials(x, y).asFunction(algebra, defaultValue)
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    data: Map<T, T>,
+    defaultValue: T,
+): (T) -> T = interpolatePolynomials(data).asFunction(algebra, defaultValue)
+
+public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
+    data: List<Pair<T, T>>,
+    defaultValue: T,
+): (T) -> T = interpolatePolynomials(data).asFunction(algebra, defaultValue)
\ No newline at end of file
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt
new file mode 100644
index 000000000..b55f16cf2
--- /dev/null
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt
@@ -0,0 +1,43 @@
+/*
+ * 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.interpolation
+
+import space.kscience.kmath.data.XYColumnarData
+import space.kscience.kmath.functions.PiecewisePolynomial
+import space.kscience.kmath.functions.ListPolynomial
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.Field
+import space.kscience.kmath.operations.invoke
+
+@OptIn(UnstableKMathAPI::class)
+internal fun <T : Comparable<T>> insureSorted(points: XYColumnarData<*, T, *>) {
+    for (i in 0 until points.size - 1)
+        require(points.x[i + 1] > points.x[i]) { "Input data is not sorted at index $i" }
+}
+
+/**
+ * Reference JVM implementation: https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/analysis/interpolation/LinearInterpolator.java
+ */
+public class LinearInterpolator<T : Comparable<T>>(override val algebra: Field<T>) : PolynomialInterpolator<T> {
+
+    @OptIn(UnstableKMathAPI::class)
+    override fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T> = algebra {
+        require(points.size > 0) { "Point array should not be empty" }
+        insureSorted(points)
+
+        PiecewisePolynomial(points.x[0]) {
+            for (i in 0 until points.size - 1) {
+                val slope = (points.y[i + 1] - points.y[i]) / (points.x[i + 1] - points.x[i])
+                val const = points.y[i] - slope * points.x[i]
+                val polynomial = ListPolynomial(const, slope)
+                putRight(points.x[i + 1], polynomial)
+            }
+        }
+    }
+}
+
+public val <T : Comparable<T>> Field<T>.linearInterpolator: LinearInterpolator<T>
+    get() = LinearInterpolator(this)
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt
new file mode 100644
index 000000000..0bcbfd0c6
--- /dev/null
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt
@@ -0,0 +1,88 @@
+/*
+ * 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.
+ */
+
+/*
+ * 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.interpolation
+
+import space.kscience.kmath.data.XYColumnarData
+import space.kscience.kmath.functions.PiecewisePolynomial
+import space.kscience.kmath.functions.ListPolynomial
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.Field
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.MutableBufferFactory
+
+/**
+ * Generic spline interpolator. Not recommended for performance critical places, use platform-specific and type
+ * specific ones.
+ *
+ * Based on
+ * https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/SplineInterpolator.java
+ */
+public class SplineInterpolator<T : Comparable<T>>(
+    override val algebra: Field<T>,
+    public val bufferFactory: MutableBufferFactory<T>,
+) : PolynomialInterpolator<T> {
+    //TODO possibly optimize zeroed buffers
+
+    @OptIn(UnstableKMathAPI::class)
+    override fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T> = algebra {
+        require(points.size >= 3) { "Can't use spline interpolator with less than 3 points" }
+        insureSorted(points)
+        // Number of intervals.  The number of data points is n + 1.
+        val n = points.size - 1
+        // Differences between knot points
+        val h = bufferFactory(n) { i -> points.x[i + 1] - points.x[i] }
+        val mu = bufferFactory(n) { zero }
+        val z = bufferFactory(n + 1) { zero }
+
+        for (i in 1 until n) {
+            val g = 2.0 * (points.x[i + 1] - points.x[i - 1]) - h[i - 1] * mu[i - 1]
+            mu[i] = h[i] / g
+            z[i] =
+                ((points.y[i + 1] * h[i - 1] - points.y[i] * (points.x[i + 1] - points.x[i - 1]) + points.y[i - 1] * h[i]) * 3.0 /
+                        (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g
+        }
+
+        // cubic spline coefficients --  b is linear, c quadratic, d is cubic (original y's are constants)
+
+        PiecewisePolynomial(points.x[points.size - 1]) {
+            var cOld = zero
+
+            for (j in n - 1 downTo 0) {
+                val c = z[j] - mu[j] * cOld
+                val a = points.y[j]
+                val b = (points.y[j + 1] - points.y[j]) / h[j] - h[j] * (cOld + 2.0 * c) / 3.0
+                val d = (cOld - c) / (3.0 * h[j])
+                val x0 = points.x[j]
+                val x02 = x0 * x0
+                val x03 = x02 * x0
+                //Shift coefficients to represent absolute polynomial instead of one with an offset
+                val polynomial = ListPolynomial(
+                    a - b * x0 + c * x02 - d * x03,
+                    b - 2 * c * x0 + 3 * d * x02,
+                    c - 3 * d * x0,
+                    d
+                )
+                cOld = c
+                putLeft(x0, polynomial)
+            }
+        }
+    }
+}
+
+
+public fun <T : Comparable<T>> Field<T>.splineInterpolator(
+    bufferFactory: MutableBufferFactory<T>,
+): SplineInterpolator<T> = SplineInterpolator(this, bufferFactory)
+
+public val DoubleField.splineInterpolator: SplineInterpolator<Double>
+    get() = SplineInterpolator(this, ::DoubleBuffer)
\ No newline at end of file
diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SplineIntegralTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SplineIntegralTest.kt
new file mode 100644
index 000000000..aae0ad017
--- /dev/null
+++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SplineIntegralTest.kt
@@ -0,0 +1,48 @@
+/*
+ * 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.integration
+
+import space.kscience.kmath.functions.ListPolynomial
+import space.kscience.kmath.functions.integrate
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.DoubleField
+import kotlin.math.PI
+import kotlin.math.sin
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+@OptIn(UnstableKMathAPI::class)
+class SplineIntegralTest {
+
+    @Test
+    fun integratePolynomial(){
+        val polynomial = ListPolynomial(1.0, 2.0, 3.0)
+        val integral = polynomial.integrate(DoubleField,1.0..2.0)
+        assertEquals(11.0, integral, 0.001)
+    }
+
+    @Test
+    fun gaussSin() {
+        val res = DoubleField.splineIntegrator.integrate(0.0..2 * PI, IntegrandMaxCalls(5)) { x ->
+            sin(x)
+        }
+        assertEquals(0.0, res.value, 1e-2)
+    }
+
+    @Test
+    fun gaussUniform() {
+        val res = DoubleField.splineIntegrator.integrate(35.0..100.0, IntegrandMaxCalls(20)) { x ->
+            if(x in 30.0..50.0){
+                1.0
+            } else {
+                0.0
+            }
+        }
+        assertEquals(15.0, res.value, 0.5)
+    }
+
+
+}
\ No newline at end of file
diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/LinearInterpolatorTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/LinearInterpolatorTest.kt
new file mode 100644
index 000000000..1143036d4
--- /dev/null
+++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/LinearInterpolatorTest.kt
@@ -0,0 +1,29 @@
+/*
+ * 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.interpolation
+
+import space.kscience.kmath.operations.DoubleField
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+internal class LinearInterpolatorTest {
+    @Test
+    fun testInterpolation() {
+        val data = listOf(
+            0.0 to 0.0,
+            1.0 to 1.0,
+            2.0 to 3.0,
+            3.0 to 4.0
+        )
+
+        //val polynomial: PiecewisePolynomial<Double> = DoubleField.linearInterpolator.interpolatePolynomials(data)
+        val function = DoubleField.linearInterpolator.interpolate(data)
+        assertEquals(null, function(-1.0))
+        assertEquals(0.5, function(0.5))
+        assertEquals(2.0, function(1.5))
+        assertEquals(3.0, function(2.0))
+    }
+}
diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/SplineInterpolatorTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/SplineInterpolatorTest.kt
new file mode 100644
index 000000000..f748535e2
--- /dev/null
+++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/SplineInterpolatorTest.kt
@@ -0,0 +1,36 @@
+/*
+ * 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.
+ */
+
+/*
+ * 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.interpolation
+
+import space.kscience.kmath.operations.DoubleField
+import kotlin.math.PI
+import kotlin.math.sin
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+internal class SplineInterpolatorTest {
+    @Test
+    fun testInterpolation() {
+        val data = (0..10).map {
+            val x = it.toDouble() / 5 * PI
+            x to sin(x)
+        }
+
+        //val polynomial: PiecewisePolynomial<Double> = DoubleField.splineInterpolator.interpolatePolynomials(data)
+
+        val function = DoubleField.splineInterpolator.interpolate(data, Double.NaN)
+
+        assertEquals(Double.NaN, function(-1.0))
+        assertEquals(sin(0.5), function(0.5), 0.1)
+        assertEquals(sin(1.5), function(1.5), 0.1)
+        assertEquals(sin(2.0), function(2.0), 0.1)
+    }
+}