diff --git a/examples/noise_fractal_2d/README.md b/examples/noise_fractal_2d/README.md
new file mode 100644
index 000000000..ea94dd58f
--- /dev/null
+++ b/examples/noise_fractal_2d/README.md
@@ -0,0 +1,16 @@
+# Example - noise_fractal_2d 📘
+
+This example demonstrates the usage of the V Scientific Library for generating pink noise.
+
+## Instructions
+
+1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io).
+2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)!
+3. Navigate to this directory.
+4. Run the example using the following command:
+
+```sh
+v run main.v
+```
+
+Enjoy exploring the capabilities of the V Scientific Library! 😊
diff --git a/examples/noise_fractal_2d/main.v b/examples/noise_fractal_2d/main.v
new file mode 100644
index 000000000..a5ef30d25
--- /dev/null
+++ b/examples/noise_fractal_2d/main.v
@@ -0,0 +1,54 @@
+module main
+
+import rand
+import vsl.noise
+import vsl.plot
+
+fn main() {
+	// Creation of the noise generator.
+	// Other noise functions and dimensions count are available.
+	rand.seed([u32(3155200429), u32(3208395956)])
+	mut generator := noise.Generator.new()
+	generator.randomize()
+
+	mut x := []f64{}
+	mut y := []f64{}
+	mut z := []f64{}
+
+	// 5 layers of simplex noise
+	octaves := 5
+	persistence := 0.5
+
+	for xx in 0 .. 200 {
+		for yy in 0 .. 200 {
+			mut total := 0.0
+			mut frequency := 1.0
+			mut amplitude := 1.0
+			mut max_value := 0.0
+
+			for _ in 0 .. octaves {
+				total += generator.simplex_2d(0.03 * xx * frequency, 0.03 * yy * frequency) * amplitude
+				max_value += amplitude
+				amplitude *= persistence
+				frequency *= 2.0
+			}
+
+			// Normalize to [-1, 1]
+			total /= max_value
+			x << xx
+			y << yy
+			z << total
+		}
+	}
+
+	mut plt := plot.Plot.new()
+	plt.heatmap(
+		x: x
+		y: y
+		z: z
+	)
+	plt.layout(
+		title: 'Pink `fractal` noise 2d example'
+	)
+	plt.show()!
+}
diff --git a/examples/noise_simplex_2d/README.md b/examples/noise_simplex_2d/README.md
new file mode 100644
index 000000000..f21dc9268
--- /dev/null
+++ b/examples/noise_simplex_2d/README.md
@@ -0,0 +1,16 @@
+# Example - noise_simple_2d 📘
+
+This example demonstrates the usage of the V Scientific Library for generating simplex noise.
+
+## Instructions
+
+1. Ensure you have the V compiler installed. You can download it from [here](https://vlang.io).
+2. Ensure you have the VSL installed. You can do it following the [installation guide](https://github.com/vlang/vsl?tab=readme-ov-file#-installation)!
+3. Navigate to this directory.
+4. Run the example using the following command:
+
+```sh
+v run main.v
+```
+
+Enjoy exploring the capabilities of the V Scientific Library! 😊
diff --git a/examples/noise_simplex_2d/main.v b/examples/noise_simplex_2d/main.v
new file mode 100644
index 000000000..031441cb3
--- /dev/null
+++ b/examples/noise_simplex_2d/main.v
@@ -0,0 +1,38 @@
+module main
+
+import rand
+import vsl.noise
+import vsl.plot
+
+fn main() {
+	// Creation of the noise generator.
+	// Other noise functions and dimensions count are available.
+	rand.seed([u32(3155200429), u32(3208395956)])
+	mut generator := noise.Generator.new()
+	generator.randomize()
+
+	mut x := []f64{}
+	mut y := []f64{}
+	mut z := []f64{}
+
+	for xx in 0 .. 40 {
+		for yy in 0 .. 40 {
+			// We need to scale the coordinates to control the frequency of the noise.
+			val := generator.simplex_2d(0.03 * xx, 0.03 * yy)
+			x << xx
+			y << yy
+			z << val
+		}
+	}
+
+	mut plt := plot.Plot.new()
+	plt.heatmap(
+		x: x
+		y: y
+		z: z
+	)
+	plt.layout(
+		title: 'Simplex noise 2d example'
+	)
+	plt.show()!
+}
diff --git a/ml/README.md b/ml/README.md
new file mode 100644
index 000000000..e9c8d40a5
--- /dev/null
+++ b/ml/README.md
@@ -0,0 +1,80 @@
+# VSL Machine Learning (vsl.ml)
+
+VSL aims to provide a robust set of tools for scientific computing with an emphasis
+on performance and ease of use. In the `vsl.ml` module, some machine learning
+models are designed as observers of data, meaning they re-train automatically when
+data changes, while others do not require this functionality.
+
+## Key Features
+
+- **Observers of Data**: Some machine learning models in VSL act as observers,
+  re-training automatically when data changes.
+- **High Performance**: Leverages V’s performance optimizations and can integrate
+  with C and Fortran libraries like Open BLAS and LAPACK.
+- **Versatile Algorithms**: Supports a variety of machine learning algorithms and
+  models.
+
+## Usage
+
+### Loading Data
+
+The `Data` struct in `vsl.ml` is designed to hold data in matrix format for machine
+learning tasks. Here's a brief overview of how to use it:
+
+#### Creating a Data Object
+
+You can create a `Data` object using the following methods:
+
+- `Data.new`: Creates a new `Data` object with specified dimensions.
+- `Data.from_raw_x`: Creates a `Data` object from raw x values (without y values).
+- `Data.from_raw_xy`: Creates a `Data` object from raw x and y values combined in a single matrix.
+- `Data.from_raw_xy_sep`: Creates a `Data` object from separate x and y raw values.
+
+### Data Methods
+
+The `Data` struct has several key methods to manage and manipulate data:
+
+- `set(x, y)`: Sets the x matrix and y vector and notifies observers.
+- `set_y(y)`: Sets the y vector and notifies observers.
+- `set_x(x)`: Sets the x matrix and notifies observers.
+- `split(ratio)`: Splits the data into two parts based on the given ratio.
+- `clone()`: Returns a deep copy of the Data object without observers.
+- `clone_with_same_x()`: Returns a deep copy of the Data object but shares the same x reference.
+- `add_observer(obs)`: Adds an observer to the data object.
+- `notify_update()`: Notifies observers of data changes.
+
+### Stat Observer
+
+The `Stat` struct is an observer of `Data`, providing statistical analysis of the
+data it observes. It automatically updates its statistics when the underlying data
+changes.
+
+## Observer Models
+
+The following machine learning models in VSL are compatible with the `Observer`
+pattern. This means they can observe data changes and automatically update
+themselves.
+
+### K-Means Clustering
+
+K-Means Clustering is used for unsupervised learning to group data points into
+clusters. As an observer model, it re-trains automatically when the data changes,
+which is useful for dynamic datasets that require continuous updates.
+
+### K-Nearest Neighbors (KNN)
+
+K-Nearest Neighbors (KNN) is used for classification tasks where the target
+variable is categorical. As an observer model, it re-trains automatically when the
+data changes, which is beneficial for datasets that are frequently updated.
+
+## Non-Observer Models
+
+The following machine learning models in VSL do not require the observer pattern
+and are trained once on a dataset without continuous updates.
+
+### Linear Regression
+
+Linear Regression is used for predicting a continuous target variable based on one
+or more predictor variables. It is typically trained once on a dataset and used to
+make predictions without requiring continuous updates. Hence, it is not implemented
+as an observer model.
diff --git a/noise/noise.v b/noise/noise.v
new file mode 100644
index 000000000..f5c1aa994
--- /dev/null
+++ b/noise/noise.v
@@ -0,0 +1,20 @@
+module noise
+
+import rand
+
+// Generator is a struct holding the permutation table used in perlin and simplex noise
+pub struct Generator {
+mut:
+	perm []int = rand.shuffle_clone(permutations) or { panic(err) }
+}
+
+// new is a function that return a new Generator struct
+pub fn Generator.new() Generator {
+	return Generator{}
+}
+
+// randomize is a function that shuffle the permutation set inside the Generator struct
+// will not shuffle if rand.seed is not changed
+pub fn (mut generator Generator) randomize() {
+	generator.perm = rand.shuffle_clone(permutations) or { panic(err) }
+}
diff --git a/noise/perlin2d.v b/noise/perlin2d.v
index 4e5d2661e..0b60ffda0 100644
--- a/noise/perlin2d.v
+++ b/noise/perlin2d.v
@@ -1,26 +1,7 @@
 module noise
 
-import rand
-
-// Perlin is a struct that hold the permutation set for perlin noise
-pub struct Perlin {
-mut:
-	perm []int = rand.shuffle_clone(permutations) or { panic(err) }
-}
-
-// new is a function that return a new Perlin struct
-pub fn Perlin.new() Perlin {
-	return Perlin{}
-}
-
-// randomize is a function that shuffle the permutation set inside the Perlin struct
-// will not shuffle if rand.seed is not changed
-pub fn (mut perlin Perlin) randomize() {
-	perlin.perm = rand.shuffle_clone(permutations) or { panic(err) }
-}
-
 // perlin2d is a function that return a single value of perlin noise for a given 2d position
-pub fn (perlin Perlin) perlin2d(x f64, y f64) f64 {
+pub fn (generator Generator) perlin2d(x f64, y f64) f64 {
 	xi := int(x) & 0xFF
 	yi := int(y) & 0xFF
 
@@ -30,13 +11,13 @@ pub fn (perlin Perlin) perlin2d(x f64, y f64) f64 {
 	u := fade(xf)
 	v := fade(yf)
 
-	pxi := perlin.perm[xi]
-	pxi1 := perlin.perm[xi + 1]
+	pxi := generator.perm[xi]
+	pxi1 := generator.perm[xi + 1]
 
-	aa := perlin.perm[pxi + yi]
-	ab := perlin.perm[pxi + yi + 1]
-	ba := perlin.perm[pxi1 + yi]
-	bb := perlin.perm[pxi1 + yi + 1]
+	aa := generator.perm[pxi + yi]
+	ab := generator.perm[pxi + yi + 1]
+	ba := generator.perm[pxi1 + yi]
+	bb := generator.perm[pxi1 + yi + 1]
 
 	x1 := lerp(grad2d(aa, xf, yf), grad2d(ba, xf - 1, yf), u)
 	x2 := lerp(grad2d(ab, xf, yf - 1), grad2d(bb, xf - 1, yf - 1), u)
diff --git a/noise/perlin2d_test.v b/noise/perlin2d_test.v
index 6b4c81d83..f63a5d3e2 100644
--- a/noise/perlin2d_test.v
+++ b/noise/perlin2d_test.v
@@ -6,7 +6,7 @@ import vsl.float.float64
 fn test_perlin2d() {
 	rand.seed([u32(3155200429), u32(3208395956)])
 
-	mut gen := Perlin.new()
+	mut gen := Generator.new()
 	gen.randomize()
 
 	result := gen.perlin2d(0.125, 0.125)
diff --git a/noise/perlin3d.v b/noise/perlin3d.v
index c32a6131e..9b77e820b 100644
--- a/noise/perlin3d.v
+++ b/noise/perlin3d.v
@@ -1,7 +1,7 @@
 module noise
 
 // perlin3d is a function that return a single value of perlin noise for a given 3d position
-pub fn (perlin Perlin) perlin3d(x f64, y f64, z f64) f64 {
+pub fn (generator Generator) perlin3d(x f64, y f64, z f64) f64 {
 	xi := int(x) & 0xFF
 	yi := int(y) & 0xFF
 	zi := int(z) & 0xFF
@@ -12,21 +12,21 @@ pub fn (perlin Perlin) perlin3d(x f64, y f64, z f64) f64 {
 	v := fade(yf)
 	w := fade(zf)
 
-	pxi := perlin.perm[xi]
-	pxi_yi := perlin.perm[pxi + yi]
-	pxi_yi1 := perlin.perm[pxi + yi + 1]
-	pxi1 := perlin.perm[xi + 1]
-	pxi1_yi := perlin.perm[pxi1 + yi]
-	pxi1_yi1 := perlin.perm[pxi1 + yi + 1]
+	pxi := generator.perm[xi]
+	pxi_yi := generator.perm[pxi + yi]
+	pxi_yi1 := generator.perm[pxi + yi + 1]
+	pxi1 := generator.perm[xi + 1]
+	pxi1_yi := generator.perm[pxi1 + yi]
+	pxi1_yi1 := generator.perm[pxi1 + yi + 1]
 
-	aaa := perlin.perm[pxi_yi + zi]
-	aba := perlin.perm[pxi_yi1 + zi]
-	aab := perlin.perm[pxi_yi + zi + 1]
-	abb := perlin.perm[pxi_yi1 + zi + 1]
-	baa := perlin.perm[pxi1_yi + zi]
-	bba := perlin.perm[pxi1_yi1 + zi]
-	bab := perlin.perm[pxi1_yi + zi + 1]
-	bbb := perlin.perm[pxi1_yi1 + zi + 1]
+	aaa := generator.perm[pxi_yi + zi]
+	aba := generator.perm[pxi_yi1 + zi]
+	aab := generator.perm[pxi_yi + zi + 1]
+	abb := generator.perm[pxi_yi1 + zi + 1]
+	baa := generator.perm[pxi1_yi + zi]
+	bba := generator.perm[pxi1_yi1 + zi]
+	bab := generator.perm[pxi1_yi + zi + 1]
+	bbb := generator.perm[pxi1_yi1 + zi + 1]
 
 	mut x1 := lerp(grad3d(aaa, xf, yf, zf), grad3d(baa, xf - 1, yf, zf), u)
 	mut x2 := lerp(grad3d(aba, xf, yf - 1, zf), grad3d(bba, xf - 1, yf - 1, zf), u)
diff --git a/noise/perlin3d_test.v b/noise/perlin3d_test.v
index 8edde204c..0f5a67aed 100644
--- a/noise/perlin3d_test.v
+++ b/noise/perlin3d_test.v
@@ -6,7 +6,7 @@ import vsl.float.float64
 fn test_perlin3d() {
 	rand.seed([u32(3155200429), u32(3208395956)])
 
-	mut gen := Perlin.new()
+	mut gen := Generator.new()
 	gen.randomize()
 
 	result := gen.perlin3d(0.125, 0.125, 0.125)
diff --git a/noise/simplex.v b/noise/simplex.v
new file mode 100644
index 000000000..176fe40a6
--- /dev/null
+++ b/noise/simplex.v
@@ -0,0 +1,389 @@
+module noise
+
+import math
+
+// skewing factors for 2D case
+const f2 = 0.5 * (math.sqrt(3.0) - 1.0)
+const g2 = (3.0 - math.sqrt(3)) / 6.0
+// skewing factors for 3D case
+const f3 = f64(1.0 / 3.0)
+const g3 = f64(1.0 / 6.0)
+// skewing factors for 4D case
+const f4 = f64((math.sqrt(5.0) - 1.0) / 4)
+const g4 = f64((5.0 - math.sqrt(5.0)) / 20.0)
+
+// vfmt off
+const simplex := [
+  [0, 1, 2, 3], [0, 1, 3, 2], [0, 0, 0, 0], [0, 2, 3, 1], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [1, 2, 3, 0],
+  [0, 2, 1, 3], [0, 0, 0, 0], [0, 3, 1, 2], [0, 3, 2, 1], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [1, 3, 2, 0],
+  [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0],
+  [1, 2, 0, 3], [0, 0, 0, 0], [1, 3, 0, 2], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [2, 3, 0, 1], [2, 3, 1, 0],
+  [1, 0, 2, 3], [1, 0, 3, 2], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [2, 0, 3, 1], [0, 0, 0, 0], [2, 1, 3, 0],
+  [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0],
+  [2, 0, 1, 3], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [3, 0, 1, 2], [3, 0, 2, 1], [0, 0, 0, 0], [3, 1, 2, 0],
+  [2, 1, 0, 3], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [3, 1, 0, 2], [0, 0, 0, 0], [3, 2, 0, 1], [3, 2, 1, 0]
+]
+// vfmt on
+
+// grad_1d returns a gradient value for a given hash and x position
+fn grad_1d(hash int, x f64) f64 {
+	h := hash & 15
+	mut grad := 1.0 + f64(h & 7)
+	grad = if h & 8 == 0 { grad } else { -grad }
+	return grad * x
+}
+
+// grad_2d returns a gradient value for a given hash and x, y position
+fn grad_2d(hash int, x f64, y f64) f64 {
+	h := hash & 7
+	u := if h < 4 { x } else { y }
+	v := if h < 4 { y } else { x }
+	return (if h & 1 == 0 {
+		u
+	} else {
+		-u
+	}) + (if h & 2 == 0 {
+		2.0 * v
+	} else {
+		-2.0 * v
+	})
+}
+
+// grad_3d returns a gradient value for a given hash and x, y, z position
+fn grad_3d(hash int, x f64, y f64, z f64) f64 {
+	h := hash & 15
+	u := if h < 8 { x } else { y }
+	v := if h < 4 {
+		y
+	} else {
+		if h == 12 || h == 14 { x } else { z }
+	}
+	return (if h & 1 == 0 {
+		u
+	} else {
+		-u
+	}) + (if h & 2 == 0 {
+		v
+	} else {
+		-v
+	})
+}
+
+// grad_4d returns a gradient value for a given hash and x, y, z, t position
+fn grad_4d(hash int, x f64, y f64, z f64, t f64) f64 {
+	h := hash & 31
+	u := if h < 24 { x } else { y }
+	v := if h < 16 { y } else { z }
+	w := if h < 8 { z } else { t }
+	return (if h & 1 == 0 {
+		u
+	} else {
+		-u
+	}) + (if h & 2 == 0 {
+		v
+	} else {
+		-v
+	}) + (if h & 4 == 0 {
+		w
+	} else {
+		-w
+	})
+}
+
+// simplex_1d returns a simplex noise value for a given x position
+pub fn (generator Generator) simplex_1d(x f64) f64 {
+	i0 := int(x)
+	i1 := i0 + 1
+	x0 := x - i0
+	x1 := x0 - 1
+
+	mut t0 := 1.0 - x0 * x0
+	t0 *= t0
+	n0 := t0 * t0 * grad_1d(generator.perm[i0 & 0xff], x0)
+
+	mut t1 := 1.0 - x1 * x1
+	t1 *= t1
+	n1 := t1 * t1 * grad_1d(generator.perm[i1 & 0xff], x1)
+
+	return 0.395 * (n0 + n1)
+}
+
+// simplex_2d returns a simplex noise value for a given x, y position
+pub fn (generator Generator) simplex_2d(x f64, y f64) f64 {
+	s := (x + y) * noise.f2
+	i := int(x + s)
+	j := int(y + s)
+
+	t := f64(i + j) * noise.g2
+	x0 := x - (i - t)
+	y0 := y - (j - t)
+
+	i1, j1 := if x0 > y0 {
+		1, 0
+	} else {
+		0, 1
+	}
+
+	x1 := x0 - f64(i1) + noise.g2
+	y1 := y0 - f64(j1) + noise.g2
+	x2 := x0 - 1.0 + noise.g2 * 2.0
+	y2 := y0 - 1.0 + noise.g2 * 2.0
+
+	ii := i & 0xff
+	jj := j & 0xff
+
+	mut t0 := 0.5 - x0 * x0 - y0 * y0
+	mut n0 := 0.0
+	if t0 < 0 {
+		n0 = 0.0
+	} else {
+		t0 *= t0
+		n0 = t0 * t0 * grad_2d(generator.perm[ii + generator.perm[jj]], x0, y0)
+	}
+
+	mut t1 := 0.5 - x1 * x1 - y1 * y1
+	mut n1 := 0.0
+	if t1 < 0 {
+		n1 = 0.0
+	} else {
+		t1 *= t1
+		n1 = t1 * t1 * grad_2d(generator.perm[ii + i1 + generator.perm[jj + j1]], x1,
+			y1)
+	}
+
+	mut t2 := 0.5 - x2 * x2 - y2 * y2
+	mut n2 := 0.0
+	if t2 < 0 {
+		n2 = 0.0
+	} else {
+		t2 *= t2
+		n2 = t2 * t2 * grad_2d(generator.perm[ii + 1 + generator.perm[jj + 1]], x2, y2)
+	}
+
+	return 40.0 * (n0 + n1 + n2)
+}
+
+// simplex_3d returns a simplex noise value for a given x, y, z position
+pub fn (generator Generator) simplex_3d(x f64, y f64, z f64) f64 {
+	s := (x + y + z) * noise.f3
+	xs := x + s
+	ys := y + s
+	zs := z + s
+	i := int(xs)
+	j := int(ys)
+	k := int(zs)
+
+	t := f64(i + j + k) * noise.g3
+	x0 := x - (i - t)
+	y0 := y - (j - t)
+	z0 := z - (k - t)
+
+	mut i1 := 0
+	mut j1 := 0
+	mut k1 := 0
+	mut i2 := 0
+	mut j2 := 0
+	mut k2 := 0
+
+	if x0 >= y0 {
+		if y0 >= z0 {
+			i1 = 1
+			i2 = 1
+			j2 = 1
+		} else if x0 >= z0 {
+			i1 = 1
+			i2 = 1
+			k2 = 1
+		} else {
+			k1 = 1
+			i2 = 1
+			k2 = 1
+		}
+	} else {
+		if y0 < z0 {
+			k1 = 1
+			j2 = 1
+			k2 = 1
+		} else if x0 < z0 {
+			j1 = 1
+			j2 = 1
+			k2 = 1
+		} else {
+			j1 = 1
+			i2 = 1
+			j2 = 1
+		}
+	}
+
+	x1 := x0 - i1 + noise.g3
+	y1 := y0 - j1 + noise.g3
+	z1 := z0 - k1 + noise.g3
+	x2 := x0 - i2 + 2.0 * noise.g3
+	y2 := y0 - j2 + 2.0 * noise.g3
+	z2 := z0 - k2 + 2.0 * noise.g3
+	x3 := x0 - 1.0 + 3.0 * noise.g3
+	y3 := y0 - 1.0 + 3.0 * noise.g3
+	z3 := z0 - 1.0 + 3.0 * noise.g3
+
+	ii := i & 0xff
+	jj := j & 0xff
+	kk := k & 0xff
+
+	mut t0 := 0.6 - x0 * x0 - y0 * y0 - z0 * z0
+	mut n0 := 0.0
+	if t0 < 0 {
+		n0 = 0.0
+	} else {
+		t0 *= t0
+		n0 = t0 * t0 * grad_3d(generator.perm[ii + generator.perm[jj + generator.perm[kk]]],
+			x0, y0, z0)
+	}
+
+	mut t1 := 0.6 - x1 * x1 - y1 * y1 - z1 * z1
+	mut n1 := 0.0
+	if t1 < 0 {
+		n1 = 0.0
+	} else {
+		t1 *= t1
+		n1 = t1 * t1 * grad_3d(generator.perm[ii + i1 + generator.perm[jj + j1 + generator.perm[kk +
+			k1]]], x1, y1, z1)
+	}
+
+	mut t2 := 0.6 - x2 * x2 - y2 * y2 - z2 * z2
+	mut n2 := 0.0
+	if t2 < 0 {
+		n2 = 0.0
+	} else {
+		t2 *= t2
+		n2 = t2 * t2 * grad_3d(generator.perm[ii + i2 + generator.perm[jj + j2 + generator.perm[kk +
+			k2]]], x2, y2, z2)
+	}
+
+	mut t3 := 0.6 - x3 * x3 - y3 * y3 - z3 * z3
+	mut n3 := 0.0
+	if t3 < 0 {
+		n3 = 0.0
+	} else {
+		t3 *= t3
+		n3 = t3 * t3 * grad_3d(generator.perm[ii + 1 + generator.perm[jj + 1 + generator.perm[kk +
+			1]]], x3, y3, z3)
+	}
+
+	return 32.0 * (n0 + n1 + n2 + n3)
+}
+
+// simplex_4d returns a simplex noise value for a given x, y, z, w position
+pub fn (generator Generator) simplex_4d(x f64, y f64, z f64, w f64) f64 {
+	s := (x + y + z + w) * noise.f4
+	xs := x + s
+	ys := y + s
+	zs := z + s
+	ws := w + s
+	i := int(xs)
+	j := int(ys)
+	k := int(zs)
+	l := int(ws)
+
+	t := f64(i + j + k + l) * noise.g4
+	x0 := x - (i - t)
+	y0 := y - (j - t)
+	z0 := z - (k - t)
+	w0 := w - (l - t)
+
+	c1 := if x0 > y0 { 32 } else { 0 }
+	c2 := if x0 > z0 { 16 } else { 0 }
+	c3 := if y0 > z0 { 8 } else { 0 }
+	c4 := if x0 > w0 { 4 } else { 0 }
+	c5 := if y0 > w0 { 2 } else { 0 }
+	c6 := if z0 > w0 { 1 } else { 0 }
+	c := c1 + c2 + c3 + c4 + c5 + c6
+
+	i1 := if noise.simplex[c][0] >= 3 { 1 } else { 0 }
+	j1 := if noise.simplex[c][1] >= 3 { 1 } else { 0 }
+	k1 := if noise.simplex[c][2] >= 3 { 1 } else { 0 }
+	l1 := if noise.simplex[c][3] >= 3 { 1 } else { 0 }
+
+	i2 := if noise.simplex[c][0] >= 2 { 1 } else { 0 }
+	j2 := if noise.simplex[c][1] >= 2 { 1 } else { 0 }
+	k2 := if noise.simplex[c][2] >= 2 { 1 } else { 0 }
+	l2 := if noise.simplex[c][3] >= 2 { 1 } else { 0 }
+
+	i3 := if noise.simplex[c][0] >= 1 { 1 } else { 0 }
+	j3 := if noise.simplex[c][1] >= 1 { 1 } else { 0 }
+	k3 := if noise.simplex[c][2] >= 1 { 1 } else { 0 }
+	l3 := if noise.simplex[c][3] >= 1 { 1 } else { 0 }
+
+	x1 := x0 - i1 + noise.g4
+	y1 := y0 - j1 + noise.g4
+	z1 := z0 - k1 + noise.g4
+	w1 := w0 - l1 + noise.g4
+	x2 := x0 - i2 + 2.0 * noise.g4
+	y2 := y0 - j2 + 2.0 * noise.g4
+	z2 := z0 - k2 + 2.0 * noise.g4
+	w2 := w0 - l2 + 2.0 * noise.g4
+	x3 := x0 - i3 + 3.0 * noise.g4
+	y3 := y0 - j3 + 3.0 * noise.g4
+	z3 := z0 - k3 + 3.0 * noise.g4
+	w3 := w0 - l3 + 3.0 * noise.g4
+	x4 := x0 - 1.0 + 4.0 * noise.g4
+	y4 := y0 - 1.0 + 4.0 * noise.g4
+	z4 := z0 - 1.0 + 4.0 * noise.g4
+	w4 := w0 - 1.0 + 4.0 * noise.g4
+
+	ii := i & 0xff
+	jj := j & 0xff
+	kk := k & 0xff
+	ll := l & 0xff
+
+	mut t0 := 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0
+	mut n0 := 0.0
+	if t0 < 0.0 {
+		n0 = 0.0
+	} else {
+		t0 *= t0
+		n0 = t0 * t0 * grad_4d(generator.perm[ii + generator.perm[jj + generator.perm[kk +
+			generator.perm[ll]]]], x0, y0, z0, w0)
+	}
+
+	mut t1 := 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1
+	mut n1 := 0.0
+	if t1 < 0.0 {
+		n1 = 0.0
+	} else {
+		t1 *= t1
+		n1 = t1 * t1 * grad_4d(generator.perm[ii + i1 + generator.perm[jj + j1 + generator.perm[kk +
+			k1 + generator.perm[ll + l1]]]], x1, y1, z1, w1)
+	}
+
+	mut t2 := 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2
+	mut n2 := 0.0
+	if t2 < 0.0 {
+		n2 = 0.0
+	} else {
+		t2 *= t2
+		n2 = t2 * t2 * grad_4d(generator.perm[ii + i2 + generator.perm[jj + j2 + generator.perm[kk +
+			k2 + generator.perm[ll + l2]]]], x2, y2, z2, w2)
+	}
+
+	mut t3 := 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3
+	mut n3 := 0.0
+	if t3 < 0.0 {
+		n3 = 0.0
+	} else {
+		t3 *= t3
+		n3 = t3 * t3 * grad_4d(generator.perm[ii + i3 + generator.perm[jj + j3 + generator.perm[kk +
+			k3 + generator.perm[ll + l3]]]], x3, y3, z3, w3)
+	}
+
+	mut t4 := 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4
+	mut n4 := 0.0
+	if t4 < 0.0 {
+		n4 = 0.0
+	} else {
+		t4 *= t4
+		n4 = t4 * t4 * grad_4d(generator.perm[ii + 1 + generator.perm[jj + 1 + generator.perm[kk +
+			1 + generator.perm[ll + 1]]]], x4, y4, z4, w4)
+	}
+
+	return 27.0 * (n0 + n1 + n2 + n3 + n4)
+}
diff --git a/noise/simplex_test.v b/noise/simplex_test.v
new file mode 100644
index 000000000..66dd946a7
--- /dev/null
+++ b/noise/simplex_test.v
@@ -0,0 +1,40 @@
+module noise
+
+import rand
+import vsl.float.float64
+
+fn test_simplex_1d() {
+	rand.seed([u32(3155200429), u32(3208395956)])
+	mut generator := Generator.new()
+	generator.randomize()
+	result := generator.simplex_1d(0.287)
+	expected := -0.3544283326507284
+	assert float64.tolerance(result, expected, 1.0e-6)
+}
+
+fn test_simplex_2d() {
+	rand.seed([u32(3075200429), u32(3094395956)])
+	mut generator := Generator.new()
+	generator.randomize()
+	result := generator.simplex_2d(0.287, 0.475)
+	expected := -0.09948661872545192
+	assert float64.tolerance(result, expected, 1.0e-6)
+}
+
+fn test_simplex_3d() {
+	rand.seed([u32(3155200429), u32(3208395956)])
+	mut generator := Generator.new()
+	generator.randomize()
+	result := generator.simplex_3d(0.287, 0.475, 1.917)
+	expected := -0.06034653476116279
+	assert float64.tolerance(result, expected, 1.0e-6)
+}
+
+fn test_simplex_4d() {
+	rand.seed([u32(3075200429), u32(3094395956)])
+	mut generator := Generator.new()
+	generator.randomize()
+	result := generator.simplex_4d(0.287, 0.475, 1.917, 0.684)
+	expected := 0.015098415881100141
+	assert float64.tolerance(result, expected, 1.0e-6)
+}
diff --git a/poly/README.md b/poly/README.md
index 2fe043dcb..0e6855662 100644
--- a/poly/README.md
+++ b/poly/README.md
@@ -1,13 +1,16 @@
 # Polynomials
 
-This chapter describes functions for evaluating and solving polynomials.
-There are routines for finding real and complex roots of quadratic and
-cubic equations using analytic methods. An iterative polynomial solver
-is also available for finding the roots of general polynomials with real
-coefficients (of any order). The functions are declared in the module `vsl.poly`.
+This chapter describes functions for evaluating and solving polynomials. There are routines
+for finding real and complex roots of quadratic and cubic equations using analytic methods.
+An iterative polynomial solver is also available for finding the roots of general polynomials
+with real coefficients (of any order). The functions are declared in the module `vsl.poly`.
 
 ## Polynomial Evaluation
 
+```v ignore
+fn eval(c []f64, x f64) f64
+```
+
 The functions described here evaluate the polynomial
 
 ```console
@@ -16,20 +19,13 @@ P(x) = c[0] + c[1] x + c[2] x^2 + . . . + c[len-1] x^(len-1)
 
 using Horner's method for stability.
 
-```v ignore
-fn eval(c []f64, x f64) f64
-```
-
-This function evaluates a polynomial with real coefficients for the real variable `x`.
-
 ```v ignore
 fn eval_derivs(c []f64, x f64, lenres u64) []f64
 ```
 
-This function evaluates a polynomial and its derivatives storing the
-results in the array `res` of size `lenres`. The output array
-contains the values of `d^k P(x)/d x^k` for the specified value of
-`x` starting with `k = 0`.
+This function evaluates a polynomial and its derivatives, storing the results in the array
+`res` of size `lenres`. The output array contains the values of `d^k P(x)/d x^k` for the
+specified value of `x`, starting with `k = 0`.
 
 ## Quadratic Equations
 
@@ -43,22 +39,17 @@ This function finds the real roots of the quadratic equation,
 a x^2 + b x + c = 0
 ```
 
-The number of real roots (either zero, one or two) is returned, and
-their locations are are returned as `[ x0, x1 ]`. If no real roots
-are found then `[]` is returned. If one real root
-is found (i.e. if `a=0`) then it is are returned as `[ x0 ]`. When two
-real roots are found they are are returned as `[ x0, x1 ]` in
-ascending order. The case of coincident roots is not considered
-special. For example `(x-1)^2=0` will have two roots, which happen
-to have exactly equal values.
+The number of real roots (either zero, one or two) is returned, and their locations are
+returned as `[ x0, x1 ]`. If no real roots are found then `[]` is returned. If one real root
+is found (i.e. if `a=0`) then it is returned as `[ x0 ]`. When two real roots are found they
+are returned as `[ x0, x1 ]` in ascending order. The case of coincident roots is not considered
+special. For example `(x-1)^2=0` will have two roots, which happen to have exactly equal values.
 
-The number of roots found depends on the sign of the discriminant
-`b^2 - 4 a c`. This will be subject to rounding and cancellation
-errors when computed in double precision, and will also be subject to
-errors if the coefficients of the polynomial are inexact. These errors
-may cause a discrete change in the number of roots. However, for
-polynomials with small integer coefficients the discriminant can always
-be computed exactly.
+The number of roots found depends on the sign of the discriminant `b^2 - 4 a c`. This will
+be subject to rounding and cancellation errors when computed in double precision, and will
+also be subject to errors if the coefficients of the polynomial are inexact. These errors may
+cause a discrete change in the number of roots. However, for polynomials with small integer
+coefficients the discriminant can always be computed exactly.
 
 ## Cubic Equations
 
@@ -72,13 +63,93 @@ This function finds the real roots of the cubic equation,
 x^3 + a x^2 + b x + c = 0
 ```
 
-with a leading coefficient of unity. The number of real roots (either
-one or three) is returned, and their locations are returned as `[ x0, x1, x2 ]`.
-If one real root is found then only `[ x0 ]`
-is returned. When three real roots are found they are returned as
-`[ x0, x1, x2 ]` in ascending order. The case of
-coincident roots is not considered special. For example, the equation
-`(x-1)^3=0` will have three roots with exactly equal values. As
-in the quadratic case, finite precision may cause equal or
-closely-spaced real roots to move off the real axis into the complex
-plane, leading to a discrete change in the number of real roots.
+with a leading coefficient of unity. The number of real roots (either one or three) is
+returned, and their locations are returned as `[ x0, x1, x2 ]`. If one real root is found
+then only `[ x0 ]` is returned. When three real roots are found they are returned as
+`[ x0, x1, x2 ]` in ascending order. The case of coincident roots is not considered special.
+For example, the equation `(x-1)^3=0` will have three roots with exactly equal values. As
+in the quadratic case, finite precision may cause equal or closely-spaced real roots to move
+off the real axis into the complex plane, leading to a discrete change in the number of real roots.
+
+## Companion Matrix
+
+```v ignore
+fn companion_matrix(a []f64) [][]f64
+```
+
+Creates a companion matrix for the polynomial
+
+```console
+P(x) = a_n * x^n + a_{n-1} * x^{n-1} + ... + a_1 * x + a_0
+```
+
+The companion matrix `C` is defined as:
+
+```
+[0 0 0 ... 0 -a_0/a_n]
+[1 0 0 ... 0 -a_1/a_n]
+[0 1 0 ... 0 -a_2/a_n]
+[. . . ... . ........]
+[0 0 0 ... 1 -a_{n-1}/a_n]
+```
+
+## Balanced Companion Matrix
+
+```v ignore
+fn balance_companion_matrix(cm [][]f64) [][]f64
+```
+
+Balances a companion matrix `C` to improve numerical stability. It uses an iterative scaling
+process to make the row and column norms as close to each other as possible. The output is
+a balanced matrix `B` such that `D^(-1)CD = B`, where `D` is a diagonal matrix.
+
+## Polynomial Operations
+
+```v ignore
+fn add(a []f64, b []f64) []f64
+```
+
+Adds two polynomials:
+
+```console
+(a_n * x^n + ... + a_0) + (b_m * x^m + ... + b_0)
+```
+
+Returns the result as `[a_0 + b_0, a_1 + b_1, ..., a_k + b_k ...]`.
+
+```v ignore
+fn subtract(a []f64, b []f64) []f64
+```
+
+Subtracts two polynomials:
+
+```console
+(a_n * x^n + ... + a_0) - (b_m * x^m + ... + b_0)
+```
+
+Returns the result as `[a_0 - b_0, a_1 - b_1, ..., a_k - b_k, ...]`.
+
+```v ignore
+fn multiply(a []f64, b []f64) []f64
+```
+
+Multiplies two polynomials:
+
+```console
+(a_n * x^n + ... + a_0) * (b_m * x^m + ... + b_0)
+```
+
+Returns the result as `[c_0, c_1, ..., c_{n+m}]` where `c_k = ∑_{i+j=k} a_i * b_j`.
+
+```v ignore
+fn divide(a []f64, b []f64) ([]f64, []f64)
+```
+
+Divides two polynomials:
+
+```console
+(a_n * x^n + ... + a_0) / (b_m * x^m + ... + b_0)
+```
+
+Uses polynomial long division algorithm. Returns `(q, r)` where `q` is the quotient and `r`
+is the remainder such that `a(x) = b(x) * q(x) + r(x)` and `degree(r) < degree(b)`.
diff --git a/poly/poly.v b/poly/poly.v
index 41926d9ab..401ce60e6 100644
--- a/poly/poly.v
+++ b/poly/poly.v
@@ -6,6 +6,10 @@ import vsl.errors
 const radix = 2
 const radix2 = (radix * radix)
 
+// eval is a function that evaluates a polynomial P(x) = a_n * x^n + a_{n-1} * x^{n-1} + ... + a_1 * x + a_0
+// using Horner's method: P(x) = (...((a_n * x + a_{n-1}) * x + a_{n-2}) * x + ... + a_1) * x + a_0
+// Input: c = [a_0, a_1, ..., a_n], x
+// Output: P(x)
 pub fn eval(c []f64, x f64) f64 {
 	if c.len == 0 {
 		errors.vsl_panic('coeficients can not be empty', .efailed)
@@ -18,6 +22,9 @@ pub fn eval(c []f64, x f64) f64 {
 	return ans
 }
 
+// eval_derivs evaluates a polynomial P(x) and its derivatives P'(x), P''(x), ..., P^(k)(x)
+// Input: c = [a_0, a_1, ..., a_n] representing P(x), x, and lenres (k+1)
+// Output: [P(x), P'(x), P''(x), ..., P^(k)(x)]
 pub fn eval_derivs(c []f64, x f64, lenres int) []f64 {
 	mut res := []f64{}
 	lenc := c.len
@@ -49,7 +56,11 @@ pub fn eval_derivs(c []f64, x f64, lenres int) []f64 {
 	return res
 }
 
-pub fn solve_quadratic(a f64, b f64, c f64) []f64 { // Handle linear case
+// solve_quadratic solves the quadratic equation ax^2 + bx + c = 0
+// using the quadratic formula: x = (-b ± √(b^2 - 4ac)) / (2a)
+// Input: a, b, c
+// Output: Array of real roots (if any)
+pub fn solve_quadratic(a f64, b f64, c f64) []f64 {
 	if a == 0 {
 		if b == 0 {
 			return []
@@ -76,6 +87,10 @@ pub fn solve_quadratic(a f64, b f64, c f64) []f64 { // Handle linear case
 	}
 }
 
+// solve_cubic solves the cubic equation x^3 + ax^2 + bx + c = 0
+// using Cardano's formula and trigonometric solution
+// Input: a, b, c
+// Output: Array of real roots
 pub fn solve_cubic(a f64, b f64, c f64) []f64 {
 	q_ := (a * a - 3.0 * b)
 	r_ := (2.0 * a * a * a - 9.0 * a * b + 27.0 * c)
@@ -88,15 +103,6 @@ pub fn solve_cubic(a f64, b f64, c f64) []f64 {
 	if r == 0.0 && q == 0.0 {
 		return [-a / 3.0, -a / 3.0, -a / 3.0]
 	} else if cr2 == cq3 {
-		/*
-		this test is actually r2 == q3, written in a form suitable
-                        for exact computation with integers
-		*/
-		/*
-		Due to finite precision some double roots may be missed, and
-                        considered to be a pair of complex roots z = x +/- epsilon i
-                        close to the real axis.
-		*/
 		sqrt_q := math.sqrt(q)
 		if r > 0.0 {
 			return [-2.0 * sqrt_q - a / 3.0, sqrt_q - a / 3.0, sqrt_q - a / 3.0]
@@ -121,11 +127,17 @@ pub fn solve_cubic(a f64, b f64, c f64) []f64 {
 	}
 }
 
+// swap_ swaps two numbers: f(a, b) = (b, a)
+// Input: a, b
+// Output: (b, a)
 @[inline]
 fn swap_(a f64, b f64) (f64, f64) {
 	return b, a
 }
 
+// sorted_3_ sorts three numbers in ascending order: f(x, y, z) = (min(x,y,z), median(x,y,z), max(x,y,z))
+// Input: x, y, z
+// Output: (min, median, max)
 @[inline]
 fn sorted_3_(x_ f64, y_ f64, z_ f64) (f64, f64, f64) {
 	mut x := x_
@@ -143,6 +155,15 @@ fn sorted_3_(x_ f64, y_ f64, z_ f64) (f64, f64, f64) {
 	return x, y, z
 }
 
+// companion_matrix creates a companion matrix for the polynomial P(x) = a_n * x^n + a_{n-1} * x^{n-1} + ... + a_1 * x + a_0
+// The companion matrix C is defined as:
+// [0 0 0 ... 0 -a_0/a_n]
+// [1 0 0 ... 0 -a_1/a_n]
+// [0 1 0 ... 0 -a_2/a_n]
+// [. . . ... . ........]
+// [0 0 0 ... 1 -a_{n-1}/a_n]
+// Input: a = [a_0, a_1, ..., a_n]
+// Output: Companion matrix C
 pub fn companion_matrix(a []f64) [][]f64 {
 	nc := a.len - 1
 	mut cm := [][]f64{len: nc, init: []f64{len: nc}}
@@ -161,6 +182,10 @@ pub fn companion_matrix(a []f64) [][]f64 {
 	return cm
 }
 
+// balance_companion_matrix balances a companion matrix C to improve numerical stability
+// Uses an iterative scaling process to make the row and column norms as close to each other as possible
+// Input: Companion matrix C
+// Output: Balanced matrix B such that D^(-1)CD = B, where D is a diagonal matrix
 pub fn balance_companion_matrix(cm [][]f64) [][]f64 {
 	nc := cm.len
 	mut m := cm.clone()
@@ -169,7 +194,7 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 {
 	mut col_norm := 0.0
 	for not_converged {
 		not_converged = false
-		for i := 0; i < nc; i++ { // column norm, excluding the diagonal
+		for i := 0; i < nc; i++ {
 			if i != nc - 1 {
 				col_norm = math.abs(m[i + 1][i])
 			} else {
@@ -177,7 +202,7 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 {
 				for j := 0; j < nc - 1; j++ {
 					col_norm += math.abs(m[j][nc - 1])
 				}
-			} // row norm, excluding the diagonal
+			}
 			if i == 0 {
 				row_norm = math.abs(m[0][nc - 1])
 			} else if i == nc - 1 {
@@ -222,86 +247,66 @@ pub fn balance_companion_matrix(cm [][]f64) [][]f64 {
 	return m
 }
 
-// Arithmetic operations on polynomials
-//
-// In the following descriptions a, b, c are polynomials of degree
-// na, nb, nc respectively.
-//
-// Each polynomial is represented by an array containing its
-// coefficients, together with a separately declared integer equal
-// to the degree of the polynomial.  The coefficients appear in
-// ascending order; that is,
-//
-// a(x)  =  a[0]  +  a[1] * x  +  a[2] * x^2   +  ...  +  a[na] * x^na  .
-//
-//
-//
-// sum = eval( a, x )    Evaluate polynomial a(t) at t = x.
-// c = add( a, b )   c = a + b, nc = max(na, nb)
-// c = sub( a, b )   c = a - b, nc = max(na, nb)
-// c = mul( a, b )   c = a * b, nc = na+nb
-//
-//
-// Division:
-//
-// c = div( a, b )       c = a / b, nc = MAXPOL
-//
-// returns i = the degree of the first nonzero coefficient of a.
-// The computed quotient c must be divided by x^i.  An error message
-// is printed if a is identically zero.
-//
-//
-// Change of variables:
-// If a and b are polynomials, and t = a(x), then
-//     c(t) = b(a(x))
-// is a polynomial found by substituting a(x) for t.  The
-// subroutine call for this is
-//
+// add adds two polynomials: (a_n * x^n + ... + a_0) + (b_m * x^m + ... + b_0)
+// Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m]
+// Output: [a_0 + b_0, a_1 + b_1, ..., max(a_k, b_k), ...]
 pub fn add(a []f64, b []f64) []f64 {
-	na := a.len
-	nb := b.len
-	nc := int(math.max(na, nb))
-	mut c := []f64{len: nc}
-	for i := 0; i < nc; i++ {
-		if i > na {
-			c[i] = b[i]
-		} else if i > nb {
-			c[i] = a[i]
-		} else {
-			c[i] = a[i] + b[i]
-		}
+	mut result := []f64{len: math.max(a.len, b.len)}
+	for i in 0 .. result.len {
+		result[i] = if i < a.len { a[i] } else { 0.0 } + if i < b.len { b[i] } else { 0.0 }
 	}
-	return c
+	return result
 }
 
-pub fn substract(a []f64, b []f64) []f64 {
-	na := a.len
-	nb := b.len
-	nc := int(math.max(na, nb))
-	mut c := []f64{len: nc}
-	for i := 0; i < nc; i++ {
-		if i > na {
-			c[i] = -b[i]
-		} else if i > nb {
-			c[i] = a[i]
-		} else {
-			c[i] = a[i] - b[i]
-		}
+// subtract subtracts two polynomials: (a_n * x^n + ... + a_0) - (b_m * x^m + ... + b_0)
+// Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m]
+// Output: [a_0 - b_0, a_1 - b_1, ..., a_k - b_k, ...]
+pub fn subtract(a []f64, b []f64) []f64 {
+	mut result := []f64{len: math.max(a.len, b.len)}
+	for i in 0 .. result.len {
+		result[i] = if i < a.len { a[i] } else { 0.0 } - if i < b.len { b[i] } else { 0.0 }
 	}
-	return c
+	return result
 }
 
+// multiply multiplies two polynomials: (a_n * x^n + ... + a_0) * (b_m * x^m + ... + b_0)
+// Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m]
+// Output: [c_0, c_1, ..., c_{n+m}] where c_k = ∑_{i+j=k} a_i * b_j
 pub fn multiply(a []f64, b []f64) []f64 {
-	na := a.len
-	nb := b.len
-	nc := na + nb
-	mut c := []f64{len: nc}
-	for i := 0; i < na; i++ {
-		x := a[i]
-		for j := 0; j < nb; j++ {
-			k := i + j
-			c[k] += x * b[j]
+	mut result := []f64{len: a.len + b.len - 1}
+	for i in 0 .. a.len {
+		for j in 0 .. b.len {
+			result[i + j] += a[i] * b[j]
+		}
+	}
+	return result
+}
+
+// divide divides two polynomials: (a_n * x^n + ... + a_0) / (b_m * x^m + ... + b_0)
+// Uses polynomial long division algorithm
+// Input: a = [a_0, ..., a_n], b = [b_0, ..., b_m]
+// Output: (q, r) where q is the quotient and r is the remainder
+// such that a(x) = b(x) * q(x) + r(x) and degree(r) < degree(b)
+pub fn divide(a []f64, b []f64) ([]f64, []f64) {
+	mut quotient := []f64{}
+	mut remainder := a.clone()
+	b_lead_coef := b[0]
+
+	for remainder.len >= b.len {
+		lead_coef := remainder[0] / b_lead_coef
+		quotient << lead_coef
+		for i in 0 .. b.len {
+			remainder[i] -= lead_coef * b[i]
+		}
+		remainder = remainder[1..]
+		for remainder.len > 0 && math.abs(remainder[0]) < 1e-10 {
+			remainder = remainder[1..]
 		}
 	}
-	return c
+
+	if remainder.len == 0 {
+		remainder = []f64{}
+	}
+
+	return quotient, remainder
 }
diff --git a/poly/poly_test.v b/poly/poly_test.v
index 8db4f0e1d..b9c25cdb2 100644
--- a/poly/poly_test.v
+++ b/poly/poly_test.v
@@ -1,5 +1,7 @@
 module poly
 
+import math
+
 fn test_eval() {
 	// ans = 2
 	// ans = 4.0 + 4 * 2 = 12
@@ -16,29 +18,38 @@ fn test_swap() {
 	assert a == 202.0 && b == 101.0
 }
 
-// fn test_sorted_3_(){
-// 	a := 5.0
-// 	b := 7.0
-// 	c := -8.0
-// 	x, y, z := sorted_3_(a, b, c)
-// 	assert y == 7.0
-// }
-
 fn test_add() {
-	a := [6.0, 777, -3]
-	b := [1.0, -755, -4]
-	assert add(a, b) == [7.0, 22, -7]
+	a := [1.0, 2.0, 3.0]
+	b := [6.0, 20.0, -10.0]
+	result := add(a, b)
+	println('Add result: ${result}')
+	assert result == [7.0, 22.0, -7.0]
 }
 
-fn test_substract() {
-	a := [6.0, 777, -3]
-	b := [1.0, -755, -4]
-	assert substract(a, b) == [5.0, 1532, 1]
+fn test_subtract() {
+	a := [6.0, 1532.0, -4.0]
+	b := [1.0, -1.0, -5.0]
+	result := subtract(a, b)
+	println('Subtract result: ${result}')
+	assert result == [5.0, 1533.0, 1.0]
 }
 
-// fn test_multiply(){
-// 	a := [9.0, -1, 5]
-// 	b := [0.0, -1, 7]
+fn test_multiply() {
+	// (2+3x+4x^2) * (-3x+2x^2) = (-6x -5x^2 -6x^3 + 8x^4)
+	a := [2.0, 3.0, 4.0]
+	b := [0.0, -3.0, 2.0]
+	result := multiply(a, b)
+	println('Multiply result: ${result}')
+	assert result == [0.0, -6.0, -5.0, -6.0, 8.0]
+}
 
-// 	assert multiply(a, b) == [0.0, -9, 73, -12, 35, 0]
-// }
+fn test_divide() {
+	// (x^2 + 2x + 1) / (x + 1) = (x + 1)
+	a := [1.0, 2.0, 1.0]
+	b := [1.0, 1.0]
+	quotient, remainder := divide(a, b)
+	println('Divide quotient: ${quotient}')
+	println('Divide remainder: ${remainder}')
+	assert quotient == [1.0, 1.0]
+	assert remainder == [] // The empty set indicates that two polynomials divide each other exactly (without remainder).
+}