diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt
new file mode 100644
index 000000000..8511ed66e
--- /dev/null
+++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt
@@ -0,0 +1,92 @@
+package space.kscience.kmath.commons.integration
+
+import org.apache.commons.math3.analysis.integration.IterativeLegendreGaussIntegrator
+import org.apache.commons.math3.analysis.integration.SimpsonIntegrator
+import space.kscience.kmath.integration.*
+import space.kscience.kmath.misc.UnstableKMathAPI
+
+/**
+ * Integration wrapper for Common-maths UnivariateIntegrator
+ */
+public class CMIntegrator(
+    private val defaultMaxCalls: Int = 200,
+    public val integratorBuilder: (Integrand) -> org.apache.commons.math3.analysis.integration.UnivariateIntegrator,
+) : UnivariateIntegrator<Double> {
+
+    public class TargetRelativeAccuracy(public val value: Double) : IntegrandFeature
+    public class TargetAbsoluteAccuracy(public val value: Double) : IntegrandFeature
+
+    public class MinIterations(public val value: Int) : IntegrandFeature
+    public class MaxIterations(public val value: Int) : IntegrandFeature
+
+    override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
+        val integrator = integratorBuilder(integrand)
+        val maxCalls = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: defaultMaxCalls
+        val remainingCalls = maxCalls - integrand.calls
+        val range = integrand.getFeature<IntegrationRange<Double>>()?.range
+            ?: error("Integration range is not provided")
+        val res = integrator.integrate(remainingCalls, integrand.function, range.start, range.endInclusive)
+
+        return integrand +
+                IntegrandValue(res) +
+                IntegrandAbsoluteAccuracy(integrator.absoluteAccuracy) +
+                IntegrandRelativeAccuracy(integrator.relativeAccuracy) +
+                IntegrandCalls(integrator.evaluations + integrand.calls)
+    }
+
+
+    public companion object {
+        /**
+         * Create a Simpson integrator based on [SimpsonIntegrator]
+         */
+        public fun simpson(defaultMaxCalls: Int = 200): CMIntegrator = CMIntegrator(defaultMaxCalls) { integrand ->
+            val absoluteAccuracy = integrand.getFeature<TargetAbsoluteAccuracy>()?.value
+                ?: SimpsonIntegrator.DEFAULT_ABSOLUTE_ACCURACY
+            val relativeAccuracy = integrand.getFeature<TargetRelativeAccuracy>()?.value
+                ?: SimpsonIntegrator.DEFAULT_ABSOLUTE_ACCURACY
+            val minIterations = integrand.getFeature<MinIterations>()?.value
+                ?: SimpsonIntegrator.DEFAULT_MIN_ITERATIONS_COUNT
+            val maxIterations = integrand.getFeature<MaxIterations>()?.value
+                ?: SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT
+
+            SimpsonIntegrator(relativeAccuracy, absoluteAccuracy, minIterations, maxIterations)
+        }
+
+        /**
+         * Create a Gauss-Legandre integrator based on [IterativeLegendreGaussIntegrator]
+         */
+        public fun legandre(numPoints: Int, defaultMaxCalls: Int = numPoints * 5): CMIntegrator =
+            CMIntegrator(defaultMaxCalls) { integrand ->
+                val absoluteAccuracy = integrand.getFeature<TargetAbsoluteAccuracy>()?.value
+                    ?: IterativeLegendreGaussIntegrator.DEFAULT_ABSOLUTE_ACCURACY
+                val relativeAccuracy = integrand.getFeature<TargetRelativeAccuracy>()?.value
+                    ?: IterativeLegendreGaussIntegrator.DEFAULT_ABSOLUTE_ACCURACY
+                val minIterations = integrand.getFeature<MinIterations>()?.value
+                    ?: IterativeLegendreGaussIntegrator.DEFAULT_MIN_ITERATIONS_COUNT
+                val maxIterations = integrand.getFeature<MaxIterations>()?.value
+                    ?: IterativeLegendreGaussIntegrator.DEFAULT_MAX_ITERATIONS_COUNT
+
+                IterativeLegendreGaussIntegrator(
+                    numPoints,
+                    relativeAccuracy,
+                    absoluteAccuracy,
+                    minIterations,
+                    maxIterations
+                )
+            }
+    }
+}
+
+@UnstableKMathAPI
+public var MutableList<IntegrandFeature>.targetAbsoluteAccuracy: Double?
+    get() = filterIsInstance<CMIntegrator.TargetAbsoluteAccuracy>().lastOrNull()?.value
+    set(value){
+        value?.let { add(CMIntegrator.TargetAbsoluteAccuracy(value))}
+    }
+
+@UnstableKMathAPI
+public var MutableList<IntegrandFeature>.targetRelativeAccuracy: Double?
+    get() = filterIsInstance<CMIntegrator.TargetRelativeAccuracy>().lastOrNull()?.value
+    set(value){
+        value?.let { add(CMIntegrator.TargetRelativeAccuracy(value))}
+    }
\ No newline at end of file
diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt
new file mode 100644
index 000000000..5a18756ac
--- /dev/null
+++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt
@@ -0,0 +1,94 @@
+/* 
+ * Copyright 2015 Alexander Nozik.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package space.kscience.kmath.commons.integration
+
+import org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator
+import org.apache.commons.math3.analysis.integration.gauss.GaussIntegratorFactory
+import space.kscience.kmath.integration.*
+
+/**
+ * A simple one-pass integrator based on Gauss rule
+ */
+public class GaussRuleIntegrator(
+    private val numpoints: Int,
+    private var type: GaussRule = GaussRule.LEGANDRE,
+) : UnivariateIntegrator<Double> {
+
+    override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
+        val range = integrand.getFeature<IntegrationRange<Double>>()?.range
+            ?: error("Integration range is not provided")
+        val integrator: GaussIntegrator = getIntegrator(range)
+        //TODO check performance
+        val res: Double = integrator.integrate(integrand.function)
+        return integrand + IntegrandValue(res) + IntegrandCalls(integrand.calls + numpoints)
+    }
+
+    private fun getIntegrator(range: ClosedRange<Double>): GaussIntegrator {
+        return when (type) {
+            GaussRule.LEGANDRE -> factory.legendre(
+                numpoints,
+                range.start,
+                range.endInclusive
+            )
+            GaussRule.LEGANDREHP -> factory.legendreHighPrecision(
+                numpoints,
+                range.start,
+                range.endInclusive
+            )
+            GaussRule.UNIFORM -> GaussIntegrator(
+                getUniformRule(
+                    range.start,
+                    range.endInclusive,
+                    numpoints
+                )
+            )
+        }
+    }
+
+    private fun getUniformRule(
+        min: Double,
+        max: Double,
+        numPoints: Int,
+    ): org.apache.commons.math3.util.Pair<DoubleArray, DoubleArray> {
+        assert(numPoints > 2)
+        val points = DoubleArray(numPoints)
+        val weights = DoubleArray(numPoints)
+        val step = (max - min) / (numPoints - 1)
+        points[0] = min
+        for (i in 1 until numPoints) {
+            points[i] = points[i - 1] + step
+            weights[i] = step
+        }
+        return org.apache.commons.math3.util.Pair<DoubleArray, DoubleArray>(points, weights)
+    }
+
+    public enum class GaussRule {
+        UNIFORM, LEGANDRE, LEGANDREHP
+    }
+
+    public companion object {
+        private val factory: GaussIntegratorFactory = GaussIntegratorFactory()
+
+        public fun integrate(
+            range: ClosedRange<Double>,
+            numPoints: Int = 100,
+            type: GaussRule = GaussRule.LEGANDRE,
+            function: (Double) -> Double,
+        ): Double = GaussRuleIntegrator(numPoints, type).integrate(
+            UnivariateIntegrand(function, IntegrationRange(range))
+        ).value!!
+    }
+}
\ No newline at end of file
diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizationProblem.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizationProblem.kt
index d655b4f61..13a10475f 100644
--- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizationProblem.kt
+++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizationProblem.kt
@@ -19,7 +19,7 @@ import kotlin.reflect.KClass
 public operator fun PointValuePair.component1(): DoubleArray = point
 public operator fun PointValuePair.component2(): Double = value
 
