diff --git a/README.md b/README.md
index ac904a3..3836a67 100644
--- a/README.md
+++ b/README.md
@@ -108,6 +108,8 @@ To facilitate the dual-license model, contributions to the project can only be a
     *  I2C slave AXI lite master
     *  MDIO master
     *  UART
+*  Math
+    *  MT19937/MT19937-64 Mersenne Twister PRNG
 *  PCI Express
     *  PCIe AXI lite master
     *  PCIe AXI lite master for Xilinx UltraScale
diff --git a/src/math/lib/taxi b/src/math/lib/taxi
new file mode 120000
index 0000000..1b20c9f
--- /dev/null
+++ b/src/math/lib/taxi
@@ -0,0 +1 @@
+../../../
\ No newline at end of file
diff --git a/src/math/rtl/taxi_mt19937.sv b/src/math/rtl/taxi_mt19937.sv
new file mode 100644
index 0000000..8dbbe3d
--- /dev/null
+++ b/src/math/rtl/taxi_mt19937.sv
@@ -0,0 +1,269 @@
+// SPDX-License-Identifier: CERN-OHL-S-2.0
+/*
+
+Copyright (c) 2014-2025 FPGA Ninja, LLC
+
+Authors:
+- Alex Forencich
+
+*/
+
+`resetall
+`timescale 1ns / 1ps
+`default_nettype none
+
+/*
+ * MT19937/MT19937-64 Mersenne Twister PRNG
+ */
+module taxi_mt19937 #(
+    parameter integer MT_W = 32,
+    parameter logic [MT_W-1:0] INIT_SEED = 5489
+)
+(
+    input  wire logic             clk,
+    input  wire logic             rst,
+
+    /*
+     * AXI output
+     */
+    taxi_axis_if.src              m_axis,
+
+    /*
+     * Status
+     */
+    output wire logic             busy,
+
+    /*
+     * Configuration
+     */
+    input  wire logic [MT_W-1:0]  seed_val,
+    input  wire logic             seed_start
+);
+
+// check configuration
+if (MT_W != 32 && MT_W != 64)
+    $fatal(0, "Error: MT_W must be 32 or 64 (instance %m)", MT_W);
+
+if (m_axis.DATA_W != MT_W)
+    $fatal(0, "Error: Interface DATA_W parameter must be %d (instance %m)", MT_W);
+
+localparam MT_N = MT_W == 64 ? 312 : 624;
+localparam MT_M = MT_W == 64 ? 156 : 397;
+
+localparam CL_MT_N = $clog2(MT_N);
+
+// state register
+localparam [0:0]
+    STATE_SEED = 1'd0,
+    STATE_RUN = 1'd1;
+
+logic [0:0] state_reg = STATE_SEED, state_next;
+
+logic [MT_W-1:0] mt_ram[MT_N];
+logic [MT_W-1:0] mt_save_reg = '0, mt_save_next;
+logic [CL_MT_N-1:0] mti_reg = '0, mti_next;
+
+logic [CL_MT_N-1:0] mt_wr_ptr;
+logic [MT_W-1:0] mt_wr_data;
+logic mt_wr_en;
+
+logic [CL_MT_N-1:0] mt_rd_a_ptr_reg = '0, mt_rd_a_ptr_next;
+logic [MT_W-1:0] mt_rd_a_data_reg = '0;
+
+logic [CL_MT_N-1:0] mt_rd_b_ptr_reg = '0, mt_rd_b_ptr_next;
+logic [MT_W-1:0] mt_rd_b_data_reg = '0;
+
+logic [MT_W-1:0] product_reg = INIT_SEED, product_next;
+logic [MT_W-1:0] factor1_reg = '0, factor1_next;
+logic [MT_W-1:0] factor2_reg = '0, factor2_next;
+logic mul_done_reg = 1'b1, mul_done_next;
+
+logic [MT_W-1:0] m_axis_tdata_reg = '0, m_axis_tdata_next;
+logic m_axis_tvalid_reg = 1'b0, m_axis_tvalid_next;
+
+logic busy_reg = 1'b0;
+
+assign m_axis.tdata = m_axis_tdata_reg;
+assign m_axis.tkeep  = '1;
+assign m_axis.tstrb  = m_axis.tkeep;
+assign m_axis.tvalid = m_axis_tvalid_reg;
+assign m_axis.tlast  = 1'b1;
+assign m_axis.tid    = '0;
+assign m_axis.tdest  = '0;
+assign m_axis.tuser  = '0;
+
+assign busy = busy_reg;
+
+wire [MT_W-1:0] y1, y2, y3, y4, y5;
+
+if (MT_W == 32) begin
+
+    assign y1 = {mt_save_reg[31], mt_rd_a_data_reg[30:0]};
+    assign y2 = mt_rd_b_data_reg ^ (y1 >> 1) ^ (y1[0] ? 32'h9908b0df : 32'h0);
+    assign y3 = y2 ^ (y2 >> 11);
+    assign y4 = y3 ^ ((y3 << 7) & 32'h9d2c5680);
+    assign y5 = y4 ^ ((y4 << 15) & 32'hefc60000);
+
+end else begin
+
+    assign y1 = {mt_save_reg[63:31], mt_rd_a_data_reg[30:0]};
+    assign y2 = mt_rd_b_data_reg ^ (y1 >> 1) ^ (y1[0] ? 64'hB5026F5AA96619E9 : 64'h0);
+    assign y3 = y2 ^ ((y2 >> 29) & 64'h5555555555555555);
+    assign y4 = y3 ^ ((y3 << 17) & 64'h71D67FFFEDA60000);
+    assign y5 = y4 ^ ((y4 << 37) & 64'hFFF7EEE000000000);
+
+end
+
+always_comb begin
+    state_next = STATE_SEED;
+
+    mt_save_next = mt_save_reg;
+    mti_next = mti_reg;
+
+    mt_wr_data = y2;
+    mt_wr_ptr = mti_reg;
+    mt_wr_en = '0;
+
+    mt_rd_a_ptr_next = mt_rd_a_ptr_reg;
+    mt_rd_b_ptr_next = mt_rd_b_ptr_reg;
+
+    product_next = product_reg;
+    factor1_next = factor1_reg;
+    factor2_next = factor2_reg;
+    mul_done_next = mul_done_reg;
+
+    m_axis_tdata_next = m_axis_tdata_reg;
+    m_axis_tvalid_next = m_axis_tvalid_reg && !m_axis.tready;
+
+    case (state_reg)
+        STATE_SEED: begin
+            mt_save_next = product_reg + MT_W'(mti_reg);
+            mt_wr_data = mt_save_next;
+            mt_rd_b_ptr_next = MT_M;
+            if (mul_done_reg) begin
+                product_next = '0;
+                if (MT_W == 32) begin
+                    factor1_next = mt_save_next ^ (mt_save_next >> 30);
+                    /* verilator lint_off WIDTHEXPAND */
+                    factor2_next = 32'd1812433253;
+                    /* verilator lint_on WIDTHEXPAND */
+                end else begin
+                    factor1_next = mt_save_next ^ (mt_save_next >> 62);
+                    /* verilator lint_off WIDTHTRUNC */
+                    factor2_next = 64'd6364136223846793005;
+                    /* verilator lint_on WIDTHTRUNC */
+                end
+                mul_done_next = 1'b0;
+                if (mti_reg < CL_MT_N'(MT_N)) begin
+                    product_next = '0;
+                    mul_done_next = 1'b0;
+                    mt_wr_data = mt_save_next;
+                    mt_wr_ptr = mti_reg;
+                    mt_wr_en = 1'b1;
+                    mti_next = mti_reg + 1;
+                    mt_rd_a_ptr_next = '0;
+                    state_next = STATE_SEED;
+                end else begin
+                    mti_next = '0;
+                    mt_rd_a_ptr_next = 1;
+                    mt_rd_b_ptr_next = MT_M;
+                    mt_save_next = mt_rd_a_data_reg;
+                    state_next = STATE_RUN;
+                end
+            end else begin
+                factor1_next = factor1_reg << 1;
+                factor2_next = factor2_reg >> 1;
+                mul_done_next = factor2_reg[8:1] == 0;
+                if (factor2_reg[0]) begin
+                    product_next = product_reg + factor1_reg;
+                end
+                state_next = STATE_SEED;
+            end
+        end
+        STATE_RUN: begin
+            // idle state
+            if (m_axis.tready) begin
+                if (mti_reg < CL_MT_N'(MT_N-1))
+                    mti_next = mti_reg + 1;
+                else
+                    mti_next = '0;
+
+                if (mt_rd_a_ptr_reg < CL_MT_N'(MT_N-1))
+                    mt_rd_a_ptr_next = mt_rd_a_ptr_reg + 1;
+                else
+                    mt_rd_a_ptr_next = '0;
+
+                if (mt_rd_b_ptr_reg < CL_MT_N'(MT_N-1))
+                    mt_rd_b_ptr_next = mt_rd_b_ptr_reg + 1;
+                else
+                    mt_rd_b_ptr_next = '0;
+
+                mt_save_next = mt_rd_a_data_reg;
+
+                if (MT_W == 32) begin
+                    m_axis_tdata_next = y5 ^ (y5 >> 18);
+                end else begin
+                    m_axis_tdata_next = y5 ^ (y5 >> 43);
+                end
+                m_axis_tvalid_next = 1'b1;
+
+                mt_wr_data = y2;
+                mt_wr_ptr = mti_reg;
+                mt_wr_en = 1'b1;
+            end
+            state_next = STATE_RUN;
+        end
+    endcase
+
+    if (seed_start) begin
+        product_next = seed_val;
+        mti_next = '0;
+        mul_done_next = 1'b1;
+        state_next = STATE_SEED;
+    end
+end
+
+always_ff @(posedge clk) begin
+    state_reg <= state_next;
+
+    mt_save_reg <= mt_save_next;
+    mti_reg <= mti_next;
+
+    mt_rd_a_ptr_reg <= mt_rd_a_ptr_next;
+    mt_rd_b_ptr_reg <= mt_rd_b_ptr_next;
+
+    product_reg <= product_next;
+    factor1_reg <= factor1_next;
+    factor2_reg <= factor2_next;
+    mul_done_reg <= mul_done_next;
+
+    m_axis_tdata_reg <= m_axis_tdata_next;
+    m_axis_tvalid_reg <= m_axis_tvalid_next;
+
+    busy_reg <= state_next != STATE_RUN;
+
+    if (mt_wr_en) begin
+        mt_ram[mt_wr_ptr] <= mt_wr_data;
+    end
+
+    mt_rd_a_data_reg <= mt_ram[mt_rd_a_ptr_next];
+    mt_rd_b_data_reg <= mt_ram[mt_rd_b_ptr_next];
+
+    if (rst) begin
+        state_reg <= STATE_SEED;
+        mti_reg <= '0;
+        mt_rd_a_ptr_reg <= '0;
+        mt_rd_b_ptr_reg <= '0;
+        product_reg <= INIT_SEED;
+        factor1_reg <= '0;
+        factor2_reg <= '0;
+        mul_done_reg <= 1'b1;
+        m_axis_tdata_reg <= '0;
+        m_axis_tvalid_reg <= 1'b0;
+        busy_reg <= 1'b0;
+    end
+end
+
+endmodule
+
+`resetall
diff --git a/src/math/tb/mt19937.py b/src/math/tb/mt19937.py
new file mode 100644
index 0000000..b7e5e8c
--- /dev/null
+++ b/src/math/tb/mt19937.py
@@ -0,0 +1,116 @@
+#!/usr/bin/env python3
+# SPDX-License-Identifier: CERN-OHL-S-2.0
+"""
+
+Copyright (c) 2014-2025 FPGA Ninja, LLC
+
+Authors:
+- Alex Forencich
+
+This is a python implementation of MT19937 from
+http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
+http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
+
+"""
+
+class mt19937(object):
+    def __init__(self):
+        self.mt = [0]*624
+        self.mti = 625
+
+    def seed(self, seed):
+        self.mt[0] = seed & 0xffffffff
+        for i in range(1,624):
+            self.mt[i] = (1812433253 * (self.mt[i-1] ^ (self.mt[i-1] >> 30)) + i) & 0xffffffff
+        self.mti = 624
+
+    def init_by_array(self, key):
+        self.seed(19650218)
+        i = 1
+        j = 0
+        k = max(624, len(key))
+        for ki in range(k):
+            self.mt[i] = ((self.mt[i] ^ ((self.mt[i-1] ^ (self.mt[i-1] >> 30)) * 1664525)) + key[j] + j) & 0xffffffff
+            i += 1
+            j += 1
+            if i >= 624:
+                self.mt[0] = self.mt[623]
+                i = 1
+            if j >= len(key):
+                j = 0
+        for ki in range(624):
+            self.mt[i] = ((self.mt[i] ^ ((self.mt[i-1] ^ (self.mt[i-1] >> 30)) * 1566083941)) - i) & 0xffffffff
+            i += 1
+            if i >= 624:
+                self.mt[0] = self.mt[623]
+                i = 1
+        self.mt[0] = 0x80000000
+
+    def int32(self):
+        if self.mti >= 624:
+            if self.mti == 625:
+                self.seed(5489)
+
+            for k in range(623):
+                y = (self.mt[k] & 0x80000000) | (self.mt[k+1] & 0x7fffffff)
+                if k < 624 - 397:
+                    self.mt[k] = self.mt[k+397] ^ (y >> 1) ^ (0x9908b0df if y & 1 else 0)
+                else:
+                    self.mt[k] = self.mt[k+397-624] ^ (y >> 1) ^ (0x9908b0df if y & 1 else 0)
+
+            y = (self.mt[623] & 0x80000000) | (self.mt[0] & 0x7fffffff)
+            self.mt[623] = self.mt[396] ^ (y >> 1) ^ (0x9908b0df if y & 1 else 0)
+            self.mti = 0
+
+        y = self.mt[self.mti]
+        self.mti += 1
+
+        y ^= (y >> 11)
+        y ^= (y << 7) & 0x9d2c5680
+        y ^= (y << 15) & 0xefc60000
+        y ^= (y >> 18)
+
+        return y
+
+    def int32b(self):
+        if self.mti == 625:
+            self.seed(5489)
+
+        k = self.mti
+
+        if k == 624:
+            k = 0
+            self.mti = 0
+
+        if k == 623:
+            y = (self.mt[623] & 0x80000000) | (self.mt[0] & 0x7fffffff)
+            self.mt[623] = self.mt[396] ^ (y >> 1) ^ (0x9908b0df if y & 1 else 0)
+        else:
+            y = (self.mt[k] & 0x80000000) | (self.mt[k+1] & 0x7fffffff)
+            if k < 624 - 397:
+                self.mt[k] = self.mt[k+397] ^ (y >> 1) ^ (0x9908b0df if y & 1 else 0)
+            else:
+                self.mt[k] = self.mt[k+397-624] ^ (y >> 1) ^ (0x9908b0df if y & 1 else 0)
+
+        y = self.mt[self.mti]
+        self.mti += 1
+
+        y ^= (y >> 11)
+        y ^= (y << 7) & 0x9d2c5680
+        y ^= (y << 15) & 0xefc60000
+        y ^= (y >> 18)
+
+        return y
+
+if __name__ == '__main__':
+    mt = mt19937()
+    mt.init_by_array([0x123, 0x234, 0x345, 0x456])
+    print("1000 outputs of int32")
+    s=''
+    for i in range(1000):
+        s += "%10lu " % mt.int32()
+        if i % 5 == 4:
+            print(s)
+            s = ''
+    if len(s) > 0:
+        print(s)
diff --git a/src/math/tb/mt19937_64.py b/src/math/tb/mt19937_64.py
new file mode 100644
index 0000000..eb7aba0
--- /dev/null
+++ b/src/math/tb/mt19937_64.py
@@ -0,0 +1,116 @@
+#!/usr/bin/env python3
+# SPDX-License-Identifier: CERN-OHL-S-2.0
+"""
+
+Copyright (c) 2014-2025 FPGA Ninja, LLC
+
+Authors:
+- Alex Forencich
+
+This is a python implementation of MT19937-64 from
+http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html
+http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c
+
+"""
+
+class mt19937_64(object):
+    def __init__(self):
+        self.mt = [0]*312
+        self.mti = 313
+
+    def seed(self, seed):
+        self.mt[0] = seed & 0xffffffffffffffff
+        for i in range(1,312):
+            self.mt[i] = (6364136223846793005 * (self.mt[i-1] ^ (self.mt[i-1] >> 62)) + i) & 0xffffffffffffffff
+        self.mti = 312
+
+    def init_by_array(self, key):
+        self.seed(19650218)
+        i = 1
+        j = 0
+        k = max(312, len(key))
+        for ki in range(k):
+            self.mt[i] = ((self.mt[i] ^ ((self.mt[i-1] ^ (self.mt[i-1] >> 62)) * 3935559000370003845)) + key[j] + j) & 0xffffffffffffffff
+            i += 1
+            j += 1
+            if i >= 312:
+                self.mt[0] = self.mt[311]
+                i = 1
+            if j >= len(key):
+                j = 0
+        for ki in range(312):
+            self.mt[i] = ((self.mt[i] ^ ((self.mt[i-1] ^ (self.mt[i-1] >> 62)) * 2862933555777941757)) - i) & 0xffffffffffffffff
+            i += 1
+            if i >= 312:
+                self.mt[0] = self.mt[311]
+                i = 1
+        self.mt[0] = 1 << 63
+
+    def int64(self):
+        if self.mti >= 312:
+            if self.mti == 313:
+                self.seed(5489)
+
+            for k in range(311):
+                y = (self.mt[k] & 0xFFFFFFFF80000000) | (self.mt[k+1] & 0x7fffffff)
+                if k < 312 - 156:
+                    self.mt[k] = self.mt[k+156] ^ (y >> 1) ^ (0xB5026F5AA96619E9 if y & 1 else 0)
+                else:
+                    self.mt[k] = self.mt[k+156-624] ^ (y >> 1) ^ (0xB5026F5AA96619E9 if y & 1 else 0)
+
+            y = (self.mt[311] & 0xFFFFFFFF80000000) | (self.mt[0] & 0x7fffffff)
+            self.mt[311] = self.mt[155] ^ (y >> 1) ^ (0xB5026F5AA96619E9 if y & 1 else 0)
+            self.mti = 0
+
+        y = self.mt[self.mti]
+        self.mti += 1
+
+        y ^= (y >> 29) & 0x5555555555555555
+        y ^= (y << 17) & 0x71D67FFFEDA60000
+        y ^= (y << 37) & 0xFFF7EEE000000000
+        y ^= (y >> 43)
+
+        return y
+
+    def int64b(self):
+        if self.mti == 313:
+            self.seed(5489)
+
+        k = self.mti
+
+        if k == 312:
+            k = 0
+            self.mti = 0
+
+        if k == 311:
+            y = (self.mt[311] & 0xFFFFFFFF80000000) | (self.mt[0] & 0x7fffffff)
+            self.mt[311] = self.mt[155] ^ (y >> 1) ^ (0xB5026F5AA96619E9 if y & 1 else 0)
+        else:
+            y = (self.mt[k] & 0xFFFFFFFF80000000) | (self.mt[k+1] & 0x7fffffff)
+            if k < 312 - 156:
+                self.mt[k] = self.mt[k+156] ^ (y >> 1) ^ (0xB5026F5AA96619E9 if y & 1 else 0)
+            else:
+                self.mt[k] = self.mt[k+156-624] ^ (y >> 1) ^ (0xB5026F5AA96619E9 if y & 1 else 0)
+
+        y = self.mt[self.mti]
+        self.mti += 1
+
+        y ^= (y >> 29) & 0x5555555555555555
+        y ^= (y << 17) & 0x71D67FFFEDA60000
+        y ^= (y << 37) & 0xFFF7EEE000000000
+        y ^= (y >> 43)
+
+        return y
+
+if __name__ == '__main__':
+    mt = mt19937_64()
+    mt.init_by_array([0x12345, 0x23456, 0x34567, 0x45678])
+    print("1000 outputs of int64")
+    s=''
+    for i in range(1000):
+        s += "%10lu " % mt.int64b()
+        if i % 5 == 4:
+            print(s)
+            s = ''
+    if len(s) > 0:
+        print(s)
diff --git a/src/math/tb/taxi_mt19937/Makefile b/src/math/tb/taxi_mt19937/Makefile
new file mode 100644
index 0000000..6203f2d
--- /dev/null
+++ b/src/math/tb/taxi_mt19937/Makefile
@@ -0,0 +1,52 @@
+# SPDX-License-Identifier: CERN-OHL-S-2.0
+#
+# Copyright (c) 2021-2025 FPGA Ninja, LLC
+#
+# Authors:
+# - Alex Forencich
+
+TOPLEVEL_LANG = verilog
+
+SIM ?= verilator
+WAVES ?= 0
+
+COCOTB_HDL_TIMEUNIT = 1ns
+COCOTB_HDL_TIMEPRECISION = 1ps
+
+RTL_DIR = ../../rtl
+LIB_DIR = ../../lib
+TAXI_SRC_DIR = $(LIB_DIR)/taxi/src
+
+DUT      = taxi_mt19937
+COCOTB_TEST_MODULES = test_$(DUT)
+COCOTB_TOPLEVEL     = test_$(DUT)
+MODULE   = $(COCOTB_TEST_MODULES)
+TOPLEVEL = $(COCOTB_TOPLEVEL)
+VERILOG_SOURCES += $(COCOTB_TOPLEVEL).sv
+VERILOG_SOURCES += $(RTL_DIR)/$(DUT).sv
+VERILOG_SOURCES += $(TAXI_SRC_DIR)/axis/rtl/taxi_axis_if.sv
+
+# handle file list files
+process_f_file = $(call process_f_files,$(addprefix $(dir $1),$(shell cat $1)))
+process_f_files = $(foreach f,$1,$(if $(filter %.f,$f),$(call process_f_file,$f),$f))
+uniq_base = $(if $1,$(call uniq_base,$(foreach f,$1,$(if $(filter-out $(notdir $(lastword $1)),$(notdir $f)),$f,))) $(lastword $1))
+VERILOG_SOURCES := $(call uniq_base,$(call process_f_files,$(VERILOG_SOURCES)))
+
+# module parameters
+export PARAM_MT_W := 32
+export PARAM_INIT_SEED := 5489
+
+ifeq ($(SIM), icarus)
+	PLUSARGS += -fst
+
+	COMPILE_ARGS += $(foreach v,$(filter PARAM_%,$(.VARIABLES)),-P $(COCOTB_TOPLEVEL).$(subst PARAM_,,$(v))=$($(v)))
+else ifeq ($(SIM), verilator)
+	COMPILE_ARGS += $(foreach v,$(filter PARAM_%,$(.VARIABLES)),-G$(subst PARAM_,,$(v))=$($(v)))
+
+	ifeq ($(WAVES), 1)
+		COMPILE_ARGS += --trace-fst
+		VERILATOR_TRACE = 1
+	endif
+endif
+
+include $(shell cocotb-config --makefiles)/Makefile.sim
diff --git a/src/math/tb/taxi_mt19937/mt19937.py b/src/math/tb/taxi_mt19937/mt19937.py
new file mode 120000
index 0000000..00a3e5c
--- /dev/null
+++ b/src/math/tb/taxi_mt19937/mt19937.py
@@ -0,0 +1 @@
+../mt19937.py
\ No newline at end of file
diff --git a/src/math/tb/taxi_mt19937/mt19937_64.py b/src/math/tb/taxi_mt19937/mt19937_64.py
new file mode 120000
index 0000000..8d473c1
--- /dev/null
+++ b/src/math/tb/taxi_mt19937/mt19937_64.py
@@ -0,0 +1 @@
+../mt19937_64.py
\ No newline at end of file
diff --git a/src/math/tb/taxi_mt19937/test_taxi_mt19937.py b/src/math/tb/taxi_mt19937/test_taxi_mt19937.py
new file mode 100644
index 0000000..7bab4a1
--- /dev/null
+++ b/src/math/tb/taxi_mt19937/test_taxi_mt19937.py
@@ -0,0 +1,168 @@
+#!/usr/bin/env python3
+# SPDX-License-Identifier: CERN-OHL-S-2.0
+"""
+
+Copyright (c) 2021-2025 FPGA Ninja, LLC
+
+Authors:
+- Alex Forencich
+
+"""
+
+import itertools
+import logging
+import os
+import sys
+import pytest
+
+import cocotb_test.simulator
+
+import cocotb
+from cocotb.clock import Clock
+from cocotb.triggers import RisingEdge
+from cocotb.regression import TestFactory
+
+from cocotbext.axi import AxiStreamBus, AxiStreamSink
+
+try:
+    import mt19937, mt19937_64
+except ImportError:
+    # attempt import from current directory
+    sys.path.insert(0, os.path.join(os.path.dirname(__file__)))
+    try:
+        import mt19937, mt19937_64
+    finally:
+        del sys.path[0]
+
+
+class TB(object):
+    def __init__(self, dut):
+        self.dut = dut
+
+        self.log = logging.getLogger("cocotb.tb")
+        self.log.setLevel(logging.DEBUG)
+
+        cocotb.start_soon(Clock(dut.clk, 10, units="ns").start())
+
+        self.sink = AxiStreamSink(AxiStreamBus.from_entity(dut.m_axis), dut.clk, dut.rst)
+
+    def set_backpressure_generator(self, generator=None):
+        if generator:
+            self.sink.set_pause_generator(generator())
+
+    async def reset(self):
+        self.dut.rst.setimmediatevalue(0)
+        await RisingEdge(self.dut.clk)
+        await RisingEdge(self.dut.clk)
+        self.dut.rst.value = 1
+        await RisingEdge(self.dut.clk)
+        await RisingEdge(self.dut.clk)
+        self.dut.rst.value = 0
+        await RisingEdge(self.dut.clk)
+        await RisingEdge(self.dut.clk)
+
+
+async def run_test(dut, seed=None, backpressure_inserter=None):
+
+    tb = TB(dut)
+
+    await tb.reset()
+
+    tb.set_backpressure_generator(backpressure_inserter)
+
+    if dut.MT_W.value == 32:
+        mt = mt19937.mt19937()
+    else:
+        mt = mt19937_64.mt19937_64()
+
+    if seed is not None:
+        mt.seed(seed)
+
+        await RisingEdge(dut.clk)
+        dut.seed_val.value = seed
+        dut.seed_start.value = 1
+        await RisingEdge(dut.clk)
+        dut.seed_start.value = 0
+        await RisingEdge(dut.clk)
+
+    tb.sink.clear()
+
+    for i in range(1000):
+        frame = await tb.sink.recv()
+        if dut.MT_W.value == 32:
+            ref = mt.int32()
+        else:
+            ref = mt.int64()
+        assert frame.tdata[0] == ref
+
+    await RisingEdge(dut.clk)
+    await RisingEdge(dut.clk)
+
+
+def cycle_pause():
+    return itertools.cycle([1, 1, 1, 0])
+
+
+if getattr(cocotb, 'top', None) is not None:
+
+    factory = TestFactory(run_test)
+    factory.add_option("backpressure_inserter", [None, cycle_pause])
+    factory.add_option("seed", [None, 0x12345678])
+    factory.generate_tests()
+
+
+# cocotb-test
+
+tests_dir = os.path.dirname(__file__)
+rtl_dir = os.path.abspath(os.path.join(tests_dir, '..', '..', 'rtl'))
+lib_dir = os.path.abspath(os.path.join(tests_dir, '..', '..', 'lib'))
+taxi_src_dir = os.path.abspath(os.path.join(lib_dir, 'taxi', 'src'))
+
+
+def process_f_files(files):
+    lst = {}
+    for f in files:
+        if f[-2:].lower() == '.f':
+            with open(f, 'r') as fp:
+                l = fp.read().split()
+            for f in process_f_files([os.path.join(os.path.dirname(f), x) for x in l]):
+                lst[os.path.basename(f)] = f
+        else:
+            lst[os.path.basename(f)] = f
+    return list(lst.values())
+
+
+@pytest.mark.parametrize("mt_w", [32, 64])
+def test_taxi_mt19937(request, mt_w):
+    dut = "taxi_mt19937"
+    module = os.path.splitext(os.path.basename(__file__))[0]
+    toplevel = module
+
+    verilog_sources = [
+        os.path.join(tests_dir, f"{toplevel}.sv"),
+        os.path.join(rtl_dir, f"{dut}.sv"),
+        os.path.join(taxi_src_dir, "axis", "rtl", "taxi_axis_if.sv"),
+    ]
+
+    verilog_sources = process_f_files(verilog_sources)
+
+    parameters = {}
+
+    parameters['MT_W'] = mt_w
+    parameters['INIT_SEED'] = 5489
+
+    extra_env = {f'PARAM_{k}': str(v) for k, v in parameters.items()}
+
+    sim_build = os.path.join(tests_dir, "sim_build",
+        request.node.name.replace('[', '-').replace(']', ''))
+
+    cocotb_test.simulator.run(
+        simulator="verilator",
+        python_search=[tests_dir],
+        verilog_sources=verilog_sources,
+        toplevel=toplevel,
+        module=module,
+        parameters=parameters,
+        sim_build=sim_build,
+        extra_env=extra_env,
+    )
diff --git a/src/math/tb/taxi_mt19937/test_taxi_mt19937.sv b/src/math/tb/taxi_mt19937/test_taxi_mt19937.sv
new file mode 100644
index 0000000..5f8372f
--- /dev/null
+++ b/src/math/tb/taxi_mt19937/test_taxi_mt19937.sv
@@ -0,0 +1,68 @@
+// SPDX-License-Identifier: CERN-OHL-S-2.0
+/*
+
+Copyright (c) 2025 FPGA Ninja, LLC
+
+Authors:
+- Alex Forencich
+
+*/
+
+`resetall
+`timescale 1ns / 1ps
+`default_nettype none
+
+/*
+ * MT19937/MT19937-64 Mersenne Twister PRNG testbench
+ */
+module test_taxi_mt19937 #(
+    /* verilator lint_off WIDTHTRUNC */
+    /* verilator lint_off WIDTHEXPAND */
+    parameter integer MT_W = 32,
+    parameter logic [MT_W-1:0] INIT_SEED = 5489
+    /* verilator lint_on WIDTHTRUNC */
+    /* verilator lint_on WIDTHEXPAND */
+)
+();
+
+logic clk;
+logic rst;
+
+taxi_axis_if #(
+    .DATA_W(MT_W),
+    .KEEP_W(1)
+) m_axis();
+
+wire busy;
+
+logic [MT_W-1:0] seed_val;
+logic seed_start;
+
+taxi_mt19937 #(
+    .MT_W(MT_W),
+    .INIT_SEED(INIT_SEED)
+)
+uut (
+    .clk(clk),
+    .rst(rst),
+
+    /*
+     * AXI output (source)
+     */
+    .m_axis(m_axis),
+
+    /*
+        * Status
+        */
+    .busy(busy),
+
+    /*
+        * Configuration
+        */
+    .seed_val(seed_val),
+    .seed_start(seed_start)
+);
+
+endmodule
+
+`resetall