diff --git a/documentation/fast_fourier_transform.md b/documentation/fast_fourier_transform.md
new file mode 100644
index 00000000..c51daa92
--- /dev/null
+++ b/documentation/fast_fourier_transform.md
@@ -0,0 +1,29 @@
+# Fast Fourier Transform (FFT)
+
+This file documents the recursive Cooley-Tukey FFT implementation added to `R/mathematics/fast_fourier_transform.r`.
+
+## Description
+
+The `fft_recursive` function computes the discrete Fourier transform (DFT) of a numeric or complex vector using a divide-and-conquer Cooley-Tukey algorithm. If the input length is not a power of two, it is zero-padded to the next power of two.
+
+## Usage
+
+In an R session:
+
+source('mathematics/fast_fourier_transform.r')
+fft_recursive(c(0, 1, 2, 3))
+
+From the command line with Rscript:
+
+Rscript -e "source('R/mathematics/fast_fourier_transform.r'); print(fft_recursive(c(0,1,2,3)))"
+
+## Complexity
+
+Time complexity: O(n log n) for inputs with length a power of two; otherwise dominated by padding to next power of two.
+
+Space complexity: O(n) additional space for recursive calls.
+
+## Notes
+
+- The function returns a complex vector of the same length (after padding) as the input.
+- This implementation is primarily educational; production code should prefer the optimized `fft` function available in base R.
diff --git a/mathematics/fast_fourier_transform.r b/mathematics/fast_fourier_transform.r
new file mode 100644
index 00000000..0ce03f96
--- /dev/null
+++ b/mathematics/fast_fourier_transform.r
@@ -0,0 +1,53 @@
+# Fast Fourier Transform (Cooley-Tukey recursive implementation)
+#
+# This implementation accepts a numeric or complex vector and returns
+# its discrete Fourier transform as a complex vector. If the input length
+# is not a power of two, the vector is zero-padded to the next power of two.
+#
+# Usage:
+#   source('mathematics/fast_fourier_transform.r')
+#   x <- c(0,1,2,3)
+#   fft_result <- fft_recursive(x)
+#   print(fft_result)
+
+next_power_of_two <- function(n) {
+  if (n <= 0) return(1)
+  p <- 1
+  while (p < n) p <- p * 2
+  p
+}
+
+fft_recursive <- function(x) {
+  # Ensure input is complex
+  x <- as.complex(x)
+  N <- length(x)
+
+  # Pad to next power of two if necessary
+  M <- next_power_of_two(N)
+  if (M != N) {
+    x <- c(x, rep(0+0i, M - N))
+    N <- M
+  }
+
+  if (N == 1) return(x)
+
+  even <- fft_recursive(x[seq(1, N, by = 2)])
+  odd  <- fft_recursive(x[seq(2, N, by = 2)])
+
+  factor <- exp(-2i * pi * (0:(N/2 - 1)) / N)
+  T <- factor * odd
+
+  c(even + T, even - T)
+}
+
+# Example usage when run directly with Rscript
+if (identical(Sys.getenv("R_SCRIPT_NAME"), "") && interactive()) {
+  # Running in interactive R session - show sample
+  x <- c(0, 1, 2, 3)
+  cat("Input:\n")
+  print(x)
+  cat("FFT result:\n")
+  print(fft_recursive(x))
+}
+
+# When running via Rscript, users can call: Rscript -e "source('R/mathematics/fast_fourier_transform.r'); print(fft_recursive(c(0,1,2,3)))"