-public class CMOptimizationProblem(override val symbols: List<Symbol>, ) :
+public class CMOptimizationProblem(override val symbols: List<Symbol>) :
     OptimizationProblem<Double>, SymbolIndexer, OptimizationFeature {
     private val optimizationData: HashMap<KClass<out OptimizationData>, OptimizationData> = HashMap()
     private var optimizatorBuilder: (() -> MultivariateOptimizer)? = null
diff --git a/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/integration/IntegrationTest.kt b/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/integration/IntegrationTest.kt
new file mode 100644
index 000000000..3693d5796
--- /dev/null
+++ b/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/integration/IntegrationTest.kt
@@ -0,0 +1,30 @@
+package space.kscience.kmath.commons.integration
+
+import org.junit.jupiter.api.Test
+import space.kscience.kmath.integration.integrate
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.RealField.sin
+import kotlin.math.PI
+import kotlin.math.abs
+import kotlin.test.assertTrue
+
+@UnstableKMathAPI
+internal class IntegrationTest {
+    private val function: (Double) -> Double = { sin(it) }
+
+    @Test
+    fun simpson() {
+        val res = CMIntegrator.simpson().integrate(0.0..2 * PI, function)
+        assertTrue { abs(res) < 1e-3 }
+    }
+
+    @Test
+    fun customSimpson() {
+        val res = CMIntegrator.simpson().integrate(0.0..PI, function) {
+            targetRelativeAccuracy = 1e-4
+            targetAbsoluteAccuracy = 1e-4
+        }
+        assertTrue { abs(res - 2) < 1e-3 }
+        assertTrue { abs(res - 2) > 1e-12 }
+    }
+}
\ No newline at end of file
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrand.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrand.kt
new file mode 100644
index 000000000..ae964b271
--- /dev/null
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrand.kt
@@ -0,0 +1,22 @@
+package space.kscience.kmath.integration
+
+import kotlin.reflect.KClass
+
+public interface IntegrandFeature
+
+public interface Integrand {
+    public fun <T : IntegrandFeature> getFeature(type: KClass<T>): T?
+}
+
+public inline fun <reified T : IntegrandFeature> Integrand.getFeature(): T? = getFeature(T::class)
+
+public class IntegrandValue<T : Any>(public val value: T) : IntegrandFeature
+
+public class IntegrandRelativeAccuracy(public val accuracy: Double) : IntegrandFeature
+
+public class IntegrandAbsoluteAccuracy(public val accuracy: Double) : IntegrandFeature
+
+public class IntegrandCalls(public val calls: Int) : IntegrandFeature
+public val Integrand.calls: Int get() = getFeature<IntegrandCalls>()?.calls ?: 0
+
+public class IntegrandMaxCalls(public val maxCalls: Int) : IntegrandFeature
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt
new file mode 100644
index 000000000..7027e62c7
--- /dev/null
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt
@@ -0,0 +1,11 @@
+package space.kscience.kmath.integration
+
+/**
+ * A general interface for all integrators
+ */
+public interface Integrator<I: Integrand> {
+    /**
+     * Run one integration pass and return a new [Integrand] with a new set of features
+     */
+    public fun integrate(integrand: I): I
+}
\ No newline at end of file
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/MultivariateIntegrand.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/MultivariateIntegrand.kt
new file mode 100644
index 000000000..4b3a6deb4
--- /dev/null
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/MultivariateIntegrand.kt
@@ -0,0 +1,27 @@
+package space.kscience.kmath.integration
+
+import space.kscience.kmath.linear.Point
+import kotlin.reflect.KClass
+
+public class MultivariateIntegrand<T : Any> internal constructor(
+    private val features: Map<KClass<*>, IntegrandFeature>,
+    public val function: (Point<T>) -> T,
+) : Integrand {
+
+    @Suppress("UNCHECKED_CAST")
+    override fun <T : IntegrandFeature> getFeature(type: KClass<T>): T? = features[type] as? T
+
+    public operator fun <F : IntegrandFeature> plus(pair: Pair<KClass<out F>, F>): MultivariateIntegrand<T> =
+        MultivariateIntegrand(features + pair, function)
+
+    public operator fun <F : IntegrandFeature> plus(feature: F): MultivariateIntegrand<T> =
+        plus(feature::class to feature)
+}
+
+@Suppress("FunctionName")
+public fun <T : Any> MultivariateIntegrand(
+    vararg features: IntegrandFeature,
+    function: (Point<T>) -> T,
+): MultivariateIntegrand<T> = MultivariateIntegrand(features.associateBy { it::class }, function)
+
+public val <T : Any> MultivariateIntegrand<T>.value: T? get() = getFeature<IntegrandValue<T>>()?.value
diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt
new file mode 100644
index 000000000..9389318e8
--- /dev/null
+++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt
@@ -0,0 +1,64 @@
+package space.kscience.kmath.integration
+
+import space.kscience.kmath.misc.UnstableKMathAPI
+import kotlin.reflect.KClass
+
+public class UnivariateIntegrand<T : Any> internal constructor(
+    private val features: Map<KClass<*>, IntegrandFeature>,
+    public val function: (T) -> T,
+) : Integrand {
+
+    @Suppress("UNCHECKED_CAST")
+    override fun <T : IntegrandFeature> getFeature(type: KClass<T>): T? = features[type] as? T
+
+    public operator fun <F : IntegrandFeature> plus(pair: Pair<KClass<out F>, F>): UnivariateIntegrand<T> =
+        UnivariateIntegrand(features + pair, function)
+
+    public operator fun <F : IntegrandFeature> plus(feature: F): UnivariateIntegrand<T> =
+        plus(feature::class to feature)
+}
+
+@Suppress("FunctionName")
+public fun <T : Any> UnivariateIntegrand(
+    function: (T) -> T,
+    vararg features: IntegrandFeature,
+): UnivariateIntegrand<T> = UnivariateIntegrand(features.associateBy { it::class }, function)
+
+public typealias UnivariateIntegrator<T> = Integrator<UnivariateIntegrand<T>>
+
+public inline class IntegrationRange<T : Comparable<T>>(public val range: ClosedRange<T>) : IntegrandFeature
+
+public val <T : Any> UnivariateIntegrand<T>.value: T? get() = getFeature<IntegrandValue<T>>()?.value
+
+/**
+ * A shortcut method to integrate a [function] in [range] with additional [features].
+ * The [function] is placed in the end position to allow passing a lambda.
+ */
+@UnstableKMathAPI
+public fun UnivariateIntegrator<Double>.integrate(
+    range: ClosedRange<Double>,
+    vararg features: IntegrandFeature,
+    function: (Double) -> Double,
+): Double = integrate(
+    UnivariateIntegrand(function, IntegrationRange(range), *features)
+).value ?: error("Unexpected: no value after integration.")
+
+/**
+ * A shortcut method to integrate a [function] in [range] with additional [features].
+ * The [function] is placed in the end position to allow passing a lambda.
+ */
+@UnstableKMathAPI
+public fun UnivariateIntegrator<Double>.integrate(
+    range: ClosedRange<Double>,
+    function: (Double) -> Double,
+    featureBuilder: (MutableList<IntegrandFeature>.() -> Unit) = {},
+): Double {
+    //TODO use dedicated feature builder class instead or add extensions to MutableList<IntegrandFeature>
+    val features = buildList {
+        featureBuilder()
+        add(IntegrationRange(range))
+    }
+    return integrate(
+        UnivariateIntegrand(function, *features.toTypedArray())
+    ).value ?: error("Unexpected: no value after integration.")
+}
\ No newline at end of file