diff --git a/brainpy/_src/math/__init__.py b/brainpy/_src/math/__init__.py
index 3128c5e67..c730029cd 100644
--- a/brainpy/_src/math/__init__.py
+++ b/brainpy/_src/math/__init__.py
@@ -44,7 +44,7 @@
 from .compat_numpy import *
 from .compat_tensorflow import *
 from .others import *
-from . import random, linalg, fft
+from . import random, linalg, fft, taichi_random
 
 # operators
 from .op_register import *
diff --git a/brainpy/_src/math/jitconn/_event_matvec_taichi.py b/brainpy/_src/math/jitconn/_event_matvec_taichi.py
index 7c402b912..9bf018720 100644
--- a/brainpy/_src/math/jitconn/_event_matvec_taichi.py
+++ b/brainpy/_src/math/jitconn/_event_matvec_taichi.py
@@ -51,19 +51,15 @@ def _event_mv_prob_homo_bool_cpu(
   
   # ti.loop_config(serialize=True)
   for i_col in range(num_col):
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_row = uniform_int_distribution(result, 1, clen_value)
     event = events[i_col]
     while i_row < num_row:
       if event:
         out[i_row] += weight_value
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_row += uniform_int_distribution(result, 1, clen_value)
 
 
@@ -87,19 +83,15 @@ def _event_mv_prob_homo_bool_gpu(
     i_col = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_row = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_row < num_row:
       if event:
         out[i_row] += weight_value
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_row += uniform_int_distribution(result, 1, clen_value) * 32
 
 
@@ -120,20 +112,16 @@ def _event_mv_prob_homo_outdim_parallel_bool_cpu(
   
   # ti.loop_config(serialize=True)
   for i_row in range(num_row):
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_col = uniform_int_distribution(result, 1, clen_value)
     event = events[i_col]
     while i_col < num_col:
       if event:
         r += weight_value
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_col += uniform_int_distribution(result, 1, clen_value)
     out[i_row] = r
 
@@ -158,20 +146,16 @@ def _event_mv_prob_homo_outdim_parallel_bool_gpu(
     i_row = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_col = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_col < num_col:
       if event:
         r += weight_value
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_col += uniform_int_distribution(result, 1, clen_value)  * 32
     out[i_row] += r
 
@@ -193,19 +177,15 @@ def _event_mv_prob_homo_cpu(
 
   # ti.loop_config(serialize=True)
   for i_col in range(num_col):
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_row = uniform_int_distribution(result, 1, clen_value)
     event = events[i_col]
     while i_row < num_row:
       if event > 0.:
         out[i_row] += weight_value
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_row += uniform_int_distribution(result, 1, clen_value)
 
 
@@ -229,19 +209,15 @@ def _event_mv_prob_homo_gpu(
     i_col = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_row = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_row < num_row:
       if event > 0.:
         out[i_row] += weight_value
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_row += uniform_int_distribution(result, 1, clen_value) * 32
 
 
@@ -262,20 +238,16 @@ def _event_mv_prob_homo_outdim_parallel_cpu(
 
   # ti.loop_config(serialize=True)
   for i_row in range(num_row):
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_col = uniform_int_distribution(result, 1, clen_value)
     event = events[i_col]
     while i_col < num_col:
       if event > 0.:
         r += weight_value
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_col += uniform_int_distribution(result, 1, clen_value)
     out[i_row] = r
 
@@ -300,20 +272,16 @@ def _event_mv_prob_homo_outdim_parallel_gpu(
     i_row = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_col = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_col < num_col:
       if event > 0.:
         r += weight_value
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_col += uniform_int_distribution(result, 1, clen_value) * 32
     out[i_row] += r
 
@@ -563,20 +531,16 @@ def _event_mv_prob_uniform_bool_cpu(
 
   # ti.loop_config(serialize=True)
   for i_col in range(num_col):
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_row = uniform_int_distribution(result, 1, clen_value)
     event = events[i_col]
     while i_row < num_row:
       if event:
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         out[i_row] += uniform_real_distribution(result, w_min_value, w_max_value)
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_row += uniform_int_distribution(result, 1, clen_value)
 
 
@@ -603,20 +567,16 @@ def _event_mv_prob_uniform_bool_gpu(
     i_col = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_row = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_row < num_row:
       if event:
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         out[i_row] += uniform_real_distribution(result, w_min_value, w_max_value)
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_row += uniform_int_distribution(result, 1, clen_value) * 32
 
 
@@ -639,21 +599,17 @@ def _event_mv_prob_uniform_outdim_parallel_bool_cpu(
 
   # ti.loop_config(serialize=True)
   for i_row in range(num_row):
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_col = uniform_int_distribution(result, 1, clen_value)
     event = events[i_col]
     while i_col < num_col:
       if event:
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         r += uniform_real_distribution(result, w_min_value, w_max_value)
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_col += uniform_int_distribution(result, 1, clen_value)
     out[i_row] = r
 
@@ -680,21 +636,17 @@ def _event_mv_prob_uniform_outdim_parallel_bool_gpu(
     i_row = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_col = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_col < num_col:
       if event:
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         r += uniform_real_distribution(result, w_min_value, w_max_value)
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_col += uniform_int_distribution(result, 1, clen_value) * 32
     out[i_row] += r
 
@@ -719,20 +671,16 @@ def _event_mv_prob_uniform_cpu(
 
   # ti.loop_config(serialize=True)
   for i_col in range(num_col):
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_row = uniform_int_distribution(result, 1, clen_value)
     event = events[i_col]
     while i_row < num_row:
       if event > 0.:
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         out[i_row] += uniform_real_distribution(result, w_min_value, w_max_value)
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_row += uniform_int_distribution(result, 1, clen_value)
 
 
@@ -759,20 +707,16 @@ def _event_mv_prob_uniform_gpu(
     i_col = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_row = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_row < num_row:
       if event > 0.:
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         out[i_row] += uniform_real_distribution(result, w_min_value, w_max_value)
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_row += uniform_int_distribution(result, 1, clen_value) * 32
 
 
@@ -795,21 +739,17 @@ def _event_mv_prob_uniform_outdim_parallel_cpu(
 
   # ti.loop_config(serialize=True)
   for i_row in range(num_row):
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_col = uniform_int_distribution(result, 1, clen_value)
     event = events[i_col]
     while i_col < num_col:
       if event > 0.:
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         r += uniform_real_distribution(result, w_min_value, w_max_value)
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_col += uniform_int_distribution(result, 1, clen_value)
     out[i_row] = r
 
@@ -836,21 +776,17 @@ def _event_mv_prob_uniform_outdim_parallel_gpu(
     i_row = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+    seeds, result = random_generator(seeds)
     i_col = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_col < num_col:
       if event > 0.:
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         r += uniform_real_distribution(result, w_min_value, w_max_value)
-      s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+      seeds, result = random_generator(seeds)
       i_col += uniform_int_distribution(result, 1, clen_value) * 32
     out[i_row] += r
 
@@ -1094,20 +1030,15 @@ def _event_mv_prob_normal_bool_cpu(
 
   # ti.loop_config(serialize=True)
   for i_col in range(num_col):
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result1 = ti.f32(0.)
-    result2 = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+    seeds, result1 = random_generator(seeds)
     i_row = uniform_int_distribution(result1, 1, clen_value)
     event = events[i_col]
     while i_row < num_row:
       if event:
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-        s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
+        seeds, result2 = random_generator(seeds)
         out[i_row] += normal_distribution(result1, result2, w_mu_value, w_sigma_value)
       i_row += uniform_int_distribution(result1, 1, clen_value)
 
@@ -1134,20 +1065,15 @@ def _event_mv_prob_normal_bool_gpu(
     i_col = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result1 = ti.f32(0.)
-    result2 = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+    seeds, result1 = random_generator(seeds)
     i_row = uniform_int_distribution(result1, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_row < num_row:
       if event:
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-        s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
+        seeds, result2 = random_generator(seeds)
         out[i_row] += normal_distribution(result1, result2, w_mu_value, w_sigma_value)
       i_row += uniform_int_distribution(result1, 1, clen_value) * 32
 
@@ -1171,21 +1097,16 @@ def _event_mv_prob_normal_outdim_parallel_bool_cpu(
 
   # ti.loop_config(serialize=True)
   for i_row in range(num_row):
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result1 = ti.f32(0.)
-    result2 = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+    seeds, result1 = random_generator(seeds)
     i_col = uniform_int_distribution(result1, 1, clen_value)
     event = events[i_col]
     while i_col < num_col:
       if event:
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-        s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
+        seeds, result2 = random_generator(seeds)
         r += normal_distribution(result1, result2, w_mu_value, w_sigma_value)
       i_col += uniform_int_distribution(result1, 1, clen_value)
     out[i_row] = r
@@ -1213,21 +1134,16 @@ def _event_mv_prob_normal_outdim_parallel_bool_gpu(
     i_row = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result1 = ti.f32(0.)
-    result2 = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+    seeds, result1 = random_generator(seeds)
     i_col = uniform_int_distribution(result1, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_col < num_col:
       if event:
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-        s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
+        seeds, result2 = random_generator(seeds)
         r += normal_distribution(result1, result2, w_mu_value, w_sigma_value)
       i_col += uniform_int_distribution(result1, 1, clen_value) * 32
     out[i_row] += r
@@ -1252,20 +1168,15 @@ def _event_mv_prob_normal_cpu(
 
   # ti.loop_config(serialize=True)
   for i_col in range(num_col):
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result1 = ti.f32(0.)
-    result2 = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+    seeds, result1 = random_generator(seeds)
     i_row = uniform_int_distribution(result1, 1, clen_value)
     event = events[i_col]
     while i_row < num_row:
       if event > 0.:
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-        s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
+        seeds, result2 = random_generator(seeds)
         out[i_row] += normal_distribution(result1, result2, w_mu_value, w_sigma_value)
       i_row += uniform_int_distribution(result1, 1, clen_value)
 
@@ -1292,20 +1203,15 @@ def _event_mv_prob_normal_gpu(
     i_col = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_col
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result1 = ti.f32(0.)
-    result2 = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
     
-    s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+    seeds, result1 = random_generator(seeds)
     i_row = uniform_int_distribution(result1, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_row < num_row:
       if event > 0.:
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-        s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
+        seeds, result2 = random_generator(seeds)
         out[i_row] += normal_distribution(result1, result2, w_mu_value, w_sigma_value)
       i_row += uniform_int_distribution(result1, 1, clen_value) * 32
 
@@ -1329,21 +1235,16 @@ def _event_mv_prob_normal_outdim_parallel_cpu(
 
   # ti.loop_config(serialize=True)
   for i_row in range(num_row):
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result1 = ti.f32(0.)
-    result2 = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+    seeds, result1 = random_generator(seeds)
     i_col = uniform_int_distribution(result1, 1, clen_value)
     event = events[i_col]
     while i_col < num_col:
       if event > 0.:
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-        s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
+        seeds, result2 = random_generator(seeds)
         r += normal_distribution(result1, result2, w_mu_value, w_sigma_value)
       i_col += uniform_int_distribution(result1, 1, clen_value)
     out[i_row] = r
@@ -1371,21 +1272,16 @@ def _event_mv_prob_normal_outdim_parallel_gpu(
     i_row = i >> 5
     index = i & 31
     
-    s1 = seed_value + 1 + i_row
-    s2 = seed_value + 7
-    s3 = seed_value + 15
-    b = ti.u32(0)
-    result1 = ti.f32(0.)
-    result2 = ti.f32(0.)
+    seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
     r = 0.
     
-    s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+    seeds, result1 = random_generator(seeds)
     i_col = uniform_int_distribution(result1, 1, clen_value) + avg_num_uniform * index
     event = events[i_col]
     while i_col < num_col:
       if event > 0.:
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-        s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
+        seeds, result2 = random_generator(seeds)
         r += normal_distribution(result1, result2, w_mu_value, w_sigma_value)
       i_col += uniform_int_distribution(result1, 1, clen_value) * 32
     out[i_row] += r
diff --git a/brainpy/_src/math/jitconn/_matvec_taichi.py b/brainpy/_src/math/jitconn/_matvec_taichi.py
index 35a7e31f8..7e3f8d29e 100644
--- a/brainpy/_src/math/jitconn/_matvec_taichi.py
+++ b/brainpy/_src/math/jitconn/_matvec_taichi.py
@@ -32,6 +32,7 @@
 # ) -> ti.i32:
 #     return uniform_int_distribution(random_generator(seed), 1, clen)
 
+
 @ti.kernel
 def _mv_prob_homo_cpu(
     vector: ti.types.ndarray(ndim=1),
@@ -49,19 +50,15 @@ def _mv_prob_homo_cpu(
 
     # ti.loop_config(serialize=True)
     for i_col in range(num_col):
-        s1 = seed_value + 1 + i_col
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
-        result = ti.f32(0.)
+        seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
 
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         i_row = uniform_int_distribution(result, 1, clen_value)
         v = vector[i_col] * weight_value
         while i_row < num_row:
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             out[i_row] += uniform_int_distribution(result, 1, clen_value) * v
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             i_row += uniform_int_distribution(result, 1, clen_value)
 
 @ti.kernel
@@ -84,19 +81,15 @@ def _mv_prob_homo_gpu(
         i_col = i >> 5
         index = i & 31
         
-        s1 = seed_value + 1 + i_col
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
-        result = ti.f32(0.)
+        seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
         
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         i_row = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
         v = vector[i_col] * weight_value
         while i_row < num_row:
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             out[i_row] += uniform_int_distribution(result, 1, clen_value) * v
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             i_row += uniform_int_distribution(result, 1, clen_value) * 32
 
 
@@ -117,18 +110,14 @@ def _mv_prob_homo_outdim_parallel_cpu(
     
     # ti.loop_config(serialize=True)
     for i_row in range(num_row):
-        s1 = seed_value + 1 + i_row
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
+        seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
         r = 0.
-        result = ti.f32(0.)
 
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         i_col = uniform_int_distribution(result, 1, clen_value)
         while i_col < num_col:
             r += vector[i_col]
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             i_col += uniform_int_distribution(result, 1, clen_value)
         out[i_row] = r * weight_value
 
@@ -152,18 +141,14 @@ def _mv_prob_homo_outdim_parallel_gpu(
         i_row = i >> 5
         index = i & 31
 
-        s1 = seed_value + 1 + i_row
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
+        seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
         r = 0.
-        result = ti.f32(0.)
         
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         i_col = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
         while i_col < num_col:
             r += vector[i_col]
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             i_col += uniform_int_distribution(result, 1, clen_value) * 32
         out[i_row] += r * weight_value
 
@@ -390,18 +375,14 @@ def _mv_prob_uniform_cpu(
 
     # ti.loop_config(serialize=True)
     for i_col in range(num_col):
-        s1 = seed_value + 1 + i_col
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
-        result = ti.f32(0.)
+        seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
 
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         i_row = uniform_int_distribution(result, 1, clen_value)
         while i_row < num_row:
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             out[i_row] += uniform_real_distribution(result, w_min_value, w_max_value) * vector[i_col]
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             i_row += uniform_int_distribution(result, 1, clen_value)
 
 @ti.kernel
@@ -426,18 +407,14 @@ def _mv_prob_uniform_gpu(
         i_col = i >> 5
         index = i & 31
         
-        s1 = seed_value + 1 + i_col
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
-        result = ti.f32(0.)
+        seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
         
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         i_row = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
         while i_row < num_row:
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             out[i_row] += uniform_real_distribution(result, w_min_value, w_max_value) * vector[i_col]
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             i_row += uniform_int_distribution(result, 1, clen_value) * 32
 
 @ti.kernel
@@ -459,19 +436,15 @@ def _mv_prob_uniform_outdim_parallel_cpu(
     
     # ti.loop_config(serialize=True)
     for i_row in range(num_row):
-        s1 = seed_value + 1 + i_row
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
-        result = ti.f32(0.)
+        seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
         r = 0.
 
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         i_col = uniform_int_distribution(result, 1, clen_value)
         while i_col < num_col:
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             r += uniform_real_distribution(result, w_min_value, w_max_value) * vector[i_col]
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             i_col += uniform_int_distribution(result, 1, clen_value)
         out[i_row] = r
 
@@ -497,19 +470,15 @@ def _mv_prob_uniform_outdim_parallel_gpu(
         i_row = i >> 5
         index = i & 31
         
-        s1 = seed_value + 1 + i_row
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
-        result = ti.f32(0.)
+        seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
         r = 0.
         
-        s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+        seeds, result = random_generator(seeds)
         i_col = uniform_int_distribution(result, 1, clen_value) + avg_num_uniform * index
         while i_col < num_col:
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             r += uniform_real_distribution(result, w_min_value, w_max_value) * vector[i_col]
-            s1, s2, s3, b, result = random_generator(s1, s2, s3, b)
+            seeds, result = random_generator(seeds)
             i_col += uniform_int_distribution(result, 1, clen_value) * 32
         out[i_row] += r
 
@@ -732,21 +701,16 @@ def _mv_prob_normal_cpu(
     
     # ti.loop_config(serialize=True)
     for i_col in range(num_col):
-        s1 = seed_value + 1 + i_col
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
+        seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
         r = 0.
-        result1 = ti.f32(0.)
-        result2 = ti.f32(0.)
 
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
         i_row = uniform_int_distribution(result1, 1, clen_value)
         while i_row < num_row:
-            s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-            s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+            seeds, result1 = random_generator(seeds)
+            seeds, result2 = random_generator(seeds)
             out[i_row] += normal_distribution(result1, result2, w_mu_value, w_sigma_value) * vector[i_col]
-            s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+            seeds, result1 = random_generator(seeds)
             i_row += uniform_int_distribution(result1, 1, clen_value)
 
 @ti.kernel
@@ -771,21 +735,16 @@ def _mv_prob_normal_gpu(
         i_col = i >> 5
         index = i & 31
         
-        s1 = seed_value + 1 + i_col
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
+        seeds = ti.math.uvec4(seed_value + 1 + i_col, seed_value + 7, seed_value + 15, ti.u32(0))
         r = 0.
-        result1 = ti.f32(0.)
-        result2 = ti.f32(0.)
 
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
         i_row = uniform_int_distribution(result1, 1, clen_value) + avg_num_uniform * index
         while i_row < num_row:
-            s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-            s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+            seeds, result1 = random_generator(seeds)
+            seeds, result2 = random_generator(seeds)
             out[i_row] += normal_distribution(result1, result2, w_mu_value, w_sigma_value) * vector[i_col]
-            s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+            seeds, result1 = random_generator(seeds)
             i_row += uniform_int_distribution(result1, 1, clen_value) * 32
 
 @ti.kernel
@@ -807,21 +766,16 @@ def _mv_prob_normal_outdim_parallel_cpu(
     
     # ti.loop_config(serialize=True)
     for i_row in range(num_row):
-        s1 = seed_value + 1 + i_row
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
+        seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
         r = 0.
-        result1 = ti.f32(0.)
-        result2 = ti.f32(0.)
 
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
         i_col = uniform_int_distribution(result1, 1, clen_value)
         while i_col < num_col:
-            s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-            s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+            seeds, result1 = random_generator(seeds)
+            seeds, result2 = random_generator(seeds)
             r += normal_distribution(result1, result2, w_mu_value, w_sigma_value) * vector[i_col]
-            s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+            seeds, result1 = random_generator(seeds)
             i_col += uniform_int_distribution(result1, 1, clen_value)
         out[i_row] = r
 
@@ -847,21 +801,16 @@ def _mv_prob_normal_outdim_parallel_gpu(
         i_row = i >> 5
         index = i & 31
         
-        s1 = seed_value + 1 + i_row
-        s2 = seed_value + 7
-        s3 = seed_value + 15
-        b = ti.u32(0)
+        seeds = ti.math.uvec4(seed_value + 1 + i_row, seed_value + 7, seed_value + 15, ti.u32(0))
         r = 0.
-        result1 = ti.f32(0.)
-        result2 = ti.f32(0.)
 
-        s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+        seeds, result1 = random_generator(seeds)
         i_col = uniform_int_distribution(result1, 1, clen_value) + avg_num_uniform * index
         while i_col < num_col:
-            s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
-            s1, s2, s3, b, result2 = random_generator(s1, s2, s3, b)
+            seeds, result1 = random_generator(seeds)
+            seeds, result2 = random_generator(seeds)
             r += normal_distribution(result1, result2, w_mu_value, w_sigma_value) * vector[i_col]
-            s1, s2, s3, b, result1 = random_generator(s1, s2, s3, b)
+            seeds, result1 = random_generator(seeds)
             i_col += uniform_int_distribution(result1, 1, clen_value) * 32
         out[i_row] += r
 
diff --git a/brainpy/_src/math/jitconn/tests/test_event_matvec_taichi.py b/brainpy/_src/math/jitconn/tests/test_event_matvec_taichi.py
index 45e2e56b7..8d8c79868 100644
--- a/brainpy/_src/math/jitconn/tests/test_event_matvec_taichi.py
+++ b/brainpy/_src/math/jitconn/tests/test_event_matvec_taichi.py
@@ -15,9 +15,9 @@
   pytest.skip('brainpy.math package may need manual tests.', allow_module_level=True)
 
 shapes = [(100, 200),
-          (10, 1000),
+          # (10, 1000),
           (2, 1000),
-          (1000, 10),
+          # (1000, 10),
           (1000, 2)]
 
 
diff --git a/brainpy/_src/math/taichi_random.py b/brainpy/_src/math/taichi_random.py
index fd27996d0..63212f685 100644
--- a/brainpy/_src/math/taichi_random.py
+++ b/brainpy/_src/math/taichi_random.py
@@ -6,10 +6,11 @@
   # taichi.func random generator implementation
   'taichi_lcg_rand', 'taichi_uniform_int_distribution',
   'taichi_uniform_real_distribution', 'taichi_normal_distribution',
-  'taichi_lfsr88',
+  'taichi_lfsr88', 'RandomGenerator',
 ]
 
 
+
 @ti.func
 def _lcg_rand(state: ti.types.ndarray(ndim=1)):
   # LCG constants
@@ -33,7 +34,7 @@ def taichi_lcg_rand(seed: ti.types.ndarray(ndim=1)):
 
 
 @ti.func
-def taichi_lfsr88(s1: ti.u32, s2: ti.u32, s3: ti.u32, b: ti.u32):
+def taichi_lfsr88(seeds):
   """
   32-bits Random number generator U[0,1): lfsr88
   Author: Pierre L'Ecuyer,
@@ -55,15 +56,14 @@ def taichi_lfsr88(s1: ti.u32, s2: ti.u32, s3: ti.u32, b: ti.u32):
   return ((s1 ^ s2 ^ s3) * 2.3283064365386963e-10);
   }
   ```
-
   """
-  b = (((s1 << 13) ^ s1) >> 19);
-  s1 = (((s1 & ti.u32(4294967294)) << 12) ^ b)
-  b = (((s2 << 2) ^ s2) >> 25)
-  s2 = (((s2 & ti.u32(4294967288)) << 4) ^ b)
-  b = (((s3 << 3) ^ s3) >> 11)
-  s3 = (((s3 & ti.u32(4294967280)) << 17) ^ b)
-  return s1, s2, s3, b, ((s1 ^ s2 ^ s3) * ti.f32(2.3283064365386963e-10))
+  b = ti.cast((((seeds[0] << 13) ^ seeds[0]) >> 19), ti.u32);
+  s1 = (((seeds[0] & ti.u32(4294967294)) << 12) ^ b)
+  b = (((seeds[1] << 2) ^ seeds[1]) >> 25)
+  s2 = (((seeds[1] & ti.u32(4294967288)) << 4) ^ b)
+  b = (((seeds[2] << 3) ^ seeds[2]) >> 11)
+  s3 = (((seeds[2] & ti.u32(4294967280)) << 17) ^ b)
+  return ti.math.uvec4(s1, s2, s3, b), ((s1 ^ s2 ^ s3) * ti.f32(2.3283064365386963e-10))
 
 
 @ti.func
diff --git a/brainpy/math/__init__.py b/brainpy/math/__init__.py
index e24d30ae0..91e3a8260 100644
--- a/brainpy/math/__init__.py
+++ b/brainpy/math/__init__.py
@@ -31,6 +31,7 @@
 from . import fft
 from . import linalg
 from . import random
+from . import taichi_random
 
 # others
 from . import sharding
diff --git a/brainpy/math/taichi_random.py b/brainpy/math/taichi_random.py
new file mode 100644
index 000000000..1f006b952
--- /dev/null
+++ b/brainpy/math/taichi_random.py
@@ -0,0 +1,10 @@
+# -*- coding: utf-8 -*-
+
+from brainpy._src.math.taichi_random import (
+    taichi_lcg_rand as taichi_lcg_rand,
+    taichi_lfsr88 as taichi_lfsr88,
+    taichi_uniform_int_distribution as taichi_uniform_int_distribution,
+    taichi_uniform_real_distribution as taichi_uniform_real_distribution,
+    taichi_normal_distribution as taichi_normal_distribution,
+    RandomGenerator as RandomGenerator,
+)
\ No newline at end of file