diff --git a/gradle.properties b/gradle.properties
index c3f070c2d..048c9c9f3 100644
--- a/gradle.properties
+++ b/gradle.properties
@@ -9,7 +9,7 @@ kotlin.native.ignoreDisabledTargets=true
 org.gradle.configureondemand=true
 org.gradle.jvmargs=-Xmx4096m
 
-toolsVersion=0.14.2-kotlin-1.8.10
+toolsVersion=0.14.3-kotlin-1.8.20-RC
 
 
 org.gradle.parallel=true
diff --git a/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/Euclidean3DSpace.kt b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/Euclidean3DSpace.kt
index cc641a3f1..38c252bc0 100644
--- a/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/Euclidean3DSpace.kt
+++ b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/Euclidean3DSpace.kt
@@ -78,8 +78,11 @@ public object Euclidean3DSpace : GeometrySpace<DoubleVector3D>, ScaleOperations<
         }
     }
 
+    public fun vector(x: Double, y: Double, z: Double): DoubleVector3D =
+        Vector3DImpl(x, y, z)
+
     public fun vector(x: Number, y: Number, z: Number): DoubleVector3D =
-        Vector3DImpl(x.toDouble(), y.toDouble(), z.toDouble())
+        vector(x.toDouble(), y.toDouble(), z.toDouble())
 
     override val zero: DoubleVector3D by lazy { vector(0.0, 0.0, 0.0) }
 
@@ -100,6 +103,24 @@ public object Euclidean3DSpace : GeometrySpace<DoubleVector3D>, ScaleOperations<
     override fun DoubleVector3D.dot(other: DoubleVector3D): Double =
         x * other.x + y * other.y + z * other.z
 
+    /**
+     * Compute vector product of [first] and [second]. The basis assumed to be right-handed.
+     */
+    public fun vectorProduct(
+        first: DoubleVector3D,
+        second: DoubleVector3D,
+    ): DoubleVector3D {
+        val (x1, y1, z1) = first
+        val (x2, y2, z2) = second
+
+        return vector(y1 * z2 - y2 * z2, z1 * x2 - z2 * x2, x1 * y2 - x2 * y2)
+    }
+
+    /**
+     * Vector product with a right basis
+     */
+    public infix fun DoubleVector3D.cross(other: DoubleVector3D): Vector3D<Double> = vectorProduct(this, other)
+
     public val xAxis: DoubleVector3D = vector(1.0, 0.0, 0.0)
     public val yAxis: DoubleVector3D = vector(0.0, 1.0, 0.0)
     public val zAxis: DoubleVector3D = vector(0.0, 0.0, 1.0)
diff --git a/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Euclidean3DSpaceTest.kt b/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Euclidean3DSpaceTest.kt
index 6d9a169eb..291b0ad47 100644
--- a/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Euclidean3DSpaceTest.kt
+++ b/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Euclidean3DSpaceTest.kt
@@ -57,23 +57,39 @@ internal class Euclidean3DSpaceTest {
     }
 
     @Test
-    fun add() {
-        with(Euclidean3DSpace) {
-            assertVectorEquals(
-                vector(1.0, -2.0, 0.001),
-                vector(1.0, -2.0, 0.001) + zero
-            )
-            assertVectorEquals(
-                vector(8.0, -3.0, 3.001),
-                vector(1.0, 2.0, 3.0) + vector(7.0, -5.0, 0.001)
-            )
-        }
+    fun add() = with(Euclidean3DSpace) {
+        assertVectorEquals(
+            vector(1.0, -2.0, 0.001),
+            vector(1.0, -2.0, 0.001) + zero
+        )
+        assertVectorEquals(
+            vector(8.0, -3.0, 3.001),
+            vector(1.0, 2.0, 3.0) + vector(7.0, -5.0, 0.001)
+        )
     }
 
     @Test
-    fun multiply() {
-        with(Euclidean3DSpace) {
-            assertVectorEquals(vector(2.0, -4.0, 0.0), vector(1.0, -2.0, 0.0) * 2)
-        }
+    fun multiply() = with(Euclidean3DSpace) {
+        assertVectorEquals(vector(2.0, -4.0, 0.0), vector(1.0, -2.0, 0.0) * 2)
     }
+
+    @Test
+    fun vectorProduct() = with(Euclidean3DSpace) {
+        assertVectorEquals(zAxis, vectorProduct(xAxis, yAxis))
+        assertVectorEquals(zAxis, xAxis cross yAxis)
+        assertVectorEquals(-zAxis, vectorProduct(yAxis, xAxis))
+        assertVectorEquals(zAxis, vectorProduct(yAxis, xAxis, rightBasis = false))
+    }
+
+    @Test
+    fun doubleVectorProduct() = with(Euclidean3DSpace) {
+        val a = vector(1, 2, -3)
+        val b = vector(-1, 0, 1)
+        val c = vector(4, 5, 6)
+
+        val res = a cross (b cross c)
+        val expected = b * (a dot c) - c * (a dot b)
+        assertVectorEquals(expected, res)
+    }
+
 }