diff --git a/source/simple/geom.hpp b/source/simple/geom.hpp
index af3a73455b6bc4a6b449eb31bd7b1a57e100cdf3..5a822d8869eb3415ac4d305eb7d41fda30fc7769 100644
--- a/source/simple/geom.hpp
+++ b/source/simple/geom.hpp
@@ -1,2 +1,4 @@
-#include "geom/vector.hpp"
+#include "geom/algorithm.h"
+#include "geom/algorithm.hpp"
 #include "geom/segment.hpp"
+#include "geom/vector.hpp"
diff --git a/source/simple/geom/algorithm.h b/source/simple/geom/algorithm.h
new file mode 100644
index 0000000000000000000000000000000000000000..431d92c84e0bf41737274ab89edd6ddd3ef6843e
--- /dev/null
+++ b/source/simple/geom/algorithm.h
@@ -0,0 +1,27 @@
+#ifndef SIMPLE_GEOM_ALGORITHM_H
+#define SIMPLE_GEOM_ALGORITHM_H
+#include "vector.hpp"
+
+namespace simple::geom
+{
+	template <typename Coordinate, size_t Dimensions, typename Order, typename Function, size_t LoopDimesntions = Dimensions>
+	constexpr void loop_on
+	(
+		vector<Coordinate, Dimensions, Order>& index,
+		const vector<Coordinate, Dimensions, Order>& lower,
+		const vector<Coordinate, Dimensions, Order>& upper,
+		const vector<Coordinate, Dimensions, Order>& step,
+		Function&& callback
+	);
+
+	template <typename Coordinate, size_t Dimensions, typename Order, typename Function>
+	constexpr void loop
+	(
+		const vector<Coordinate, Dimensions, Order>& lower,
+		const vector<Coordinate, Dimensions, Order>& upper,
+		const vector<Coordinate, Dimensions, Order>& step,
+		Function&& callback
+	);
+} // namespace simple::geom
+
+#endif /* end of include guard */
diff --git a/source/simple/geom/algorithm.hpp b/source/simple/geom/algorithm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..862a33258d354cc4acd0025cd5cc99d2b6890bdd
--- /dev/null
+++ b/source/simple/geom/algorithm.hpp
@@ -0,0 +1,55 @@
+#ifndef SIMPLE_GEOM_ALGORITHM_HPP
+#define SIMPLE_GEOM_ALGORITHM_HPP
+#include "algorithm.h"
+
+#include <functional>
+
+namespace simple::geom
+{
+
+	template <typename Coordinate, size_t Dimensions, typename Order, typename Function, size_t LoopDimesntions = Dimensions>
+	constexpr void loop_on
+	(
+		vector<Coordinate, Dimensions, Order>& index,
+		const vector<Coordinate, Dimensions, Order>& lower,
+		const vector<Coordinate, Dimensions, Order>& upper,
+		const vector<Coordinate, Dimensions, Order>& step,
+		Function&& callback
+	)
+	{
+		if constexpr (LoopDimesntions == 0)
+			std::invoke(std::forward<Function>(callback),index);
+		else
+		{
+			constexpr size_t I = LoopDimesntions-1;
+			for
+			(
+				index[I] = lower[I];
+				index[I] < upper[I];
+				index[I] += step[I]
+			)
+			{
+				loop_on<Coordinate,Dimensions,Order,Function,I>
+					(index, lower, upper, step,
+					 std::forward<Function>(callback));
+			}
+		}
+	};
+
+	template <typename Coordinate, size_t Dimensions, typename Order, typename Function>
+	constexpr void loop
+	(
+		const vector<Coordinate, Dimensions, Order>& lower,
+		const vector<Coordinate, Dimensions, Order>& upper,
+		const vector<Coordinate, Dimensions, Order>& step,
+		Function&& callback
+	)
+	{
+		vector<Coordinate, Dimensions, Order> index{};
+		loop_on(index, lower, upper, step,
+			 std::forward<Function>(callback));
+	}
+
+} // namespace simple::geom
+
+#endif /* end of include guard */
diff --git a/unit_tests/algorithm.cpp b/unit_tests/algorithm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6cf5c38dbcc49870b29660a0eaf2992df0d4d5fb
--- /dev/null
+++ b/unit_tests/algorithm.cpp
@@ -0,0 +1,138 @@
+#include "simple/geom/algorithm.hpp"
+#include <cmath>
+#include <cassert>
+#include <vector>
+
+using namespace simple;
+using geom::vector;
+using float4 = vector<float, 4>;
+using int4 = vector<int, 4>;
+
+void Basic()
+{
+
+	int4 p (1, 2, 3, 4);
+	int4 p2 (4, 3, 2, 1);
+	assert( int4(1, 2, 2, 1) == geom::min(p, p2) );
+	assert( int4(4, 3, 3, 4) == geom::max(p, p2) );
+
+	int4 p3;
+	(p3 = p).min(p2);
+	assert( int4(1, 2, 2, 1) == p3 );
+	(p3 = p).max(p2);
+	assert( int4(4, 3, 3, 4) == p3 );
+
+	p = -p;
+	assert( int4(-1, 2, 2, -3) == geom::clamp(int4(-10, 2, 3, -3), p, p2) );
+	assert( int4(0, -2, -1, 1) == geom::clamp(int4(0, -3, -1, 10), p, p2) );
+	assert( int4(3, -1, -2, 0) == geom::clamp(int4(3, -1, -2, 0), p, p2) );
+	assert( int4(-1, 3, -3, 1) == geom::clamp(int4(-3, 7, -5, 3), p, p2) );
+
+	(p3 = int4(-10, 2, 3, -3)).clamp(p,p2);
+	assert( int4(-1, 2, 2, -3) == p3 );
+	(p3 = int4(0, -3, -1, 10)).clamp(p,p2);
+	assert( int4(0, -2, -1, 1) == p3);
+	(p3 = int4(3, -1, -2, 0)).clamp(p,p2);
+	assert( int4(3, -1, -2, 0) == p3);
+	(p3 = int4(-3, 7, -5, 3)).clamp(p,p2);
+	assert( int4(-1, 3, -3, 1) == p3);
+
+	float4 fp (1.1f, 3.4f, 4.5f, 8.9f);
+	assert( float4(1.f, 3.f, 4.f, 8.f) == trunc(fp) );
+	assert( float4(1.f, 3.f, 4.f, 8.f) == floor(fp) );
+	assert( float4(2.f, 4.f, 5.f, 9.f) == ceil(fp) );
+	assert( float4(1.f, 3.f, 5.f, 9.f) == round(fp) );
+
+	float4 fp2;
+	(fp2 = fp).trunc();
+	assert( float4(1.f, 3.f, 4.f, 8.f) == fp2 );
+	(fp2 = fp).floor();
+	assert( float4(1.f, 3.f, 4.f, 8.f) == fp2 );
+	(fp2 = fp).ceil();
+	assert( float4(2.f, 4.f, 5.f, 9.f) == fp2 );
+	(fp2 = fp).round();
+	assert( float4(1.f, 3.f, 5.f, 9.f) == fp2 );
+
+	fp = float4(1.1f, -3.4f, 4.5f, -8.9f);
+	assert( float4(1.f, -3.f, 4.f, -8.f) == trunc(fp) );
+	assert( float4(1.f, -4.f, 4.f, -9.f) == floor(fp) );
+	assert( float4(2.f, -3.f, 5.f, -8.f) == ceil(fp) );
+	assert( float4(1.f, -3.f, 5.f, -9.f) == round(fp) );
+
+	(fp2 = fp).trunc();
+	assert( float4(1.f, -3.f, 4.f, -8.f) == fp2 );
+	(fp2 = fp).floor();
+	assert( float4(1.f, -4.f, 4.f, -9.f) == fp2 );
+	(fp2 = fp).ceil();
+	assert( float4(2.f, -3.f, 5.f, -8.f) == fp2 );
+	(fp2 = fp).round();
+	assert( float4(1.f, -3.f, 5.f, -9.f) == fp2 );
+
+	int4 p4 (2, 6, 3, 0);
+	assert( p4.magnitude() ==  49);
+	assert( p4.quadrance() ==  49);
+	assert( p4.length() == 7 );
+	assert( magnitude(p4) ==  49);
+	assert( quadrance(p4) ==  49);
+	assert( length(p4) == 7 );
+
+
+	assert( signum(int4(-12, 34, 0, 1)) == int4(-1, 1, 0, 1) );
+	assert( signum(int4(9, -13, -77, 0)) == int4(1, -1, -1, 0) );
+	assert( signum(float4(-.3f, .12f, .0f, +.0f)) == float4(-1.f, 1.f, 0.f, 0.f) );
+
+	(p3 = int4(-12, 34, 0, 1)).signum();
+	assert( int4(-1, 1, 0, 1) == p3 );
+	(p3 = int4(9, -13, -77, 0)).signum();
+	assert( int4(1, -1, -1, 0) == p3 );
+	(fp2 = float4(-.3f, .12f, .0f, +.0f)).signum();
+	assert( float4(-1.f, 1.f, 0.f, 0.f) == fp2 );
+}
+
+void MultidimentionalIteration()
+{
+	using int3 = vector<int, 3>;
+	const auto lower = int3{13,3,-20};
+	const auto upper = int3{45,32,12};
+
+	std::vector<int3> test_data;
+	test_data.reserve(10000);
+	std::vector<int3> data;
+	data.reserve(10000);
+
+
+	for(int k = lower[2]; k < upper[2]; ++k)
+		for(int j = lower[1]; j < upper[1]; ++j)
+			for(int i = lower[0]; i < upper[0]; ++i)
+				test_data.push_back({i,j,k});
+
+	loop(lower, upper, int3::one(), [&data](auto& i)
+	{
+		data.push_back(i);
+	});
+
+	assert(data == test_data);
+
+	test_data.clear();
+	data.clear();
+
+	auto step = int3{1,2,3};
+
+	for(int k = lower[2]; k < upper[2]; k += step[2])
+		for(int j = lower[1]; j < upper[1]; j += step[1])
+			for(int i = lower[0]; i < upper[0]; i += step[0])
+				test_data.push_back({i,j,k});
+
+	loop(lower, upper, step, [&data](auto& i)
+	{
+		data.push_back(i);
+	});
+
+	assert(data == test_data);
+}
+
+int main()
+{
+	Basic();
+	MultidimentionalIteration();
+}
diff --git a/unit_tests/point.cpp b/unit_tests/point.cpp
index 1008913f3184ac70c8aa059c2e16b4d0645be8e9..130f9c535411f8116ec38849ec7f896bcc5a9d83 100644
--- a/unit_tests/point.cpp
+++ b/unit_tests/point.cpp
@@ -415,87 +415,6 @@ void ComparisonOperators()
 	assert(p * int4(p2 == int4::zero() | p2 == int4::one()) == int4(1, 0, 3, 4));
 }
 
-void Algorithms()
-{
-
-	int4 p (1, 2, 3, 4);
-	int4 p2 (4, 3, 2, 1);
-	assert( int4(1, 2, 2, 1) == geom::min(p, p2) );
-	assert( int4(4, 3, 3, 4) == geom::max(p, p2) );
-
-	int4 p3;
-	(p3 = p).min(p2);
-	assert( int4(1, 2, 2, 1) == p3 );
-	(p3 = p).max(p2);
-	assert( int4(4, 3, 3, 4) == p3 );
-
-	p = -p;
-	assert( int4(-1, 2, 2, -3) == geom::clamp(int4(-10, 2, 3, -3), p, p2) );
-	assert( int4(0, -2, -1, 1) == geom::clamp(int4(0, -3, -1, 10), p, p2) );
-	assert( int4(3, -1, -2, 0) == geom::clamp(int4(3, -1, -2, 0), p, p2) );
-	assert( int4(-1, 3, -3, 1) == geom::clamp(int4(-3, 7, -5, 3), p, p2) );
-
-	(p3 = int4(-10, 2, 3, -3)).clamp(p,p2);
-	assert( int4(-1, 2, 2, -3) == p3 );
-	(p3 = int4(0, -3, -1, 10)).clamp(p,p2);
-	assert( int4(0, -2, -1, 1) == p3);
-	(p3 = int4(3, -1, -2, 0)).clamp(p,p2);
-	assert( int4(3, -1, -2, 0) == p3);
-	(p3 = int4(-3, 7, -5, 3)).clamp(p,p2);
-	assert( int4(-1, 3, -3, 1) == p3);
-
-	float4 fp (1.1f, 3.4f, 4.5f, 8.9f);
-	assert( float4(1.f, 3.f, 4.f, 8.f) == trunc(fp) );
-	assert( float4(1.f, 3.f, 4.f, 8.f) == floor(fp) );
-	assert( float4(2.f, 4.f, 5.f, 9.f) == ceil(fp) );
-	assert( float4(1.f, 3.f, 5.f, 9.f) == round(fp) );
-
-	float4 fp2;
-	(fp2 = fp).trunc();
-	assert( float4(1.f, 3.f, 4.f, 8.f) == fp2 );
-	(fp2 = fp).floor();
-	assert( float4(1.f, 3.f, 4.f, 8.f) == fp2 );
-	(fp2 = fp).ceil();
-	assert( float4(2.f, 4.f, 5.f, 9.f) == fp2 );
-	(fp2 = fp).round();
-	assert( float4(1.f, 3.f, 5.f, 9.f) == fp2 );
-
-	fp = float4(1.1f, -3.4f, 4.5f, -8.9f);
-	assert( float4(1.f, -3.f, 4.f, -8.f) == trunc(fp) );
-	assert( float4(1.f, -4.f, 4.f, -9.f) == floor(fp) );
-	assert( float4(2.f, -3.f, 5.f, -8.f) == ceil(fp) );
-	assert( float4(1.f, -3.f, 5.f, -9.f) == round(fp) );
-
-	(fp2 = fp).trunc();
-	assert( float4(1.f, -3.f, 4.f, -8.f) == fp2 );
-	(fp2 = fp).floor();
-	assert( float4(1.f, -4.f, 4.f, -9.f) == fp2 );
-	(fp2 = fp).ceil();
-	assert( float4(2.f, -3.f, 5.f, -8.f) == fp2 );
-	(fp2 = fp).round();
-	assert( float4(1.f, -3.f, 5.f, -9.f) == fp2 );
-
-	int4 p4 (2, 6, 3, 0);
-	assert( p4.magnitude() ==  49);
-	assert( p4.quadrance() ==  49);
-	assert( p4.length() == 7 );
-	assert( magnitude(p4) ==  49);
-	assert( quadrance(p4) ==  49);
-	assert( length(p4) == 7 );
-
-
-	assert( signum(int4(-12, 34, 0, 1)) == int4(-1, 1, 0, 1) );
-	assert( signum(int4(9, -13, -77, 0)) == int4(1, -1, -1, 0) );
-	assert( signum(float4(-.3f, .12f, .0f, +.0f)) == float4(-1.f, 1.f, 0.f, 0.f) );
-
-	(p3 = int4(-12, 34, 0, 1)).signum();
-	assert( int4(-1, 1, 0, 1) == p3 );
-	(p3 = int4(9, -13, -77, 0)).signum();
-	assert( int4(1, -1, -1, 0) == p3 );
-	(fp2 = float4(-.3f, .12f, .0f, +.0f)).signum();
-	assert( float4(-1.f, 1.f, 0.f, 0.f) == fp2 );
-}
-
 int PotentiallyAmbigous(geom::vector<int, 2>){ return 2; }
 int PotentiallyAmbigous(geom::vector<int, 3>){ return 3; }
 void DisabiguateOnConstructorParameterCount()
@@ -560,7 +479,6 @@ int main()
 	Arithmetic();
 	DiscreteArithmetic();
 	ComparisonOperators();
-	Algorithms();
 	DisabiguateOnConstructorParameterCount();
 	NumericLimits();
 	CoordinateOrder();