diff --git a/.github/workflows/acceptancetest_mg7.yml b/.github/workflows/acceptancetest_mg7.yml
index bc334e43b..3366206ec 100644
--- a/.github/workflows/acceptancetest_mg7.yml
+++ b/.github/workflows/acceptancetest_mg7.yml
@@ -201,6 +201,92 @@ jobs:
             export PATH="$HOME/.cache/HEPtools/bin:$PATH"
             ./tests/test_manager.py test_group_subprocess_mg7 -pA -t0 -l INFO
 
+  acceptancetest_mg7_vegas_reproducibility:
+    needs: build_madspace
+    # mg7 seeded-generation reproducibility (p p > t t~, plain VEGAS-optimized
+    # run): the same run_card seed must give a byte-identical LHE file across
+    # independent runs, a different seed must not. Self-skips if the madspace
+    # + LHAPDF(NNPDF23) stack is unavailable.
+    runs-on: ubuntu-24.04
+    if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true
+    steps:
+      - uses: actions/checkout@v4
+      - uses: ./.github/actions/checkout_mg5
+      - uses: ./.github/actions/install_madspace
+      - uses: ./.github/actions/restore-pip-cache
+      - uses: ./.github/actions/restore_heptools
+      - name: test one of the test test_vegas_reproducibility_mg7
+        run: |
+            cd $GITHUB_WORKSPACE
+            export PATH="$HOME/.cache/HEPtools/bin:$PATH"
+            ./tests/test_manager.py test_vegas_reproducibility_mg7 -pA -t0 -l INFO
+
+  acceptancetest_mg7_madnis_reproducible_mode:
+    needs: build_madspace
+    # mg7 madnis.reproducible = True (p p > t t~): madnis training itself must be
+    # byte-identical for the same seed independent of the cpu thread pool size,
+    # both with online-only training and with buffered (off-policy replay)
+    # training enabled. Self-skips if the madspace + LHAPDF(NNPDF23) stack is
+    # unavailable.
+    runs-on: ubuntu-24.04
+    if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true
+    steps:
+      - uses: actions/checkout@v4
+      - uses: ./.github/actions/checkout_mg5
+      - uses: ./.github/actions/install_madspace
+      - uses: ./.github/actions/restore-pip-cache
+      - uses: ./.github/actions/restore_heptools
+      - name: test one of the test test_madnis_reproducible_mode_mg7
+        run: |
+            cd $GITHUB_WORKSPACE
+            export PATH="$HOME/.cache/HEPtools/bin:$PATH"
+            ./tests/test_manager.py test_madnis_reproducible_mode_mg7 -pA -t0 -l INFO
+
+  acceptancetest_mg7_gridpack_reproducibility:
+    needs: build_madspace
+    # mg7 seeded-generation reproducibility for a gridpack trained with
+    # madnis (p p > t t~): the gridpack's own bin/generate_events --seed must
+    # give a byte-identical LHE file across independent runs with the same
+    # seed, and a different result for a different seed. Training itself is
+    # not required to be deterministic. Self-skips if the madspace +
+    # LHAPDF(NNPDF23) stack is unavailable.
+    runs-on: ubuntu-24.04
+    if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true
+    steps:
+      - uses: actions/checkout@v4
+      - uses: ./.github/actions/checkout_mg5
+      - uses: ./.github/actions/install_madspace
+      - uses: ./.github/actions/restore-pip-cache
+      - uses: ./.github/actions/restore_heptools
+      - name: test one of the test test_gridpack_reproducibility_mg7
+        run: |
+            cd $GITHUB_WORKSPACE
+            export PATH="$HOME/.cache/HEPtools/bin:$PATH"
+            ./tests/test_manager.py test_gridpack_reproducibility_mg7 -pA -t0 -l INFO
+
+  acceptancetest_mg7_gridpack_reproducibility_vegas:
+    needs: build_madspace
+    # mg7 seeded-generation reproducibility for a plain VEGAS-optimized
+    # gridpack (p p > t t~, no madnis training): the gridpack's own
+    # bin/generate_events --seed must give a byte-identical LHE file across
+    # independent runs with the same seed, and a different result for a
+    # different seed. Companion to acceptancetest_mg7_gridpack_reproducibility,
+    # which covers the madnis-trained case. Self-skips if the madspace +
+    # LHAPDF(NNPDF23) stack is unavailable.
+    runs-on: ubuntu-24.04
+    if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true
+    steps:
+      - uses: actions/checkout@v4
+      - uses: ./.github/actions/checkout_mg5
+      - uses: ./.github/actions/install_madspace
+      - uses: ./.github/actions/restore-pip-cache
+      - uses: ./.github/actions/restore_heptools
+      - name: test one of the test test_gridpack_reproducibility_vegas_mg7
+        run: |
+            cd $GITHUB_WORKSPACE
+            export PATH="$HOME/.cache/HEPtools/bin:$PATH"
+            ./tests/test_manager.py test_gridpack_reproducibility_vegas_mg7 -pA -t0 -l INFO
+
   acceptancetest_mg7_merged_flavor_uq:
     needs: build_madspace
     # mg7 cross-section for the merged-flavor u q > u q (q = u d), pinned to the
diff --git a/madgraph/iolibs/template_files/mg7/gridpack.py b/madgraph/iolibs/template_files/mg7/gridpack.py
index 62626f327..ccc71999b 100644
--- a/madgraph/iolibs/template_files/mg7/gridpack.py
+++ b/madgraph/iolibs/template_files/mg7/gridpack.py
@@ -31,7 +31,6 @@
 import tomllib
 import argparse
 
-
 def resolve_verbosity(verbosity: str) -> str:
     """Resolve the run_card "auto" verbosity to "pretty"/"log" depending on
     whether stdout is attached to a terminal; other values pass through
@@ -69,6 +68,12 @@ def main() -> None:
     # parse command line arguments
     parser = argparse.ArgumentParser()
     parser.add_argument("--run_name", type=str, default=run_args["run_name"])
+    parser.add_argument(
+        "--seed", type=int, default=run_args.get("seed", 0),
+        help="0 draws a fresh random seed; a positive integer makes the run "
+             "reproducible at any --cpu_thread_pool_size (for --output_format lhe "
+             "also set combine_thread_pool_size = 1)"
+    )
     parser.add_argument("--device", type=str, nargs="*")
     parser.add_argument(
         "--cpu_thread_pool_size", type=int, default=run_args["cpu_thread_pool_size"]
@@ -173,6 +178,7 @@ def main() -> None:
         channels=channel_generators,
         status_file=os.path.join(run_path, "info.json"),
         config=config,
+        seed=args.seed,
     )
 
     # run generation
diff --git a/madgraph/iolibs/template_files/mg7/madevent.py b/madgraph/iolibs/template_files/mg7/madevent.py
index 6114a4871..002013bae 100644
--- a/madgraph/iolibs/template_files/mg7/madevent.py
+++ b/madgraph/iolibs/template_files/mg7/madevent.py
@@ -212,7 +212,7 @@ def init_context(self) -> None:
         self.device_types = []
         self.devices = []
         self.pool_sizes = []
-        for i, device_name in enumerate(device_names):
+        for device_name in device_names:
             if ":" in device_name:
                 device_type, device_index_str = device_name.split(":")
                 device_index = int(device_index_str)
@@ -426,6 +426,9 @@ def build_event_generator(self, phasespaces: list[PhaseSpace]) -> ms.EventGenera
             channels=channel_generators,
             status_file=os.path.join(self.run_path, "info.json"),
             config=self.event_generator_config,
+            # 0 conventionally means "no seed" in the run card; EventGenerator takes
+            # an actual seed or None (non-deterministic), not a 0 sentinel.
+            seed=self.run_card["run"]["seed"] or None,
         )
         unused_globals = (
             set(self.contexts[0].global_names()) - event_generator.used_globals()
@@ -436,36 +439,42 @@ def build_event_generator(self, phasespaces: list[PhaseSpace]) -> ms.EventGenera
         return event_generator
 
     def survey_phasespaces(
-        self, phasespaces: list[PhaseSpace | None]
+        self, phasespaces: list[PhaseSpace | None], survey_pass: int = 0
     ) -> ms.EventGenerator | None:
         ps_filtered = [ps for ps in phasespaces if ps is not None]
         if len(ps_filtered) == 0:
             return None
         event_generator = self.build_event_generator(ps_filtered)
-        event_generator.survey()
+        event_generator.survey(survey_pass)
         return event_generator
 
     def survey(self) -> None:
+        # survey_pass distinguishes the survey() calls below: "both" mode can
+        # re-survey a channel carried over unchanged from the multichannel pass
+        # into the final (simplified) pass, and both passes schedule jobs on the
+        # same underlying ChannelEventGenerator. The explicit pass index keeps each
+        # pass's job seeds independent of the other passes' job counts, rather than
+        # depending on call history.
         phasespace_mode = self.run_card["phasespace"]["mode"]
         if phasespace_mode == "multichannel":
             self.phasespaces = [
                 subproc.build_multichannel_phasespace()
                 for subproc in self.subprocesses
             ]
-            self.event_generator = self.survey_phasespaces(self.phasespaces)
+            self.event_generator = self.survey_phasespaces(self.phasespaces, 0)
         elif phasespace_mode == "flat":
             self.phasespaces = [
                 subproc.build_flat_phasespace()
                 for subproc in self.subprocesses
             ]
-            self.event_generator = self.survey_phasespaces(self.phasespaces)
+            self.event_generator = self.survey_phasespaces(self.phasespaces, 0)
         elif phasespace_mode == "both":
             kept_count = self.run_card["phasespace"]["simplified_channel_count"]
             phasespaces_multi = [
                 subproc.build_multichannel_phasespace()
                 for subproc in self.subprocesses
             ]
-            evgen_multi = self.survey_phasespaces(phasespaces_multi)
+            evgen_multi = self.survey_phasespaces(phasespaces_multi, 0)
 
             phasespaces_flat = [
                 subproc.build_flat_phasespace()
@@ -473,7 +482,7 @@ def survey(self) -> None:
                 None
                 for subproc in self.subprocesses
             ]
-            #evgen_flat = self.survey_phasespaces(phasespaces_flat, "flat")
+            #evgen_flat = self.survey_phasespaces(phasespaces_flat, 1)
 
             channel_status = evgen_multi.channel_status()
             cross_sections = []
@@ -494,7 +503,7 @@ def survey(self) -> None:
                     self.subprocesses, phasespaces_multi, phasespaces_flat, cross_sections
                 )
             ]
-            self.event_generator = self.survey_phasespaces(self.phasespaces)
+            self.event_generator = self.survey_phasespaces(self.phasespaces, 2)
         else:
             raise ValueError("Unknown phasespace mode")
 
@@ -534,6 +543,7 @@ def train_madnis(self) -> None:
         config.buffer_unweighting_quantile = madnis_args["buffer_unweighting_quantile"]
         config.fixed_cwnet_fraction = madnis_args["fixed_cwnet_fraction"]
         config.softclip_threshold = madnis_args["softclip_threshold"]
+        config.reproducible = madnis_args["reproducible"]
         madnis_phasespaces = []
         integrands = []
         cwnets = []
@@ -559,6 +569,11 @@ def train_madnis(self) -> None:
             config=config,
             integrands=integrands,
             cwnets=cwnets,
+            # Reuses the run_card's own seed (also used by build_event_generator()).
+            # Only the single-channel CPU sample-generation path is currently seeded
+            # -- buffered training and GPU multi-channel batches are still
+            # non-deterministic.
+            seed=run_args["seed"] or None,
         )
         madnis_training.train()
         for phasespace, active_channels in zip(
@@ -1301,6 +1316,9 @@ def simplify_phasespace(
 
     def build_madnis(self, phasespace: PhaseSpace) -> PhaseSpace:
         madnis_args = self.process.run_card["madnis"]
+        # Shared across all networks below: initialize_globals() derives an
+        # independent, non-colliding stream per tensor from this one base seed.
+        seed = self.process.run_card["run"]["seed"] or None
         channels = []
         for channel_id, channel in enumerate(phasespace.channels):
             prefix = f"subproc{self.subproc_id}.channel{channel_id}"
@@ -1318,7 +1336,7 @@ def build_madnis(self, phasespace: PhaseSpace) -> PhaseSpace:
                     subnet_layers=madnis_args["discrete_layers"],
                     subnet_activation=self.activation(madnis_args["discrete_activation"]),
                 )
-                discrete_before.initialize_globals(self.process.contexts[0])
+                discrete_before.initialize_globals(self.process.contexts[0], seed)
                 cond_dim += perm_count
 
             flow_dim = channel.phasespace_mapping.random_dim()
@@ -1333,10 +1351,11 @@ def build_madnis(self, phasespace: PhaseSpace) -> PhaseSpace:
                 invert_spline=madnis_args["flow_invert_spline"],
             )
             if channel.adaptive_mapping is None:
-                flow.initialize_globals(self.process.contexts[0])
+                flow.initialize_globals(self.process.contexts[0], seed)
             else:
                 flow.initialize_from_vegas(
-                    self.process.contexts[0], channel.adaptive_mapping.grid_name()
+                    self.process.contexts[0], channel.adaptive_mapping.grid_name(),
+                    seed
                 )
             cond_dim += flow_dim
 
@@ -1351,7 +1370,7 @@ def build_madnis(self, phasespace: PhaseSpace) -> PhaseSpace:
                     subnet_layers=madnis_args["discrete_layers"],
                     subnet_activation=self.activation(madnis_args["discrete_activation"]),
                 )
-                discrete_after.initialize_globals(self.process.contexts[0])
+                discrete_after.initialize_globals(self.process.contexts[0], seed)
 
             channels.append(Channel(
                 phasespace_mapping = channel.phasespace_mapping,
@@ -1420,7 +1439,9 @@ def build_cwnet(self, channel_count: int) -> ms.ChannelWeightNetwork:
             activation=self.activation(madnis_args["cwnet_activation"]),
             prefix=f"subproc{self.subproc_id}.cwnet",
         )
-        cwnet.initialize_globals(self.process.contexts[0])
+        cwnet.initialize_globals(
+            self.process.contexts[0], self.process.run_card["run"]["seed"] or None
+        )
         return cwnet
 
     def t_channel_mode(self, name: str) -> ms.PhaseSpaceMapping.TChannelMode:
diff --git a/madgraph/iolibs/template_files/mg7/run_card.toml b/madgraph/iolibs/template_files/mg7/run_card.toml
index 28410c262..35b9bc173 100644
--- a/madgraph/iolibs/template_files/mg7/run_card.toml
+++ b/madgraph/iolibs/template_files/mg7/run_card.toml
@@ -1,5 +1,11 @@
 [run]
 run_name = %(run.run_name)s
+# 0 draws a fresh random seed each run; a positive integer makes the run
+# reproducible: the same (seed, thread-pool sizes) reproduces bit-identically and
+# distinct seeds are statistically independent. Generation is reproducible at any
+# cpu_thread_pool_size; for output_format = lhe also set combine_thread_pool_size = 1
+# (compact_npy and lhe_npy are reproducible at any combine size).
+seed = %(run.seed)s
 # options: cuda, hip, cpp, cppnone, cppsse4, cppavx2, cpp512y, cpp512z, cppauto
 devices = %(run.devices)s
 # options:
@@ -136,3 +142,4 @@ batch_size_threshold = %(madnis.batch_size_threshold)s
 channel_grouping_mode = %(madnis.channel_grouping_mode)s # options: none, uniform, learned
 fixed_cwnet_fraction = %(madnis.fixed_cwnet_fraction)s
 softclip_threshold = %(madnis.softclip_threshold)s
+reproducible = %(madnis.reproducible)s # make training fully reproducible for a given seed, independent of thread count (CPU only, some throughput/memory cost)
diff --git a/madgraph/various/banner.py b/madgraph/various/banner.py
index 5d5b0c63d..474d6465e 100755
--- a/madgraph/various/banner.py
+++ b/madgraph/various/banner.py
@@ -6453,6 +6453,12 @@ def default_setup(self):
 
         # ----------------------------- [run] --------------------------
         self.add_toml_param('run', 'run_name', "run", gridpack=True)
+        self.add_toml_param('run', 'seed', 0, gridpack=True,
+            comment="0 draws a fresh random seed each run; any positive integer "
+                    "makes the run reproducible: the same (seed, thread-pool sizes) "
+                    "reproduces bit-identically. Generation is reproducible at any "
+                    "cpu_thread_pool_size; for output_format = lhe also set "
+                    "combine_thread_pool_size = 1")
         self.add_toml_param('run', 'devices', ["cppnone"], typelist=str, gridpack=True,
             comment="options: cuda, hip, cpp, cppnone, cppsse4, cppavx2, cpp512y, cpp512z, cppauto")
         self.add_toml_param('run', 'simd_vector_size', -1,
@@ -6588,6 +6594,7 @@ def default_setup(self):
             allowed=['none', 'uniform', 'learned'])
         self.add_toml_param('madnis', 'fixed_cwnet_fraction', 0.33)
         self.add_toml_param('madnis', 'softclip_threshold', 30.0)
+        self.add_toml_param('madnis', 'reproducible', False)
 
         # ----------------- dynamic (free-form) sections ---------------
         self.dynamic_sections['multiparticles'] = collections.OrderedDict([
diff --git a/madspace/.clang-format-ignore b/madspace/.clang-format-ignore
index 2ff5b5bfb..9343bc477 100644
--- a/madspace/.clang-format-ignore
+++ b/madspace/.clang-format-ignore
@@ -7,3 +7,5 @@ src/cpu/runtime_mixin.inc
 src/gpu/runtime_mixin.inc
 src/cpu/runtime_backward_mixin.inc
 src/gpu/runtime_backward_mixin.inc
+include/madspace/mixmax/mixmax.hpp
+include/madspace/mixmax/mixmax_skip_N17.c
diff --git a/madspace/CMakeLists.txt b/madspace/CMakeLists.txt
index c7a6dd13d..38f3382c0 100644
--- a/madspace/CMakeLists.txt
+++ b/madspace/CMakeLists.txt
@@ -180,6 +180,7 @@ add_library(
     src/driver/backend.cpp
     src/driver/lhe_output.cpp
     src/driver/madnis_training.cpp
+    src/driver/random.cpp
     src/compgraphs/type.cpp
     src/compgraphs/function.cpp
     src/compgraphs/instruction.cpp
@@ -232,6 +233,7 @@ add_library(
     include/madspace/driver/lhe_output.hpp
     include/madspace/driver/logger.hpp
     include/madspace/driver/madnis_training.hpp
+    include/madspace/driver/random.hpp
     include/madspace/compgraphs.hpp
     include/madspace/compgraphs/type.hpp
     include/madspace/compgraphs/function.hpp
@@ -270,6 +272,8 @@ add_library(
     include/madspace/phasespace/matrix_element.hpp
     include/madspace/phasespace/cross_section.hpp
     include/madspace/phasespace/scale.hpp
+    include/madspace/mixmax/mixmax.hpp
+    include/madspace/mixmax/mixmax_skip_N17.c
 )
 
 #target_compile_options(madspace PRIVATE -Wall -Wextra -Wpedantic -Werror)
diff --git a/madspace/include/madspace/compgraphs/instruction.hpp b/madspace/include/madspace/compgraphs/instruction.hpp
index 184c83948..0b4125483 100644
--- a/madspace/include/madspace/compgraphs/instruction.hpp
+++ b/madspace/include/madspace/compgraphs/instruction.hpp
@@ -28,6 +28,9 @@ class Instruction {
     const std::string& name() const { return _name; }
     int opcode() const { return _opcode; }
     bool differentiable() const { return _differentiable; }
+    // draws from a shared RNG stream: never CSE'd, always scheduled on the main
+    // GPU stream (see FunctionBuilder::instruction)
+    virtual bool is_random() const { return false; }
 
 protected:
     void check_arg_count(const ValueVec& args, std::size_t count) const;
@@ -197,6 +200,7 @@ class RandomInstruction : public Instruction {
     RandomInstruction(int opcode, bool differentiable) :
         Instruction("random", opcode, differentiable) {}
     TypeVec signature(const ValueVec& args) const override;
+    bool is_random() const override { return true; }
 };
 
 class RandomIntInstruction : public Instruction {
@@ -204,6 +208,7 @@ class RandomIntInstruction : public Instruction {
     RandomIntInstruction(int opcode, bool differentiable) :
         Instruction("random_int", opcode, differentiable) {}
     TypeVec signature(const ValueVec& args) const override;
+    bool is_random() const override { return true; }
 };
 
 class UnweightInstruction : public Instruction {
@@ -211,6 +216,7 @@ class UnweightInstruction : public Instruction {
     UnweightInstruction(int opcode, bool differentiable) :
         Instruction("unweight", opcode, differentiable) {}
     TypeVec signature(const ValueVec& args) const override;
+    bool is_random() const override { return true; }
 };
 
 class MatrixElementInstruction : public Instruction {
diff --git a/madspace/include/madspace/driver/backend.hpp b/madspace/include/madspace/driver/backend.hpp
index 792d5828b..7dbbeeada 100644
--- a/madspace/include/madspace/driver/backend.hpp
+++ b/madspace/include/madspace/driver/backend.hpp
@@ -1,7 +1,11 @@
 #pragma once
 
+#include <cstdint>
+#include <optional>
+
 #include "madspace/compgraphs.hpp"
 #include "madspace/driver/context.hpp"
+#include "madspace/driver/random.hpp"
 #include "madspace/driver/tensor.hpp"
 
 namespace madspace {
@@ -9,6 +13,8 @@ namespace madspace {
 class Runtime {
 public:
     virtual ~Runtime() = default;
+    // `seed`, when set, makes this call's random numbers a deterministic function of
+    // `seed` alone. Backends that can't support this should reject it, not ignore it.
     virtual TensorVec run(const TensorVec& inputs) = 0;
     virtual std::tuple<TensorVec, TensorVec, std::vector<bool>> run_with_grad(
         const TensorVec& inputs, const std::vector<bool>& input_requires_grad
@@ -19,6 +25,8 @@ class Runtime {
         const std::vector<bool>& eval_grad,
         bool return_contiguous_grads = false
     ) = 0;
+    virtual void set_seed(DerivedSeed seed) = 0;
+
     friend std::unique_ptr<Runtime>
     build_runtime(const Function& function, ContextPtr context, bool concurrent);
 
diff --git a/madspace/include/madspace/driver/channel_generator.hpp b/madspace/include/madspace/driver/channel_generator.hpp
index 957037d15..8b6005d95 100644
--- a/madspace/include/madspace/driver/channel_generator.hpp
+++ b/madspace/include/madspace/driver/channel_generator.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include <random>
+#include <optional>
 #include <unordered_set>
 #include <vector>
 
@@ -13,6 +13,7 @@
 #include "madspace/driver/discrete_optimizer.hpp"
 #include "madspace/driver/generator_data.hpp"
 #include "madspace/driver/io.hpp"
+#include "madspace/driver/random.hpp"
 #include "madspace/driver/vegas_optimizer.hpp"
 #include "madspace/phasespace.hpp"
 
@@ -49,7 +50,9 @@ class ChannelEventGenerator {
     bool needs_optimization() const {
         return (_vegas_optimizer || _discrete_optimizer) && !_status.optimized;
     }
-    void set_target_count(double target_count) { _status.count_target = target_count; }
+    void set_target_count(std::size_t target_count) {
+        _status.count_target = target_count;
+    }
     const std::unordered_set<std::string>& used_globals() const {
         return _used_globals;
     }
@@ -57,11 +60,24 @@ class ChannelEventGenerator {
     int particle_layout_extra_flags() const { return _particle_layout_extra_flags; }
     const DataLayout& event_file_layout() const { return _event_file_layout; }
 
-    void unweight_file(std::mt19937& rand_gen);
+    void unweight_file(MixMaxRandom& rand_gen);
     void integrate(const GeneratorBatchJob& job);
     void optimize_vegas(const GeneratorBatchJob& job);
     double channel_weight_sum(std::size_t event_count);
-    void start_job(GeneratorBatchJob& job, ResultQueue& result_queue);
+    void start_job(
+        GeneratorBatchJob& job,
+        ResultQueue& result_queue,
+        std::optional<std::uint64_t> seed,
+        bool is_survey,
+        std::size_t survey_pass
+    );
+    // Snapshots max_weight into job.max_weight, at a fixed point independent of
+    // when submit_unweight_job() actually dispatches it.
+    void prepare_unweight_job(GeneratorBatchJob& job) const;
+    // Submits job's unweighting to the thread pool of its own generation context
+    // (required on GPU: unweighting does a device-to-host copy).
+    void submit_unweight_job(GeneratorBatchJob& job, ResultQueue& result_queue);
+    // prepare_unweight_job() + submit_unweight_job() combined.
     void start_unweight_job(GeneratorBatchJob& job, ResultQueue& result_queue);
     std::size_t next_vegas_batch_size();
     void clear_events();
@@ -127,7 +143,17 @@ class ChannelEventGenerator {
     std::optional<Function> _histogram_function;
     RunningIntegral _cross_section;
     double _max_weight = 0.;
+    // Monotonic per-job counters keying each job's deterministic random stream;
+    // never reset. Separate for survey/generate so neither depends on the other.
+    std::size_t _survey_rng_seq = 0;
+    std::size_t _generate_rng_seq = 0;
+    // Progress of unweight_file()'s final pass over _weight_file: _unweighted_count
+    // is the index up to which it's been scanned under the *current* _max_weight,
+    // and _unweighted_accept_count the number of those accepted. Both reset to 0
+    // whenever _max_weight changes (forcing a full rescan), since a changed
+    // max_weight invalidates every prior accept/reject decision.
     std::size_t _unweighted_count = 0;
+    std::size_t _unweighted_accept_count = 0;
     std::size_t _iters_without_improvement = 0;
     double _best_rsd = std::numeric_limits<double>::max();
     std::vector<double> _large_weights;
diff --git a/madspace/include/madspace/driver/context.hpp b/madspace/include/madspace/driver/context.hpp
index 13e438084..7eaa4c91d 100644
--- a/madspace/include/madspace/driver/context.hpp
+++ b/madspace/include/madspace/driver/context.hpp
@@ -48,35 +48,43 @@ class MatrixElementApi {
     std::size_t index() const { return _index; }
     const std::string& file_name() const { return _file_name; }
     std::vector<bool> supported_inputs() const {
-        bool const* data; int count;
+        bool const* data;
+        int count;
         check_umami_status(_supported_inputs(&data, &count));
         std::vector<bool> result(UMAMI_INPUT_KEY_COUNT, false);
-        for (int i = 0; i < count && i < UMAMI_INPUT_KEY_COUNT; ++i)
+        for (int i = 0; i < count && i < UMAMI_INPUT_KEY_COUNT; ++i) {
             result[i] = data[i];
+        }
         return result;
     }
     std::vector<bool> required_inputs() const {
-        bool const* data; int count;
+        bool const* data;
+        int count;
         check_umami_status(_required_inputs(&data, &count));
         std::vector<bool> supported = supported_inputs();
         for (int i = 0; i < count && i < UMAMI_INPUT_KEY_COUNT; ++i) {
             if (data[i] && !supported[i]) {
-                throw_error(std::format(
-                    "input key {} is reported as required but not as supported", i
-                ));
+                throw_error(
+                    std::format(
+                        "input key {} is reported as required but not as supported", i
+                    )
+                );
             }
         }
         std::vector<bool> result(UMAMI_INPUT_KEY_COUNT, false);
-        for (int i = 0; i < count && i < UMAMI_INPUT_KEY_COUNT; ++i)
+        for (int i = 0; i < count && i < UMAMI_INPUT_KEY_COUNT; ++i) {
             result[i] = data[i];
+        }
         return result;
     }
     std::vector<bool> supported_outputs() const {
-        bool const* data; int count;
+        bool const* data;
+        int count;
         check_umami_status(_supported_outputs(&data, &count));
         std::vector<bool> result(UMAMI_OUTPUT_KEY_COUNT, false);
-        for (int i = 0; i < count && i < UMAMI_OUTPUT_KEY_COUNT; ++i)
+        for (int i = 0; i < count && i < UMAMI_OUTPUT_KEY_COUNT; ++i) {
             result[i] = data[i];
+        }
         return result;
     }
 
@@ -170,6 +178,7 @@ class Context {
     void load_globals(const std::string& dir);
     DevicePtr device() { return _device; }
     ThreadPool& thread_pool() { return *_thread_pool; }
+    std::size_t unique_seed_index() { return _seed_index++; }
 
 private:
     DevicePtr _device;
@@ -177,6 +186,7 @@ class Context {
     std::unordered_map<std::string, std::pair<Tensor, bool>> _globals;
     std::vector<std::unique_ptr<MatrixElementApi>> _matrix_elements;
     std::vector<std::string> _param_card_paths;
+    std::size_t _seed_index = 0;
 };
 
 using ContextPtr = std::shared_ptr<Context>;
diff --git a/madspace/include/madspace/driver/event_generator.hpp b/madspace/include/madspace/driver/event_generator.hpp
index 0e9428dff..623c046c7 100644
--- a/madspace/include/madspace/driver/event_generator.hpp
+++ b/madspace/include/madspace/driver/event_generator.hpp
@@ -2,7 +2,7 @@
 
 #include <chrono>
 #include <optional>
-#include <random>
+#include <set>
 #include <vector>
 
 #include <nlohmann/json.hpp>
@@ -31,14 +31,18 @@ class EventGenerator {
         const std::vector<ContextPtr>& contexts,
         const std::vector<std::shared_ptr<ChannelEventGenerator>>& channels,
         const std::string& status_file = "",
-        const GeneratorConfig& config = default_config
+        const GeneratorConfig& config = default_config,
+        std::optional<std::uint64_t> seed = std::nullopt
     );
-    void survey();
+    // `survey_pass` salts job seeds so repeated survey() calls on the same
+    // channel (e.g. re-survey after simplification) don't share a seed stream.
+    void survey(std::size_t survey_pass = 0);
     void generate();
     void combine_to_compact_npy(const std::string& file_name);
     void combine_to_lhe_npy(const std::string& file_name, LHECompleter& lhe_completer);
     void combine_to_lhe(
-        const std::string& file_name, LHECompleter& lhe_completer,
+        const std::string& file_name,
+        LHECompleter& lhe_completer,
         const LHEMeta& meta = {}
     );
     GeneratorStatus status() const { return _status; }
@@ -73,8 +77,42 @@ class EventGenerator {
     std::vector<bool> _channel_optimizing;
     std::vector<double> _channel_integral_fractions;
     std::vector<std::size_t> _context_job_counts;
+    // True while a channel has a steady-state batch dispatched but not yet fully
+    // committed; keeps next_batch_job_count() from double-counting in-flight work.
+    std::vector<bool> _channel_batch_pending;
+    // generate_deterministic() only: per-channel commit cursor/buffer, analogous
+    // to _ready_gen/_commit_cursor but ordered per channel instead of globally.
+    std::vector<std::set<std::size_t>> _channel_ready_gen;
+    std::vector<std::size_t> _channel_commit_cursor;
+    std::vector<bool> _channel_cursor_set;
+    // generate_deterministic() only: same as above, for a job's unweight-stage
+    // completion (tracked separately since it's a distinct completion event).
+    std::vector<std::set<std::size_t>> _channel_unweight_ready;
+    std::vector<std::size_t> _channel_unweight_cursor;
+    // generate_deterministic() only: per-context queue of job ids awaiting
+    // unweight-stage dispatch, drained with priority by start_jobs().
+    std::vector<std::vector<std::size_t>> _context_unweight_queue;
     ResultQueue _result_queue;
 
+    // Base seed for reproducible event generation; nullopt means non-deterministic.
+    std::optional<std::uint64_t> _seed;
+
+    // unweight_all() may run more than once per generate() (a channel's target can
+    // grow after it looked done, un-finishing it and triggering another round).
+    // Salted by this counter so repeated calls don't replay the same stream.
+    std::size_t _unweight_call_index = 0;
+
+    // Scheduling context for the running survey()/generate() call, read by
+    // start_jobs() to derive job seeds.
+    bool _survey_job = false;
+    std::size_t _survey_pass = 0;
+
+    // Deterministic-path state, set from _seed. Generate completions are
+    // committed in ascending job id, with _commit_cursor as the next id due.
+    bool _deterministic = false;
+    std::set<std::size_t> _ready_gen;
+    std::size_t _commit_cursor = 0;
+
     std::chrono::time_point<std::chrono::steady_clock> _start_time;
     std::size_t _start_cpu_microsec;
     std::chrono::time_point<std::chrono::steady_clock> _last_print_time;
@@ -84,8 +122,17 @@ class EventGenerator {
     std::string _status_file;
     std::unordered_map<std::string, TimingData> _timing_data;
 
-    bool start_jobs();
+    void survey_deterministic();
+    void generate_deterministic();
+    void commit_generate_job(GeneratorBatchJob& job);
+    void commit_unweight_job(GeneratorBatchJob& job);
+    void finish_channel_job(const GeneratorBatchJob& job);
+    void register_dispatched_ids(std::size_t first_id, std::size_t end_id);
+    std::size_t next_batch_job_count(std::size_t channel_index) const;
+    std::size_t start_jobs();
     void update_integral();
+    void update_integral_status();
+    void update_integral_fractions();
     void update_counts();
     void reset_start_time();
     void add_timing_data(const std::string& key);
@@ -94,14 +141,15 @@ class EventGenerator {
     void read_and_combine(
         std::vector<CombineChannelData>& channel_data,
         EventBuffer& buffer,
-        double norm_factor
+        double norm_factor,
+        MixMaxRandom& rand_gen
     );
     void fill_lhe_event(
         LHECompleter& lhe_completer,
         LHEEvent& lhe_event,
         EventBuffer& buffer,
         std::size_t event_index,
-        std::mt19937& rand_gen
+        MixMaxRandom& rand_gen
     );
 
     void init_status(const std::string& status);
diff --git a/madspace/include/madspace/driver/generator_data.hpp b/madspace/include/madspace/driver/generator_data.hpp
index bedf6bb2c..17a6c058b 100644
--- a/madspace/include/madspace/driver/generator_data.hpp
+++ b/madspace/include/madspace/driver/generator_data.hpp
@@ -2,6 +2,7 @@
 
 #include <cmath>
 #include <cstddef>
+#include <optional>
 #include <string>
 #include <vector>
 
@@ -63,6 +64,9 @@ struct GeneratorConfig {
     int combine_thread_count = -1;
     double cut_efficiency_threshold = 0.7;
     std::size_t max_cut_repetitions = 100;
+    // Cap/floor on a steady-state batch's job count; see next_batch_job_count().
+    double generation_batch_fraction = 0.5;
+    std::size_t min_batch_jobs = 1;
 };
 
 struct GeneratorStatus {
@@ -76,7 +80,7 @@ struct GeneratorStatus {
     std::size_t count_after_cuts;
     std::size_t count_after_cuts_opt;
     double count_unweighted;
-    double count_target;
+    std::size_t count_target;
     std::size_t iterations;
     bool optimized;
     bool done;
@@ -93,6 +97,7 @@ struct Histogram {
 struct GeneratorBatchJob {
     std::size_t channel_index;
     bool unweight;
+    // Nonzero: start_jobs() splits this into sub-jobs and dispatches them atomically.
     std::size_t vegas_batch_size;
     std::size_t split_job_count;
     Tensor weights;
@@ -104,6 +109,16 @@ struct GeneratorBatchJob {
     std::size_t context_index;
     std::size_t job_id;
     double max_weight;
+    // Top-level seed plus job identity, used to derive this job's DerivedSeed(s)
+    // independently for generate vs unweight (see start_job()/submit_unweight_job()).
+    std::optional<std::uint64_t> rng_seed;
+    // Per-channel dispatch sequence, assigned once at start_job() time.
+    std::size_t rng_job_index = 0;
+    bool rng_is_survey = false;
+    std::size_t rng_survey_pass = 0;
+    // True for a VEGAS-grid-optimization batch, false for a steady-state generation
+    // batch -- both use vegas_batch_size > 0 for atomic whole-batch dispatch.
+    bool is_vegas_batch = false;
 };
 
 void to_json(nlohmann::json& j, const GeneratorStatus& status);
diff --git a/madspace/include/madspace/driver/lhe_output.hpp b/madspace/include/madspace/driver/lhe_output.hpp
index 630175a22..92eaaf42b 100644
--- a/madspace/include/madspace/driver/lhe_output.hpp
+++ b/madspace/include/madspace/driver/lhe_output.hpp
@@ -1,13 +1,13 @@
 #pragma once
 
 #include <fstream>
-#include <random>
 #include <string>
 #include <unordered_map>
 #include <vector>
 
 #include <nlohmann/json.hpp>
 
+#include "madspace/driver/random.hpp"
 #include "madspace/driver/thread_pool.hpp"
 #include "madspace/phasespace/topology.hpp"
 #include "madspace/util.hpp"
@@ -86,7 +86,7 @@ class LHECompleter {
         int color_index,
         int flavor_index,
         int helicity_index,
-        std::mt19937& rand_gen
+        MixMaxRandom& rand_gen
     );
     std::size_t max_particle_count() const { return _max_particle_count; }
     void save(const std::string& file) const;
diff --git a/madspace/include/madspace/driver/madnis_training.hpp b/madspace/include/madspace/driver/madnis_training.hpp
index 932760119..8e8179028 100644
--- a/madspace/include/madspace/driver/madnis_training.hpp
+++ b/madspace/include/madspace/driver/madnis_training.hpp
@@ -1,6 +1,8 @@
 #pragma once
 
 #include <chrono>
+#include <cstdint>
+#include <optional>
 
 #include "madspace/compgraphs.hpp"
 #include "madspace/driver/adam_optimizer.hpp"
@@ -41,13 +43,22 @@ class MadnisTraining {
         double buffer_unweighting_quantile = 0.99;
         double fixed_cwnet_fraction = 0.33;
         double softclip_threshold = 0.0;
+        // Makes generator-job scheduling (CPU only) deterministic given the same
+        // seed, independent of thread count, at some cost to throughput/memory.
+        // GPU multi-channel batches are unaffected and stay non-deterministic.
+        bool reproducible = false;
     };
     MadnisTraining(
         ContextPtr generator_context,
         ContextPtr optimizer_context,
         const Config& config,
         const std::vector<std::shared_ptr<Integrand>>& integrands,
-        const std::optional<ChannelWeightNetwork>& cwnet
+        const std::optional<ChannelWeightNetwork>& cwnet,
+        std::optional<std::uint64_t> seed = std::nullopt,
+        // Offset added to this subprocess's local channel indices when deriving
+        // seeds, so MultiMadnisTraining's subprocesses (which all share the same
+        // top-level seed) get non-overlapping DerivedSeed channel_index ranges.
+        std::size_t channel_index_offset = 0
     );
     void train_step(std::size_t batch_index);
     std::vector<std::size_t> active_channels() const;
@@ -65,6 +76,9 @@ class MadnisTraining {
     struct SampleJob {
         SampleBatch samples;
         SampleBatch unweighted_samples;
+        // dispatch sequence used to commit in order: per-channel for single-channel
+        // jobs, global (see _multi_job_next_dispatch_seq) for multi-channel jobs
+        std::size_t dispatch_seq = 0;
     };
     struct ChannelData {
         std::size_t index;
@@ -77,6 +91,13 @@ class MadnisTraining {
         RuntimePtr generator_runtime = nullptr;
         RuntimePtr unweighter_runtime = nullptr;
         SampleBatch buffer;
+        // commit-ordering state for single-channel generator jobs (see
+        // process_job_results)
+        std::size_t next_dispatch_seq = 0;
+        std::size_t commit_cursor = 0;
+        std::unordered_map<std::size_t, std::size_t> ready_job_ids;
+        // reproducible mode: samples staged here, flushed into buffer next round
+        std::vector<SampleBatch> pending_buffer_samples;
     };
 
     inline static std::function<void(void)> _abort_check_function = [] {};
@@ -84,6 +105,9 @@ class MadnisTraining {
     void build_runtimes_and_optimizer();
     std::vector<std::size_t> compute_channel_sizes();
     void start_generator_jobs(const std::vector<std::size_t>& channel_fractions);
+    void maybe_start_generator_jobs(
+        const std::vector<std::size_t>& channel_fractions, bool is_online_attempt
+    );
     TensorVec permute_tensors(const TensorVec& tensors) const;
     void start_single_job(std::size_t channel_index, std::size_t batch_size);
     void start_multi_job(const std::vector<std::size_t> batch_sizes);
@@ -92,6 +116,7 @@ class MadnisTraining {
     TensorVec build_online_training_batch(const std::vector<size_t>& counts);
     TensorVec build_buffered_training_batch(const std::vector<size_t>& counts);
     void process_job_results(const std::vector<std::size_t>& job_ids);
+    void process_all_jobs();
     void buffer_store(ChannelData& channel, SampleBatch& samples);
     void
     update_history(const TensorVec& results, const std::vector<std::size_t>& counts);
@@ -115,6 +140,15 @@ class MadnisTraining {
     std::vector<std::size_t> _arg_permutation;
     bool _buffer_ready = false;
     std::vector<std::size_t> _active_flavors_count;
+    std::optional<std::uint64_t> _seed;
+    std::size_t _channel_index_offset;
+    // sequence for seeding build_buffered_training_batch's BatchSampler::run() call
+    std::size_t _buffered_batch_seq = 0;
+    // commit-ordering state for multi-channel (GPU) generator jobs (see
+    // process_job_results): global, since one job spans multiple channels at once
+    std::size_t _multi_job_next_dispatch_seq = 0;
+    std::size_t _multi_job_commit_cursor = 0;
+    std::unordered_map<std::size_t, std::size_t> _multi_job_ready_job_ids;
 };
 
 class MultiMadnisTraining {
@@ -124,7 +158,8 @@ class MultiMadnisTraining {
         ContextPtr optimizer_context,
         const MadnisTraining::Config& config,
         const nested_vector2<std::shared_ptr<Integrand>>& integrands,
-        const std::vector<std::optional<ChannelWeightNetwork>>& cwnets
+        const std::vector<std::optional<ChannelWeightNetwork>>& cwnets,
+        std::optional<std::uint64_t> seed = std::nullopt
     );
     void train();
     nested_vector2<std::size_t> active_channels() const;
diff --git a/madspace/include/madspace/driver/random.hpp b/madspace/include/madspace/driver/random.hpp
new file mode 100644
index 000000000..2a4d6d39b
--- /dev/null
+++ b/madspace/include/madspace/driver/random.hpp
@@ -0,0 +1,72 @@
+#pragma once
+
+#include <algorithm>
+#include <array>
+#include <cstdint>
+#include <optional>
+
+#include "madspace/mixmax/mixmax.hpp"
+
+namespace madspace {
+
+struct DerivedSeed {
+    static constexpr std::size_t max_channel_count = 1ULL << 12;
+    static constexpr std::size_t max_job_count = 1ULL << 32;
+    static constexpr std::size_t max_stream_count = 1ULL << 16;
+    enum SeedType {
+        none,
+        first_survey_generate,
+        first_survey_unweight,
+        second_survey_generate,
+        second_survey_unweight,
+        generator_generate,
+        generator_unweight,
+        combine_select,
+        lhe_complete,
+        unweight_pass,
+        madnis_generate,
+        madnis_unweight,
+        madnis_sample_buffer,
+        global_init,
+    };
+
+    std::array<std::uint32_t, 4> seed_parts;
+
+    DerivedSeed(
+        const std::optional<std::uint64_t>& seed = std::nullopt,
+        SeedType seed_type = none,
+        std::size_t job_index = 0,
+        std::size_t channel_index = 0,
+        std::size_t stream_index = 0
+    );
+};
+
+class MixMaxRandom {
+public:
+    MixMaxRandom() : MixMaxRandom(DerivedSeed()) {}
+    MixMaxRandom(DerivedSeed seed) :
+        _mixmax(
+            seed.seed_parts[0],
+            seed.seed_parts[1],
+            seed.seed_parts[2],
+            seed.seed_parts[3]
+        ) {}
+    explicit MixMaxRandom(std::uint64_t seed) : MixMaxRandom(DerivedSeed(seed)) {}
+    void set_seed(DerivedSeed seed) {
+        _mixmax.seed_uniquestream(
+            seed.seed_parts[0],
+            seed.seed_parts[1],
+            seed.seed_parts[2],
+            seed.seed_parts[3]
+        );
+    }
+    double generate_double() { return _mixmax.flat(); }
+    std::size_t generate_int(std::size_t max_int) {
+        return std::min<std::size_t>(_mixmax.flat() * max_int, max_int - 1);
+    }
+
+private:
+    mixmax_engine _mixmax;
+};
+
+} // namespace madspace
diff --git a/madspace/include/madspace/mixmax/mixmax.hpp b/madspace/include/madspace/mixmax/mixmax.hpp
new file mode 100644
index 000000000..5d3a67b31
--- /dev/null
+++ b/madspace/include/madspace/mixmax/mixmax.hpp
@@ -0,0 +1,375 @@
+/*
+ * mixmax.hpp
+ *
+ * C++ implementation of the MIXMAX random number generator.
+ *
+ *  Copyright (2008-2023) by Konstantin Savvidy.
+ *
+ *  Free to use, academic or commercial. Do not redistribute without permission.
+ *
+ *  G.K.Savvidy and N.G.Ter-Arutyunian,
+ *  On the Monte Carlo simulation of physical systems,
+ *  J.Comput.Phy 97, 566 (1991);
+ *  Preprint EPI-865-16-86, Yerevan, Jan. 1986
+ *
+ *  K.Savvidy
+ *  The MIXMAX random number generator
+ *  Comp. Phy Commun. 196 (2015), pp 161-165
+ *  http://dx.doi.org/10.1016/j.cpc.2015.06.003
+ *
+ *  K.Savvidy and G.Savvidy
+ *  Spectrum and Entropy of C-system MIXMAX random number generator
+ *  Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33-38
+ *  http://dx.doi.org/10.1016/j.chao2016.05.003
+ *
+ */
+
+#ifndef __MIXMAX_H
+#define __MIXMAX_H
+
+#include <cstdint>
+
+#if (defined(__CUDACC__) || defined (__CUDA_ARCH__) || defined(__HIPCC__))
+#define _dev __host__ __device__
+#else
+#define _dev
+#include <array>
+#include <iostream>
+#endif
+
+class alignas(128) mixmax_engine
+{
+static const int N = 17;
+static constexpr int BITS=61;
+static constexpr uint64_t M61=2305843009213693951ULL;
+static constexpr uint64_t MERSBASE=M61;
+static constexpr double INV_MERSBASE=(0.43368086899420177360298E-18);
+
+
+    static constexpr long long int SPECIAL   = 0;
+    static constexpr long long int SPECIALMUL= 36;
+    // Note the potential for confusion...
+
+private: // state
+    uint64_t V[N] = {0}; // init with all zero's - not the same as seeding with seed=0!
+    uint64_t sumtot = {0};
+    int counter = N;
+
+public:
+using result_type = uint64_t;                                 // should it be double?
+    static constexpr uint64_t min() {return 0;}
+    static constexpr uint64_t max() {return M61;}
+
+    _dev inline uint64_t get_next() ;
+    _dev inline double flat();
+
+
+    _dev inline mixmax_engine(); // Constructor, no seeds
+    _dev inline mixmax_engine(uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID );       // Constructor with four 32-bit seeds
+    _dev inline mixmax_engine(uint64_t seedval){seed_uniquestream( 0, 0, (uint32_t)(seedval>>32), (uint32_t)seedval );}; // Constructor with one 64bit seed
+    _dev inline void seed(uint64_t seedval){seed_uniquestream( 0, 0, (uint32_t)(seedval>>32), (uint32_t)seedval );} // seed with one 64-bit seed
+    _dev inline void seed_uniquestream( uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID );
+
+    _dev inline mixmax_engine Branch();
+    _dev inline void BranchInplace();
+
+
+    _dev inline mixmax_engine& operator=(const mixmax_engine& other );
+    _dev inline operator double() { return flat(); }
+    _dev inline operator float() { return float( flat() ); }
+    _dev inline operator unsigned int() { return static_cast<unsigned int>(get_next()); }
+
+    _dev inline uint64_t operator()()
+    {
+        return get_next();
+    }
+
+private:
+    _dev static inline uint64_t MOD_MULSPEC(uint64_t k);
+    _dev inline void seed_vielbein(); // seeds with the unit vector {1,0,0,...}
+    _dev static inline uint64_t iterate_raw_vec(uint64_t* Y, uint64_t sumtotOld);
+    _dev inline uint64_t apply_bigskip(uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID );
+    _dev static inline uint64_t modadd(uint64_t foo, uint64_t bar){return MOD_MERSENNE(foo+bar);};
+    _dev static inline uint64_t fmodmulM61(uint64_t cum, uint64_t s, uint64_t a);
+    _dev static inline uint64_t MOD_MERSENNE(uint64_t k) {return ((((k)) & MERSBASE) + (((k)) >> BITS) );}
+    _dev static inline uint64_t MULWU(uint64_t k) {return (( (k)<<(SPECIALMUL) & M61) | ( (k) >> (BITS-SPECIALMUL))  );}
+    void inline print_state();
+
+public:
+
+#if (!defined(__CUDACC__) && !defined (__CUDA_ARCH__))
+  template <class CharT, class Traits>
+  friend std::basic_ostream<CharT, Traits> &
+  operator<<(std::basic_ostream<CharT, Traits> &ost, const mixmax_engine &me) {
+    // save the state of RNG to stream
+    int j;
+    ost << "mixmax state, file version 1.0\n";
+    ost << "N=" << N << "; V[N]={";
+    for (j = 0; (j < (N - 1)); j++) {
+      ost << (std::uint64_t)me.V[j] << ", ";
+    }
+    ost << me.V[N - 1];
+    ost << "}; ";
+    ost << "counter=" << (std::uint64_t)me.counter << "; ";
+    ost << "sumtot=" << (std::uint64_t)me.sumtot << ";\n";
+    ost.flush();
+    return ost;
+  }
+
+  template <class CharT, class Traits>
+  friend std::basic_istream<CharT, Traits> &
+  operator>>(std::basic_istream<CharT, Traits> &in, mixmax_engine &me) {
+    // will set std::ios::failbit and throw an exception if format is not right
+    std::array<std::uint64_t, N> vec;
+    std::uint64_t sum = 0, sumtmp = 0, counter = 0;
+    in.ignore(150, '='); // eat chars up to N=
+    CharT xxxchar;
+    try {
+      std::basic_string<CharT> line;
+      if (std::getline(in, line)) {
+        std::basic_istringstream<CharT> iss(line);
+        std::basic_string<CharT> token;
+        getline(iss, token, xxxchar = (';'));
+        int i = std::stoi(token);
+        if (i != N) {
+          std::cerr
+              << "ERROR: Wrong dimension of the MIXMAX RNG state on input, "
+              << i << " vs " << N << "\n";
+        }                    // else{std::cerr <<"Dim ok "<< i << "\n";}
+        iss.ignore(15, '{'); // eat chars up to V[N]={
+        for (int j = 0; j < N - 1; j++) {
+          std::getline(iss, token, xxxchar = ',');
+          if (!iss.fail()) {
+            vec[j] = std::stoull(token);
+            sum = me.modadd(sum, vec[j]);
+          } // std::cerr <<vec[j]<<" + ";}
+        }
+        std::getline(iss, token, xxxchar = '}');
+        if (!iss.fail()) {
+          vec[N - 1] = std::stoull(token);
+          sum = me.modadd(sum, vec[N - 1]);
+        } else {
+          std::cerr << "ERROR: last elem empty\n";
+        }
+        iss.ignore(10, '=');
+        std::getline(iss, token, xxxchar = ';');
+        counter = std::stoi(token);
+        iss.ignore(10, '=');
+        std::getline(iss, token, xxxchar = ';');
+        sumtmp = std::stoull(token);
+      } else {
+        std::cerr << "ERROR: nothing to read\n";
+      }
+    } catch (const std::exception &e) {
+      std::cerr << "ERROR: Exception caught in MIXMAX RNG while reading state: "
+                << e.what() << "'\n"; // "\nfrom string='" << token <<
+      in.setstate(std::ios::failbit);
+      throw;
+    }
+    if (sum == sumtmp && counter > 0 && counter < N) {
+      for (int i=0; i < N; i++){ me.V[i] = vec[i];};
+      me.counter = counter;
+      me.sumtot = sumtmp;
+    } else {
+      std::cerr << "ERROR: incorrect checksum or out of range counter while "
+                   "reading  MIXMAX RNG state.\n"
+                << counter << " " << sumtmp << " vs " << sum << "\n";
+      in.setstate(std::ios::failbit);
+    }
+    return in;
+  }
+#endif // #indef __CUDA_ARCH__
+
+};
+
+
+_dev mixmax_engine  ::mixmax_engine()
+// constructor, with no params, fast and seeds with a unit vector
+{
+    seed_vielbein();
+}
+
+_dev mixmax_engine  ::mixmax_engine(uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID)
+// constructor, no need to allocate, just seed
+{
+    seed_uniquestream( clusterID,  machineID,  runID,  streamID );
+}
+
+
+_dev uint64_t mixmax_engine  ::iterate_raw_vec(uint64_t* Y, uint64_t sumtotOld)
+{
+    // operates with a raw vector, uses known sum of elements of Y
+    uint64_t  tempP, tempV;
+    Y[0] = ( tempV = sumtotOld);
+    uint64_t sumtot = Y[0], ovflow = 0; // will keep a running sum of all new elements
+    tempP = 0;              // will keep a partial sum of all old elements
+    for (int i=1; i<N; ++i){
+            uint64_t tempPO = MULWU(tempP);
+            tempP = modadd(tempP, Y[i]);
+            tempV = MOD_MERSENNE(tempV+tempP+tempPO); // new Y[i] = old Y[i] + old partial * m
+    Y[i] = tempV;
+    sumtot += tempV; if (sumtot < tempV) {++ovflow;}
+    }
+    return MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
+}
+
+_dev uint64_t mixmax_engine  ::get_next()
+{
+   int i = counter;
+
+   if ((i<=(N-1)) )
+   {
+     ++counter;
+     return V[i];
+   }
+   else
+   {
+     sumtot = iterate_raw_vec(V, sumtot);
+     counter=2;
+     return V[1];
+   }
+}
+
+_dev double mixmax_engine  ::flat()             // Returns a random double with all 53 bits random, in the range (0,1]
+{    /* cast to signed int trick suggested by Andrzej Gorlich     */
+    int64_t Z=(int64_t)get_next();
+    return INV_MERSBASE*static_cast<double>(Z);
+}
+
+_dev void mixmax_engine  ::seed_vielbein()
+{
+    for (int i=0; i < N; i++){
+        V[i] = 0;
+    }
+    V[0] = 1;
+    counter = N;  // set the counter to N if iteration should happen right away
+    sumtot = 1;
+}
+
+
+_dev void mixmax_engine  ::seed_uniquestream( uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID ){
+    seed_vielbein();
+    //print_state();
+    sumtot = apply_bigskip(clusterID,  machineID,  runID,   streamID );
+    //std::cerr << "seeding mixmax with: " << clusterID << ", " << machineID << ", " <<   runID <<  ", " <<  streamID << "\n";
+    counter = 1;
+}
+
+
+_dev uint64_t mixmax_engine  ::apply_bigskip( uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID ){
+    /*
+     makes a derived state vector, Vout, from the mother state vector Vin
+     by skipping a large number of steps, determined by the given seeding ID's
+
+     it is mathematically guaranteed that the substreams derived in this way from the SAME (!!!) Vin will not collide provided
+     1) at least one bit of ID is different
+     2) less than 10^100 numbers are drawn from the stream
+     (this is good enough : a single CPU will not exceed this in the lifetime of the universe, 10^19 sec,
+     even if it had a clock cycle of Planch time, 10^44 Hz )
+
+     Caution: never apply this to a derived vector, just choose some mother vector Vin, for example the unit vector by seed_vielbein(X,0),
+     and use it in all your runs, just change runID to get completely nonoverlapping streams of random numbers on a different day.
+
+     clusterID and machineID are provided for the benefit of large organizations who wish to ensure that a simulation
+     which is running in parallel on a large number of  clusters and machines will have non-colliding source of random number
+
+     did i repeat it enough times? the non-collision guarantee is absolute, not probabilistic
+
+     */
+
+
+    const   uint64_t skipMat17[128][17] =
+#include "mixmax_skip_N17.c"
+    ;
+
+    const uint64_t* skipMat[128];
+    for (int i=0; i<128; i++) { skipMat[i] = skipMat17[i];}
+
+    uint32_t IDvec[4] = {streamID, runID, machineID, clusterID};
+    uint64_t Y[N], cum[N];
+    uint64_t insumtot=0;
+
+    for (int i=0; i<N; i++) { Y[i] = V[i]; insumtot = modadd( insumtot, V[i]); } ;
+    for (int IDindex=0; IDindex<4; IDindex++) { // go from lower order to higher order ID
+        uint32_t id=IDvec[IDindex];
+        uint32_t r = IDindex*8*sizeof(uint32_t);
+        while (id){
+            if (id & 1) {
+                for (int i=0; i<N; i++){ cum[i] = 0; }
+                for (int j=0; j<N; j++){              // j is lag, enumerates terms of the poly
+                    // for zero lag Y is already given
+                    for (int i =0; i<N; i++){
+                        cum[i] =  fmodmulM61( cum[i], skipMat[r][j] ,  Y[i] ) ;
+                    }
+                    insumtot = iterate_raw_vec(Y, insumtot);
+                }
+                insumtot=0;
+                for (int i=0; i<N; i++){ Y[i] = cum[i]; insumtot = modadd( insumtot, cum[i]); } ;
+            }
+            id = (id >> 1); r++; // bring up the r-th bit in the ID
+        }
+    }
+    insumtot=0;
+    for (int i=0; i<N; i++){ V[i] = Y[i]; insumtot = modadd( insumtot, Y[i]); } ;  // returns sumtot, and copy the vector over to Vout
+    return (insumtot) ;
+}
+
+
+
+_dev inline uint64_t mixmax_engine  ::fmodmulM61(uint64_t cum, uint64_t s, uint64_t a)
+{
+    uint64_t o,ph,pl,ah,al;
+    uint64_t MASK32=0xFFFFFFFFULL;
+    o=(s)*a;
+    ph = ((s)>>32);
+    pl = (s) & MASK32;
+    ah = a>>32;
+    al = a & MASK32;
+    o = (o & M61) + ((ph*ah)<<3) + ((ah*pl+al*ph + ((al*pl)>>32))>>29) ;
+    o += cum;
+    o = (o & M61) + ((o>>61));
+    return o;
+}
+
+#if  !defined(__CUDA_ARCH__)
+#include <cstdio>
+ void mixmax_engine  ::print_state(){
+    int j;
+    fprintf(stdout, "mixmax state, file version 1.0\n" );
+    fprintf(stdout, "N=%u; V[N]={", N );
+    for (j=0; (j< (N-1) ); j++) {
+        fprintf(stdout, "%llu, ", (unsigned long long)V[j] );
+    }
+    fprintf(stdout, "%llu", (unsigned long long)V[N-1] );
+    fprintf(stdout, "}; " );
+    fprintf(stdout, "counter=%u; ", counter );
+    fprintf(stdout, "sumtot=%llu;\n", (unsigned long long)sumtot );
+}
+#endif
+
+_dev mixmax_engine  mixmax_engine  ::Branch(){
+    sumtot = iterate_raw_vec(V, sumtot); counter = N-1;
+    mixmax_engine tmp=*this;
+    tmp.BranchInplace();
+    return tmp;
+}
+
+_dev mixmax_engine  & mixmax_engine  ::operator=(const mixmax_engine& other ){
+    for (int i=0; i<N; i++) { V[i] = other.V[i];}
+    sumtot = other.sumtot;
+    counter = other.counter;
+    return *this;
+}
+
+_dev void mixmax_engine  ::BranchInplace(){
+    // Dont forget to iterate the mother, when branching the daughter, or else will have collisions!
+    // a 64-bit LCG from Knuth line 26, is used to mangle a vector component
+    constexpr uint64_t MULT64=6364136223846793005ULL;
+    uint64_t tmp=V[1];
+    V[1] *= MULT64; V[1] &= MERSBASE;
+    sumtot = modadd( sumtot , V[1] - tmp + MERSBASE);
+    sumtot = iterate_raw_vec(V, sumtot);
+    counter = 1;
+}
+
+#endif      // __MIXMAX_H
diff --git a/madspace/include/madspace/mixmax/mixmax_skip_N17.c b/madspace/include/madspace/mixmax/mixmax_skip_N17.c
new file mode 100644
index 000000000..e1cdbd18a
--- /dev/null
+++ b/madspace/include/madspace/mixmax/mixmax_skip_N17.c
@@ -0,0 +1,285 @@
+/*
+ skipping coefficients for N=17
+ three-parameter generator
+
+ *  MIXMAX
+ *  A Pseudo-Random Number Generator
+ *
+ *
+ *	G.K.Savvidy and N.G.Ter-Arutyunian,
+ *  On the Monte Carlo simulation of physical systems,
+ *	J.Comput.Phys. 97, 566 (1991);
+ *  Preprint EPI-865-16-86, Yerevan, Jan. 1986
+ *
+ *  K.Savvidy
+ *  The MIXMAX random number generator
+ *  Comp. Phys. Commun. 196 (2015), pp 161–165
+ *  http://dx.doi.org/10.1016/j.cpc.2015.06.003
+ *
+ */
+
+
+/*
+ skipping coefficients for N=17, m=2^36+1 and s=0
+ */
+
+
+{
+
+{ 159701761865468764 ,  126615459087983936 ,  971935980882445987 ,  1266425933430727482 ,  330954324138038300 ,  1494390944182608462 ,  985567525599556189 ,  1444303754750856766 ,  1592836147366511740 ,  944566361294427577 ,  6336044263941214 ,  686240801639476028 ,  1446998785307329212 ,  2077656654628400897 ,  458344447554229916 ,  1731527985347537415 ,  1856587133043405783 }, // x^2^512
+
+{ 582338566698838802 ,  782718157798125637 ,  684912925900277712 ,  1271246213246651012 ,  564133275057152122 ,  1254731637174240011 ,  899802366531028440 ,  203611955528219606 ,  2146610763813957982 ,  115972700754870533 ,  2092933746397577991 ,  1136406621159800137 ,  1067897330754616454 ,  1722611875894083236 ,  1630948851615767512 ,  989604548662398765 ,  421680943972696470 }, // x^2^513
+
+{ 617097029144871694 ,  2275108162235384877 ,  1363862259022425560 ,  153831723325577176 ,  1962207167612439146 ,  2011895465932033980 ,  1007329645826843776 ,  1300917526580470528 ,  1171539246226411634 ,  836665263844149790 ,  1026633245683781179 ,  32389807716467901 ,  2296276234884428109 ,  409297319724863974 ,  1911430032287806834 ,  916469329336553330 ,  65677079513594107 }, // x^2^514
+
+{ 1656997175938683235 ,  334384582636879447 ,  2289018900497662792 ,  1053313217223819648 ,  2529561867064419 ,  2085279416854675588 ,  2246944213774249775 ,  917902382070457615 ,  505594560851411438 ,  1748690703652028000 ,  1078212339021286156 ,  170401217082007747 ,  1552855376842911219 ,  2209415219208922654 ,  1174711397611724375 ,  1342069058820567972 ,  577671676013094771 }, // x^2^515
+
+{ 200864199029196267 ,  647584785983829858 ,  449579673559383642 ,  2210795481121719067 ,  1185405708665067088 ,  310118907549337558 ,  1112237647478646030 ,  907316892605650051 ,  1241087052423417724 ,  1441515240516850614 ,  13684659557896782 ,  1932376613205111310 ,  60357690290969298 ,  403792989584565525 ,  102986206885970367 ,  1838760170077135151 ,  727651972830521938 }, // x^2^516
+
+{ 2112835707304810889 ,  583170948234442322 ,  504636752552038260 ,  1287047182580144501 ,  2197579855548731362 ,  1121897458052508302 ,  1952770053431596080 ,  772611988382900586 ,  1985521487785539447 ,  1011576935509998521 ,  367571328103323233 ,  1395044686007468872 ,  398356602519830546 ,  1118155658788636138 ,  1242807300672421078 ,  783603574838596371 ,  2129902332243215295 }, // x^2^517
+
+{ 476494511866938851 ,  1225043666051930516 ,  965084907571074277 ,  48495729686479132 ,  668653785847993343 ,  1110923567965684699 ,  1480594060141743070 ,  1780050366926533029 ,  1286523699078185217 ,  669624786360998719 ,  1865875391054214711 ,  1887676069452920579 ,  9982283268123603 ,  1231842656474052119 ,  1913212077924560289 ,  1062511441187205188 ,  698250642012057138 }, // x^2^518
+
+{ 1052793057808933925 ,  1488114650898155736 ,  1809655216561165737 ,  1766730333219340872 ,  1142295800552777206 ,  1197008472054994109 ,  288311604729845194 ,  1252198368803793342 ,  1700393396484664348 ,  1561218011609178732 ,  1039496900875671159 ,  1199061431134041566 ,  217182459685479427 ,  1965769716664277780 ,  883713133563637492 ,  855967538957706448 ,  182946949964136989 }, // x^2^519
+
+{ 521221301273005220 ,  732177595172654650 ,  1177004708676171981 ,  491026482834387582 ,  1475491185407701571 ,  338551745124704898 ,  322679017279562868 ,  1102042763149957642 ,  1781637400817890058 ,  2230976132808504329 ,  1647749103867940848 ,  363540249136304168 ,  1060964616691451642 ,  856789607085849123 ,  470641990892341750 ,  1281772529585552540 ,  603420272349730830 }, // x^2^520
+
+{ 1803571173151157726 ,  1944596483605820771 ,  213094595441099858 ,  901693591156501692 ,  861915355835676099 ,  635260099576893571 ,  1982413183327998212 ,  1120049474254535333 ,  1505283430373136398 ,  284189413831487026 ,  1271922039532043918 ,  154314556649628344 ,  52549783574158204 ,  584612521732421173 ,  1901458131142776860 ,  1937288889567974526 ,  1212147240165209394 }, // x^2^521
+
+{ 2034564475748864204 ,  2101470381626210311 ,  2288088283864667246 ,  1914963188887983022 ,  1095385391183929148 ,  945106352371129676 ,  1400087958043804470 ,  1843694373265941018 ,  389168898546209530 ,  624173071047154320 ,  1627380173012800657 ,  1094171310315703468 ,  1857806715885218084 ,  1326770924244540642 ,  2245083043975660260 ,  807823933843072560 ,  2287413039178247285 }, // x^2^522
+
+{ 1438462501157028417 ,  1178282907367658045 ,  2151871713069448741 ,  2216564753928829884 ,  1135646926807917971 ,  873250056357106251 ,  1665621526869344532 ,  1732452289298951455 ,  1127204916768538887 ,  1653631650320969340 ,  1992482205762138263 ,  1836865404020205588 ,  1595463001888299444 ,  958431756408955134 ,  2110530659740056614 ,  215688460376624015 ,  1158251593711172802 }, // x^2^523
+
+{ 1733925877085130548 ,  2093555366416251473 ,  2234101572813223254 ,  27150917859587271 ,  2120208241420965005 ,  1740681603888968422 ,  344809595448125388 ,  2297225515944255048 ,  2259042601422705288 ,  872774782224758302 ,  1492625124653803821 ,  1721307140915414448 ,  331124824482108658 ,  120601336292842643 ,  1057772943542838223 ,  904281382213867256 ,  2268149550994572836 }, // x^2^524
+
+{ 1308619451049024014 ,  362949499662291350 ,  507935491259228365 ,  1578906684563693579 ,  1147252812024062738 ,  170426700812984671 ,  2292369007840245901 ,  2278676904939144735 ,  1661943824935578983 ,  1151791566444763239 ,  244664648985115502 ,  788605939290857339 ,  185329244981575100 ,  2234030485122694225 ,  1354886281367682733 ,  1172982108216683874 ,  1673999648423090090 }, // x^2^525
+
+{ 12148567885654147 ,  1372231622618064454 ,  1512411127803286398 ,  91088917975722538 ,  1479647890243287328 ,  720784302152223867 ,  796410783806366571 ,  581618991323966327 ,  850700716339425563 ,  930903615199005243 ,  1411478967795549943 ,  1395258853733301256 ,  1287449265209562277 ,  2138600707136677508 ,  1952689678998032749 ,  1536157936222739421 ,  1365983300736022067 }, // x^2^526
+
+{ 1614268848379872789 ,  1977966741412095904 ,  1868426217239032474 ,  1603841034660414664 ,  1538146046227581261 ,  686243317632590721 ,  1825936888550743483 ,  1880896341340126997 ,  368542569942961260 ,  127908944176116453 ,  2240186396348061575 ,  1299238346503882153 ,  2102660872371583004 ,  1938050823525276998 ,  1611779572790548487 ,  923339421522954197 ,  1640317376822349283 }, // x^2^527
+
+{ 1068274751244808317 ,  478761882933885805 ,  146343578721548466 ,  2122787414678627286 ,  478616655821124071 ,  1223652533398739090 ,  1178267270852093552 ,  1585682559495334861 ,  448375108812347778 ,  28390211070681947 ,  1690010627725723003 ,  1631513459519602506 ,  344109887337730484 ,  1228367175232395654 ,  1245714930611054419 ,  2076109931139359104 ,  1028725609142702754 }, // x^2^528
+
+{ 922695624671435383 ,  1212395461257377598 ,  2004462039540699849 ,  1640073156934071824 ,  2048256174899844882 ,  1784439921745628237 ,  1181704167380381729 ,  1002268784211250265 ,  1203700422140227637 ,  414802973276697203 ,  1506620899880171251 ,  1653959500223714911 ,  761577027275676953 ,  1394862681050805826 ,  283665158959068873 ,  1674581989169655504 ,  2170268137063855658 }, // x^2^529
+
+{ 1531823228490667160 ,  655178204437606577 ,  1143440975987174633 ,  1469225472530244216 ,  1459815855107895954 ,  169158503015419167 ,  896992338984130145 ,  2224392167406479582 ,  334969241330017474 ,  1210934779672254914 ,  958852377915238487 ,  916143882424925025 ,  1228314367757567546 ,  576985231081272085 ,  69708611488908557 ,  419353864719896559 ,  2063012113607020176 }, // x^2^530
+
+{ 1561213794756615274 ,  1035561623111560834 ,  236042828876112400 ,  2227997009011010087 ,  69191692217035803 ,  680016809348368644 ,  304720676536884446 ,  402216461956911069 ,  579355410310536007 ,  2067077610395244217 ,  850246508783344249 ,  1267387033405622409 ,  381166642361949681 ,  745560785470896920 ,  475789434723626541 ,  624301081669895306 ,  542440096530351799 }, // x^2^531
+
+{ 292947975657821923 ,  1975118192564031841 ,  1096919547621493467 ,  1963319303199452777 ,  1497672437671054944 ,  587569660431161322 ,  1640139998254479090 ,  2075259177632397124 ,  1660504006066529017 ,  2098448192586763026 ,  2110557029504033173 ,  1649560687483066215 ,  1687307417474238546 ,  1953341749861443552 ,  1238894388766834138 ,  1741784409005434700 ,  961090363500770371 }, // x^2^532
+
+{ 1383869753498066511 ,  1233178641086644718 ,  983298541932470614 ,  670427917276364832 ,  1594590175512954791 ,  1215937832423358959 ,  1903793454103261000 ,  2172773248760515308 ,  536919802769335489 ,  1192052568730958579 ,  1911240856897478532 ,  474283681972418461 ,  1641036102716815339 ,  1411044907503618945 ,  1001175747401094733 ,  301211985765609746 ,  87149012356549448 }, // x^2^533
+
+{ 261502789562703795 ,  2039292976033461220 ,  1383214847400482168 ,  617091957131872918 ,  11536068257558752 ,  1594678644021354556 ,  75861344789992516 ,  571299992729856381 ,  1893030021771949746 ,  955316100521264931 ,  1984661955054505780 ,  2011703935837291735 ,  1088722461445557653 ,  1179450192264585869 ,  2187858293377977903 ,  925337903456050826 ,  646310246691765740 }, // x^2^534
+
+{ 1067110047296065055 ,  1681438952663914065 ,  1912868211773513420 ,  157328871056550869 ,  1606539524193585050 ,  1826960759022391567 ,  333004198356937619 ,  1654504867491237582 ,  1675004459630487926 ,  71936368060300 ,  1918748055540842111 ,  1831638571780226898 ,  1740919230167947613 ,  1865657023192411640 ,  1938852440014266019 ,  294008211664586131 ,  1233347209473542392 }, // x^2^535
+
+{ 1793671892574625082 ,  1130285389461734712 ,  2256576488166743893 ,  499935214652909272 ,  1154848893191704610 ,  2081455429764661190 ,  132180046418477201 ,  2017092125714317442 ,  1267352468781603191 ,  1320951561119166755 ,  1819380377385314624 ,  1968657260504504792 ,  795175392323167167 ,  138949840583543195 ,  2205578352394323161 ,  2162095450973436967 ,  144991897502595141 }, // x^2^536
+
+{ 683024508151477474 ,  1048650848431764111 ,  222063671742581704 ,  1033879240711051787 ,  2065492192337124448 ,  1040278918101037386 ,  1671617116839141169 ,  1142933505587090141 ,  578155580871550072 ,  1397916094430897857 ,  1428297170374867004 ,  343607952975225879 ,  2212099097846398741 ,  659520961168888597 ,  1099474129656084832 ,  163872983722101162 ,  469023960460927874 }, // x^2^537
+
+{ 1156093888444545360 ,  477156936151507247 ,  340730541722292373 ,  1631554541005175777 ,  186957487417408624 ,  926739184035859811 ,  1357420454563306968 ,  1018117004392525830 ,  363668834887049598 ,  562598157996426449 ,  1477447398167639388 ,  570432714234334068 ,  797165670329651828 ,  927307345530529305 ,  1474634744709399927 ,  2245157369991218355 ,  586948150323806327 }, // x^2^538
+
+{ 267980979185924701 ,  1565212575331069567 ,  2036070690586106209 ,  850842928652802398 ,  1084085447966709199 ,  53832993778516816 ,  2236093215844622791 ,  361324823691972021 ,  673008249565884673 ,  1820092500335360408 ,  2155516028757275328 ,  1127096462606937572 ,  1358543969239663842 ,  2250855440900689232 ,  1750872323705159361 ,  1362418718801185396 ,  766489969361488653 }, // x^2^539
+
+{ 259726699384885201 ,  1899046140654817590 ,  1306363085361051998 ,  2081068431945330598 ,  85108856372668378 ,  1324347710267860858 ,  637751183869623188 ,  2248344050276940170 ,  2150579538601288819 ,  1105694420952032875 ,  952722029405414348 ,  1555403505154017419 ,  2124427035047642831 ,  2043566203745579024 ,  2297745588334532887 ,  383567419144690253 ,  295720540923262008 }, // x^2^540
+
+{ 1861487178452566675 ,  1408373324284705278 ,  1763971337054095719 ,  830694774025045496 ,  1222854782641346372 ,  1367429629279521808 ,  646871018425974561 ,  594901312199197009 ,  1347745885562466060 ,  1872975180162717691 ,  264288217725438485 ,  1587326205070728608 ,  636368619535248259 ,  1345295891447307881 ,  147325724921249282 ,  451151671184877984 ,  369152850758367779 }, // x^2^541
+
+{ 57601454029724418 ,  8965969940832902 ,  1015753881352776595 ,  1272955433813211634 ,  1500099915165437510 ,  1505908015386217991 ,  2200135170747962774 ,  898929787588318880 ,  94174721104731682 ,  2001843482426425606 ,  1296185933469804867 ,  2067771350105959635 ,  244437767883759705 ,  115304829559788030 ,  1570793303180703663 ,  2075829501594571540 ,  1142789736568430590 }, // x^2^542
+
+{ 1666339398200583535 ,  352498459116081565 ,  1296181205638634316 ,  943151467845475937 ,  1576590058032732201 ,  991197898685844101 ,  172929644755996237 ,  348073447626095407 ,  1119639777406826568 ,  569438225411593708 ,  2284730591697077199 ,  62008449878226744 ,  237477582274396799 ,  284434691435252140 ,  1765995421169394732 ,  1594311192946755506 ,  604156017409724002 }, // x^2^543
+
+{ 495281382117458539 ,  205408453347821044 ,  1681558300394231421 ,  768174766423516475 ,  2178619873339022076 ,  765997180179026402 ,  102658969211254460 ,  895017413823772592 ,  415845558640901473 ,  601933504588026198 ,  969782745266307970 ,  855044719798584000 ,  2201199786666131647 ,  1734600819975899129 ,  1006570536931532428 ,  1130105133424164395 ,  1127782802378551477 }, // x^2^544
+
+{ 2279222264094136096 ,  918140621636386578 ,  764051053229027870 ,  284578302438407812 ,  18442247946635388 ,  31558521455257101 ,  2294641003511473164 ,  1698376578572281576 ,  1937426563738038511 ,  809127821329680664 ,  1605067175404678954 ,  189108927488132646 ,  2042036507825496222 ,  311247532621290928 ,  1403615118897576721 ,  1715491077248039654 ,  1885408148480015203 }, // x^2^545
+
+{ 1303703719061282226 ,  2273657250242045723 ,  1357599208936374367 ,  635087130862478000 ,  1079733394689912801 ,  1118552166016645994 ,  606420152917155530 ,  1461390879031614972 ,  104855646611743091 ,  1155586802792637483 ,  1332607042374593381 ,  485954066105103318 ,  755622945227398111 ,  100703703472763330 ,  2104688739715121722 ,  884188331727454949 ,  634077400952306173 }, // x^2^546
+
+{ 951241189979241574 ,  605016052445461627 ,  2017048237067736326 ,  989577001336697190 ,  729360432477738078 ,  1891686993869218661 ,  1178287733434616934 ,  2257768800220793298 ,  386183991328298626 ,  1256576692082166247 ,  1989644154430757909 ,  677803830671194177 ,  233612802004492344 ,  1711459201958972728 ,  1493606889466009749 ,  2046914649687863407 ,  621460642415971297 }, // x^2^547
+
+{ 1751267799283060659 ,  2256655950470869985 ,  984875549537277532 ,  2146590585318069342 ,  1539686266680012185 ,  25924698613741241 ,  633436588936544176 ,  1549785231759087109 ,  934442226797007335 ,  1754392592879566030 ,  1106511892592859350 ,  2133384541112231515 ,  483183511756818991 ,  1106565501063528265 ,  1107187069130962112 ,  2075632961004138802 ,  1789791698038574768 }, // x^2^548
+
+{ 2168119975504999102 ,  12966540192424488 ,  681037686628849373 ,  194933206197471208 ,  1333929840850796338 ,  2213746149905142918 ,  1379421507825630221 ,  672612248824101656 ,  335325159518504203 ,  1891726571269569489 ,  838501400504296092 ,  1301803583046081493 ,  1750320150554937404 ,  67761357879073926 ,  1985524641824792115 ,  1623610528663151019 ,  1386176498075642165 }, // x^2^549
+
+{ 1725373491573591435 ,  1047415675966614123 ,  1624486786194930396 ,  1270371958501591477 ,  915493804007886656 ,  2289688909829684654 ,  2049603525909314516 ,  148613066066958793 ,  726433066490138797 ,  1727275655898645687 ,  1975538226366019203 ,  394922847880149584 ,  462814805660060890 ,  1087564105835469854 ,  293119996590319709 ,  561334965874644379 ,  1159197305565422760 }, // x^2^550
+
+{ 1167941940931666302 ,  186514058829683028 ,  284351732075544532 ,  1512524808600673716 ,  43881750227659649 ,  1609741131150076188 ,  1125658924019416678 ,  1011534154384997338 ,  732738003772919984 ,  277064895441615462 ,  1622677204527889015 ,  1938173807800443399 ,  2063908410477010918 ,  48488552879373011 ,  1181372505578038205 ,  135655409305671122 ,  1938059649851535706 }, // x^2^551
+
+{ 337206759525681069 ,  1755948102219815566 ,  678765415558248151 ,  1298324546232590789 ,  484850096725079792 ,  1257105102233963714 ,  1905953290815779223 ,  1631095014301272964 ,  27718500263208979 ,  850594049252190048 ,  583977122492931140 ,  1805594056761021711 ,  1719584503850087741 ,  86179004867644487 ,  1270985715490800843 ,  644779453772114088 ,  1760476987252969926 }, // x^2^552
+
+{ 853995884080236938 ,  369111089848434712 ,  559743835711047174 ,  1594116816433474826 ,  35842475217144553 ,  2171235383636811428 ,  765712905092749171 ,  941173958645971457 ,  608455896272321575 ,  560604720551238557 ,  1845649540169134916 ,  1718920359984151079 ,  1401755957661962268 ,  1423907959444825866 ,  1098311071256242778 ,  1528916661846841494 ,  1508764511499232241 }, // x^2^553
+
+{ 1508791960279504205 ,  258737323991262817 ,  2078642398527055752 ,  2047350895402792396 ,  1386159783711778842 ,  2027047763369412816 ,  58440814829283402 ,  771612188948247815 ,  1023214498156083887 ,  844220082620919469 ,  868367506978956802 ,  460654633908296364 ,  2166638508094415965 ,  2152060517892743275 ,  2195460306854592144 ,  541094667496643303 ,  2293587356561486979 }, // x^2^554
+
+{ 1214656342962583818 ,  757702747146927395 ,  2189480375257771121 ,  1816527465277679192 ,  138815664117127512 ,  1431350214961325326 ,  1187978151676236367 ,  1241806570983467832 ,  1284450519130025998 ,  1613889771871333389 ,  1061608822774622375 ,  381417471969594147 ,  2180713817066758069 ,  1573268990076834701 ,  905173682016321866 ,  1950470428221045886 ,  295485711949393017 }, // x^2^555
+
+{ 819753035519543635 ,  538715919221404699 ,  2219429537481971437 ,  752393850271742164 ,  77092035906327214 ,  1114501521607776562 ,  1781083659997559466 ,  2154696941256077128 ,  1218985834766672241 ,  211616781192221457 ,  216411965392177616 ,  1670594910477530341 ,  714505138692032780 ,  2103941244229463221 ,  2088317230298558512 ,  1289168057502589286 ,  1125637208200814412 }, // x^2^556
+
+{ 1962227768978994180 ,  1690314396114180321 ,  375465114688955302 ,  1798066679445937030 ,  462818969271150710 ,  1000261197673181741 ,  449380920277416447 ,  2291003141210787002 ,  1295168501047080160 ,  96762399469569918 ,  673992173479089483 ,  22635411461259319 ,  690222719954237584 ,  336145523606337534 ,  934574030654348792 ,  1843383288751365671 ,  1895192556383603959 }, // x^2^557
+
+{ 22617970229773822 ,  1482887891197581503 ,  190614907712811083 ,  693572772014158846 ,  2013130701640461017 ,  2080394974487399438 ,  1946114308346089259 ,  990495636788197392 ,  1782544219386492119 ,  661191434765959775 ,  2044311906562246966 ,  1069940178504200830 ,  1955224049483069684 ,  2032802635089013840 ,  1312150198687724551 ,  1790320974584854699 ,  1722451636719440405 }, // x^2^558
+
+{ 844270734204704640 ,  144569156478286025 ,  458434355726727522 ,  54107521005723860 ,  1370238437704355388 ,  635959790870056500 ,  580190401330281055 ,  1190738125435083361 ,  1508034736034112394 ,  384356207961213908 ,  2225422022608450401 ,  1307978291279169425 ,  1836303719023016910 ,  1688394455914172152 ,  1145318805385872977 ,  28214091532568837 ,  2271521727551942534 }, // x^2^559
+
+{ 707823463534947110 ,  1322853985703915624 ,  1283965375707685336 ,  1235846354925013921 ,  2241994185339373382 ,  1492018342988994366 ,  354756557339909458 ,  2038382297798252220 ,  1998256038116472547 ,  1191608355457210551 ,  431624398742182860 ,  979989445521359703 ,  1874604718333228361 ,  863699730728041575 ,  303002850561444625 ,  351750811939826058 ,  2086150517074464097 }, // x^2^560
+
+{ 938650030581340927 ,  1484495848437960060 ,  906987191271901948 ,  602861889366282296 ,  829217328598139608 ,  1170659179365524445 ,  136537691170306827 ,  2263710078069587303 ,  1916941180571676310 ,  601979432897153599 ,  2075793575506860756 ,  323119178978530887 ,  7962577560787577 ,  959899426480760280 ,  465731965934497666 ,  485390457639864074 ,  182578072602393068 }, // x^2^561
+
+{ 845056639071703739 ,  2022501832369820141 ,  1551830514420859659 ,  712651911384419505 ,  1901754482233053508 ,  771255057924942569 ,  1624974410054352307 ,  1095555020243807180 ,  48292506949982043 ,  1924628272072369018 ,  354506579070195883 ,  2074354629523999876 ,  481327135146565647 ,  1980060303227111264 ,  374025770449242754 ,  2197437559867642882 ,  1347989415958383967 }, // x^2^562
+
+{ 2284531938332082408 ,  2081446697857700129 ,  407946796335100046 ,  2284742117207953792 ,  310471079731850760 ,  1290351582175218317 ,  356273760688651989 ,  1910139122966702303 ,  380631662657835885 ,  2287374039009394440 ,  485797501634087358 ,  900166727954042323 ,  1169554169155638307 ,  1802101457472636718 ,  474217730229244850 ,  1004410789585399243 ,  588639270649921383 }, // x^2^563
+
+{ 1502830393741044692 ,  502885788914291282 ,  1000641636223469702 ,  599101320428395707 ,  2179406366865940840 ,  358024368445093522 ,  1564063948136384394 ,  2159233836362830609 ,  1080529492271670836 ,  737489746809287059 ,  779539060254650511 ,  567680696243734902 ,  803303944361973640 ,  1778879099887965284 ,  2245121666542231208 ,  280989822576167540 ,  2256030916132704537 }, // x^2^564
+
+{ 1678854909602034917 ,  1362544615060522442 ,  1338371892984320209 ,  1236016995662152262 ,  2266586199968442015 ,  1429545066557982811 ,  1057272578263945968 ,  1902362138019567296 ,  2027249515853796780 ,  1379458255463795129 ,  715456741316196206 ,  483096191538955481 ,  562452275970679447 ,  1899368762554121096 ,  1095218941132187764 ,  1182616069429316286 ,  2041204243971349283 }, // x^2^565
+
+{ 1470384206682188393 ,  993056282930013636 ,  2186884113482342484 ,  1568683072711048719 ,  2141576651319216540 ,  1545634359639056795 ,  1451267712100583056 ,  583449252763690840 ,  743076343569515666 ,  880735478838481309 ,  1139615309806234587 ,  1405631565601779814 ,  1162183658026694832 ,  758880015783224342 ,  186132385266916884 ,  299496669995968613 ,  2115484609912800200 }, // x^2^566
+
+{ 2269185097674681740 ,  368408303414772947 ,  2272053720551708375 ,  1632770988466577149 ,  2005375285078385494 ,  1280199055327212680 ,  1821048605091132853 ,  1567765420710739271 ,  1887232479893207919 ,  2230675230841783167 ,  1361156385383252011 ,  87414071150394174 ,  611216773833872434 ,  1860713875917257321 ,  1409775498386462397 ,  1690326679721511537 ,  83275546754626406 }, // x^2^567
+
+{ 790112276274407445 ,  769228027816192054 ,  2205289911751012082 ,  696806322595942722 ,  1027622205316503460 ,  1843245508230762915 ,  181348331184190272 ,  1685443775553823790 ,  765291022909566823 ,  1159224061954534775 ,  1753833798485424233 ,  1131135065567349204 ,  1004980355827594316 ,  1518506610289278726 ,  1077654190664918382 ,  2049836184231549381 ,  186534172155365157 }, // x^2^568
+
+{ 2214978522047075314 ,  1523958231722460737 ,  940540633818853829 ,  1502207690061290977 ,  1973166365192858632 ,  834167590580214713 ,  193198187070116104 ,  1200481096199012868 ,  1045382432042256669 ,  1043046955924667651 ,  1764037894194441387 ,  326214257353715662 ,  847610851336013740 ,  821380410810508534 ,  118611704753506829 ,  212870724809513245 ,  1049669518188822039 }, // x^2^569
+
+{ 1186686984731968114 ,  492485282093702769 ,  2084974668148141701 ,  1954670681979820725 ,  158025799762360442 ,  1546024887967767275 ,  321250855908273132 ,  1839134851204623669 ,  2276195450093740401 ,  216550545714992720 ,  1147405869335269365 ,  1691809800568592195 ,  92303605959988217 ,  1444340652242832966 ,  1276065726666195728 ,  324774933653588943 ,  1312430778251072843 }, // x^2^570
+
+{ 1914401401876468669 ,  1648868062260221854 ,  1464059991570664252 ,  875976384406088755 ,  859422856357056602 ,  1864309596721585546 ,  516857141291982496 ,  464095730580902783 ,  934796258134400377 ,  1817127289097008540 ,  655971001935077689 ,  680288930067826216 ,  229080197392071732 ,  921859076031383827 ,  2242366335784005876 ,  2231689867120356867 ,  102494160827714957 }, // x^2^571
+
+{ 1047877418576672693 ,  850258032575291687 ,  545055359259387307 ,  939787324567586678 ,  588678011951547088 ,  397057856974487431 ,  226066134643513040 ,  2082177109290521202 ,  866921181848612589 ,  1377167370849823007 ,  923546628885112186 ,  1369468709683575581 ,  163347741879505736 ,  63898233712570047 ,  1837149445977275412 ,  81917523037053887 ,  1724871657333737359 }, // x^2^572
+
+{ 2302595391708576272 ,  160985864775068781 ,  764451814383377581 ,  277491655221007957 ,  2065357052112797123 ,  2022995433719782905 ,  1992649094929659709 ,  534121146273082282 ,  2156280075941140984 ,  1374921692284132634 ,  7490180772772763 ,  724537764516397075 ,  1458937057231627989 ,  2128708058250265224 ,  1343524011236415024 ,  1632380201084339158 ,  218723970886479505 }, // x^2^573
+
+{ 2117760628529948776 ,  2132935432038525560 ,  768781595656590901 ,  1400974401324013175 ,  1245919916200029631 ,  1432559837530198229 ,  1163713351897028899 ,  1923146399637847839 ,  1213782762495170404 ,  541538612171639467 ,  1998922401826244596 ,  1231856100509046580 ,  579861832524317272 ,  342657654818204433 ,  1849901794740458260 ,  2279916971309950571 ,  760499761682595396 }, // x^2^574
+
+{ 2133951096918945939 ,  245427172150052114 ,  981217690063201632 ,  2240543031382098229 ,  72038350069308602 ,  597184866817343336 ,  328644676877024543 ,  183469578860761428 ,  1484808128953021465 ,  1223798115302938004 ,  622162521478669994 ,  1190522948279756621 ,  1770221081272204495 ,  300642277489993017 ,  323299607262829861 ,  1121130797903609296 ,  1341151672362966624 }, // x^2^575
+
+{ 978460885885818525 ,  1015163128148908540 ,  1977562332100546610 ,  2245843532918136978 ,  52755975506265961 ,  105385627927894157 ,  1072178687697178180 ,  2189119972792037580 ,  1846365892857759522 ,  509344149142820345 ,  1764914861664878049 ,  1812095995141658018 ,  972953132412161570 ,  1195708733082367698 ,  49856206632166450 ,  1968416938611340445 ,  2224505860268693924 }, // x^2^576
+
+{ 327862488024780357 ,  1457219179227842324 ,  2230322439686492524 ,  1014759634574897984 ,  268418936182947756 ,  1502956125465763704 ,  61942509184664792 ,  1760829907047336797 ,  1772759420678440112 ,  213631264477227981 ,  1260050930099550098 ,  1356809535018841668 ,  2084415297676571737 ,  392165946698433004 ,  823317575518425201 ,  1619251312183556482 ,  1475671506190286260 }, // x^2^577
+
+{ 1056862459352327285 ,  650761739286013201 ,  711866770850654745 ,  995772669095949341 ,  638628037366816307 ,  319488467144926010 ,  146863433379310830 ,  1101916435400125848 ,  1403815239053246876 ,  1038453114803570080 ,  467259301797897322 ,  1758832454145478316 ,  359553557196445664 ,  1487375835680896666 ,  76291938912745016 ,  844090192831690228 ,  793442078170731733 }, // x^2^578
+
+{ 1952837235717320563 ,  445430339419636400 ,  122525557445898274 ,  735595189125576324 ,  1842210194389261500 ,  1795223307042139024 ,  751864575701713084 ,  1119184159833325999 ,  700766534388037595 ,  393549636253210822 ,  1662606888175396198 ,  777088695141990512 ,  1949046976963730503 ,  659714914447882016 ,  1701619511100158747 ,  2106648518620441962 ,  2230177471130012298 }, // x^2^579
+
+{ 2242830595691081287 ,  2153702821027923023 ,  2056432480338738486 ,  290410908425161907 ,  1981631243730724167 ,  1568909835655310667 ,  1868498486537085636 ,  443377197549410236 ,  1989171328893245017 ,  832914678880127633 ,  122754022462592655 ,  2074425172477633432 ,  450810432342594695 ,  1304742577811847983 ,  1915105538864887632 ,  398240897735172221 ,  854257484180168406 }, // x^2^580
+
+{ 2200340348805749851 ,  459904691507133993 ,  1398407253135138579 ,  2223256024161740360 ,  605045994980733172 ,  1376873737656823533 ,  2154922394517147235 ,  1733951063198734927 ,  972687456101926553 ,  901838082782677970 ,  786315965699770073 ,  220890454000213916 ,  1082354731585990169 ,  1575625362271244736 ,  1934345059883646897 ,  610599809596399831 ,  342994941095873735 }, // x^2^581
+
+{ 135197059374104208 ,  763606053289129943 ,  1762314800446038038 ,  1958652087897421237 ,  1759073139864458242 ,  777882068437388702 ,  1929274536720142390 ,  911476487406644691 ,  2214565556314441473 ,  52312220897893321 ,  1846728531947717379 ,  2200466005983704615 ,  262302790665678959 ,  2226903910074888712 ,  2217869614731681025 ,  56279709804267946 ,  2290324007417928442 }, // x^2^582
+
+{ 2033888752861355042 ,  1032905065687727870 ,  1228320602981295572 ,  507759480115032778 ,  2286473996360684281 ,  2214758482416116553 ,  1809366838165194032 ,  1058372692666396235 ,  1752374720633908247 ,  783225856075389703 ,  1873729037970411918 ,  515708137390539649 ,  1758253278936485521 ,  2232321656689971091 ,  54123739319322403 ,  1645021409194515326 ,  128999663975261693 }, // x^2^583
+
+{ 1548696340275626821 ,  1058994430734061174 ,  975903143090188182 ,  1753221425369032354 ,  642770173718849415 ,  72871801728857389 ,  1532366481913853262 ,  1685179220297983661 ,  221740024932865141 ,  572304397314497403 ,  1337436598601163606 ,  615711199196774965 ,  1856428633974911262 ,  2091160220746353073 ,  1268041182491440797 ,  1823672822413606168 ,  85882960870940069 }, // x^2^584
+
+{ 1596386067990045239 ,  1082049082776347514 ,  964167595689777765 ,  2196584899289686937 ,  1775478739112084197 ,  668585706672723348 ,  148931594902619061 ,  1989674837193253753 ,  548537612282770315 ,  1196600058755599499 ,  40994948219245722 ,  585997221548702784 ,  1852642327512527056 ,  1962368180715517868 ,  647388617327737622 ,  1259769006990328711 ,  760195312392461746 }, // x^2^585
+
+{ 718683317139573302 ,  320255468933257119 ,  1602065008785819027 ,  811539625884985377 ,  473380658706433597 ,  1912262500486946501 ,  470045053505719659 ,  2044632105828108925 ,  1926951883673524133 ,  2227261383404374216 ,  486603446225142281 ,  663210447975562680 ,  939660586393225605 ,  285150104962906352 ,  1736759590302217217 ,  1557356149606075095 ,  202888193082920775 }, // x^2^586
+
+{ 1704252131283387210 ,  749862216125654467 ,  1151691390735685357 ,  1645735627001577084 ,  1704474438618968650 ,  748171802247682061 ,  2192661098187053709 ,  687285052039352483 ,  854812886608086510 ,  1254425642986060904 ,  1755448133247766409 ,  2011170119775364505 ,  2238140979402216154 ,  430682510031291698 ,  1344094596771993123 ,  1487316607310729340 ,  572027053290063200 }, // x^2^587
+
+{ 490145643447149312 ,  139932608980010944 ,  655464257025259880 ,  394296759589777370 ,  505632773723135101 ,  1661411471151623719 ,  2151487638856430824 ,  1405022191175603197 ,  1752924822934387289 ,  540705285689857192 ,  321557534615604425 ,  357668350120260760 ,  1216461776875955281 ,  227959942379372288 ,  1843268809411056088 ,  1278964120834620296 ,  1475031348508152943 }, // x^2^588
+
+{ 77196441284986546 ,  2116432279354126192 ,  477307664737255679 ,  2054284248076058866 ,  1902595507518183399 ,  892309383355817071 ,  632707825296532637 ,  2068144493364120654 ,  2269464712931165434 ,  1869550090717623995 ,  1230295724219128689 ,  265213902138950946 ,  2035223329405422145 ,  1942656580078323795 ,  1193709919134967859 ,  2129017786957860238 ,  1381902932849007959 }, // x^2^589
+
+{ 1250515786673063664 ,  1170029995563455371 ,  1875333685750981198 ,  202301072852130868 ,  2043478147118479243 ,  1430654884019451218 ,  1932432698651776626 ,  906545901346634684 ,  300667563272986237 ,  1419422183404003464 ,  1991379549654933852 ,  647455449394880739 ,  694191473419776108 ,  276654368069638649 ,  327518818778853777 ,  776612031552014882 ,  754333675021774044 }, // x^2^590
+
+{ 1494337748978015003 ,  2195753570178835409 ,  1756395833112881603 ,  1391513303657053978 ,  1469720991718755668 ,  1385546023774706695 ,  2174495928889154879 ,  1845916193012260560 ,  1508599521784450080 ,  1646177877887412686 ,  1600214236327823295 ,  481492506494972498 ,  607939032270197735 ,  1035787083329757142 ,  1867038566495963823 ,  1327899656420714797 ,  120455193991568778 }, // x^2^591
+
+{ 1272270033586499361 ,  1355938084262029033 ,  1130547642931602057 ,  1571210037203668042 ,  109838764087692662 ,  665089929034619866 ,  1477929841071085716 ,  354404882207113353 ,  485207799804145872 ,  295271887824129191 ,  1434571796278931709 ,  1640884996563085247 ,  349639709072954139 ,  2048198598052923555 ,  524893036895320686 ,  2241615913147991183 ,  1476041495950183974 }, // x^2^592
+
+{ 1191478803597513079 ,  179928617987606724 ,  1893767020776627979 ,  2279764014883753334 ,  2000307147025567411 ,  1317337023680470054 ,  2103432586600253524 ,  72635103884479033 ,  2120393133065940248 ,  630577942926446005 ,  184690001576938472 ,  724734720983686380 ,  2222559033240707306 ,  957823503700662442 ,  1530455664164510916 ,  1897774219578829063 ,  1565275869726572277 }, // x^2^593
+
+{ 802169516563427070 ,  2266388954150660250 ,  1016109430849964502 ,  425155158957108184 ,  1509177586061010059 ,  735042889547641844 ,  435841014763623437 ,  552466077216086185 ,  1380304476103110982 ,  67287297735634201 ,  447647206919187346 ,  2068383112520459072 ,  1700683878491979663 ,  876911219219450197 ,  1724440678224620967 ,  2161646349907542457 ,  1646103618347380320 }, // x^2^594
+
+{ 612487302971671307 ,  1375201488802854426 ,  381927427343438951 ,  522937542366715499 ,  344542446499511062 ,  2207818055791032517 ,  1589757712162343309 ,  800577551503529820 ,  1901158862710281177 ,  282526012842555873 ,  978266672216571827 ,  851384948740770681 ,  1847198461333883874 ,  1701120025203313330 ,  1579726260607682225 ,  425398491832428623 ,  470525020962949424 }, // x^2^595
+
+{ 2088490725567484178 ,  485853289338907617 ,  1031877816384950332 ,  2204970857389602283 ,  1742289513696508421 ,  1887097317304814932 ,  611739697091118719 ,  620524462574164777 ,  866744533296715680 ,  113249287824674575 ,  1094771911813698402 ,  1113977650240159227 ,  597674782116109822 ,  745608676499086192 ,  1626698293359327222 ,  1428826314680818479 ,  2215603110047911506 }, // x^2^596
+
+{ 376585604887633526 ,  1048922386135384477 ,  2142424598885808255 ,  675251024248975850 ,  1545330669815065990 ,  1857893574217764638 ,  2048982207690498307 ,  1533653988510679070 ,  335130892749828677 ,  104131490905494608 ,  2194266840041437356 ,  683616020969961024 ,  2038987316895996656 ,  917297429323765875 ,  821082630191764403 ,  1666057951671669384 ,  1130032032943836952 }, // x^2^597
+
+{ 1041829149523127110 ,  1951902012535484719 ,  2152661804693326326 ,  1269908279296181874 ,  971953632796866049 ,  1266540697681616989 ,  2071865642155459708 ,  845692533418042777 ,  1241457687074422224 ,  161412923983560739 ,  390632925491962906 ,  353948016710720772 ,  2026700469979616031 ,  96483949026756707 ,  1773219759585896547 ,  993715984447404907 ,  1829150519300709735 }, // x^2^598
+
+{ 877798478622581516 ,  942092329022409512 ,  1965158005364305044 ,  1441028099309753513 ,  500071448451573682 ,  1840680567469132321 ,  1468971950339879193 ,  1438547156759802726 ,  549595496973813971 ,  390734351826704403 ,  1073287705054298510 ,  48924576386470253 ,  1989283157579117827 ,  187450184005400512 ,  1949224458829795298 ,  1951496711030456157 ,  2167858042677531808 }, // x^2^599
+
+{ 1996597802411685974 ,  114349701898949545 ,  1084220018985830117 ,  27462757590845683 ,  1081095686730339174 ,  2059987003957425282 ,  1321806358033656174 ,  1828830847664378831 ,  1455949950124722259 ,  477041639297859939 ,  2295506231622178555 ,  1152831764620264515 ,  396805345192155178 ,  2211812614000847756 ,  1962229985869550296 ,  848640427896833485 ,  548986687197667440 }, // x^2^600
+
+{ 491594119422354073 ,  1742453895833966377 ,  1525602929914324406 ,  2076597071019376631 ,  865550929971279566 ,  997883929562138029 ,  610577039659314246 ,  1003718639233780963 ,  1951174507028809908 ,  1185710764340213666 ,  1693015712372426191 ,  317455835845832474 ,  1848795541252464763 ,  401332635199514101 ,  224456822593358236 ,  742409296289707819 ,  2118807906613325507 }, // x^2^601
+
+{ 1293009127273644989 ,  1124274480353253667 ,  1810129434176838306 ,  259729045897854950 ,  1468101939413453436 ,  1111947300756487949 ,  501989537109533070 ,  1388491866061828170 ,  1312270603637319867 ,  1615629996017192524 ,  2134426516814571448 ,  548391958129447766 ,  1852805539182613949 ,  562903002408545356 ,  1582149854427086182 ,  1136405868045684314 ,  679089614531440081 }, // x^2^602
+
+{ 786115001012233217 ,  572934901393961574 ,  906662352794537253 ,  1146052338661137694 ,  1717117579480047687 ,  2238979311119472866 ,  1282072454711511764 ,  913074328950615201 ,  1508062462715721091 ,  585179548016207318 ,  2263240538685112497 ,  212602680696443031 ,  2286609014980268736 ,  2160223566289961256 ,  457719410118216561 ,  892534882541613344 ,  1802928045269113919 }, // x^2^603
+
+{ 412244890005542545 ,  882317329387396531 ,  1194216517008066508 ,  2113679728861695547 ,  754826510247388584 ,  143906396701128982 ,  1624499088946075077 ,  1517734354505524658 ,  1319665375831884375 ,  108068462219720639 ,  2061550944874121143 ,  1387192717449700004 ,  1784206536018821508 ,  1531349906234275665 ,  37851629185782003 ,  1964408593925350957 ,  675543575437141874 }, // x^2^604
+
+{ 2224871078409190905 ,  1109593689354222314 ,  2282026642284422235 ,  1526142908234895159 ,  1185903712917801891 ,  2095637833179633450 ,  1684061887934064807 ,  1459637973270975973 ,  2080524575243813761 ,  948583036372588019 ,  1196499685761340006 ,  1838560132753716756 ,  1035066941117701308 ,  1633636158414982077 ,  2082933534230116545 ,  161677865488580416 ,  936609732430805584 }, // x^2^605
+
+{ 409257345609329461 ,  1002655167281875678 ,  18470199587775000 ,  2184022032975947208 ,  148234944683729763 ,  2162572938835480438 ,  1216127754790813623 ,  1303762224282752882 ,  2064136757082740176 ,  518955562885803190 ,  1212819681300030394 ,  486562075491162613 ,  1283290324562294532 ,  285537074615846742 ,  321224661976450451 ,  2089560323722387545 ,  238275251292120492 }, // x^2^606
+
+{ 2214799429964260023 ,  694747922604332887 ,  312465276764308387 ,  314561248153726635 ,  1765153767754864756 ,  1092506825381016367 ,  1475330654269567868 ,  1379961995474423138 ,  1380037933141733384 ,  767686246763038167 ,  2008926192530430712 ,  304944109134708306 ,  261183406799178779 ,  326684916374054817 ,  1565573428651338229 ,  1920871811367768904 ,  2089282340964983664 }, // x^2^607
+
+{ 1133198628406257078 ,  977759521511265809 ,  1677025924443289061 ,  1729641585066757341 ,  2008693257322528571 ,  1421311051680525446 ,  2030362323209952329 ,  1921889928974731930 ,  1093849407549420400 ,  1836335521578497369 ,  407868381069619222 ,  10208083363410986 ,  1550647306879485613 ,  1032274398117330349 ,  889271835949247699 ,  508828554264478304 ,  313811787637759042 }, // x^2^608
+
+{ 265345699270236054 ,  432199859212969150 ,  874558322199992458 ,  1131965915049654978 ,  298519106167257998 ,  10561017697598846 ,  1036255447357471766 ,  990617823928118898 ,  183634792806558724 ,  123066866696319933 ,  1303845119890513396 ,  1915539994626985074 ,  2074006245341867746 ,  442443726074289720 ,  1398504520177556356 ,  2177822133301482183 ,  1618315631184181685 }, // x^2^609
+
+{ 2281731104267965972 ,  1565833165249376146 ,  2084133230658243672 ,  1466077014346623711 ,  2100113102745347532 ,  1509772133714278653 ,  1072220688498313491 ,  698970180110383152 ,  1702069119117461364 ,  786447122909266339 ,  1700725262288581680 ,  1201605067747742747 ,  1285909410275446050 ,  465944984489850383 ,  1581258227975449663 ,  1468892983675794809 ,  698369184018114847 }, // x^2^610
+
+{ 721614432883106726 ,  1168570529527401436 ,  1558741633956623588 ,  625821206916677993 ,  1729484582639558993 ,  202818544256349027 ,  942584755749892262 ,  384853745237896372 ,  1878778990418394018 ,  2006136039057883787 ,  2000492489726132701 ,  1628427185964600209 ,  456608645446900334 ,  559232459658214239 ,  1425457030239632449 ,  287429815653732691 ,  1537646113987862352 }, // x^2^611
+
+{ 168932144896498705 ,  1584765261048819233 ,  505170343346470483 ,  1661727371938733976 ,  58225995502531476 ,  1069662977595621678 ,  1534921440358104004 ,  1665430445118051718 ,  1421063952495178964 ,  1964407550412120431 ,  271153105959670492 ,  1303954744475330405 ,  1402417619008570585 ,  2067748336717081462 ,  635378162740815380 ,  1966388298672109925 ,  2276135507117268900 }, // x^2^612
+
+{ 209516738545946051 ,  126364818200849051 ,  1934171883678238548 ,  964383175854976730 ,  1426910927989095855 ,  1225737711393069252 ,  545817142750996161 ,  1941470598322211846 ,  1888603351110522964 ,  1693849545801590503 ,  1145960696683862080 ,  1544154409134561092 ,  2057079562609173728 ,  333449282037994372 ,  2138527354092198252 ,  1857839034788608299 ,  96149649256385712 }, // x^2^613
+
+{ 1680829537551411728 ,  2134567358027244766 ,  2233241443836144884 ,  882397883456946081 ,  1272913220199061746 ,  7737651019535212 ,  554065101789910033 ,  1237958145783195446 ,  1349598693130928045 ,  1212195573003581110 ,  1392362203075907328 ,  913027022784782751 ,  411026031675105306 ,  872967435734543202 ,  294712992119516029 ,  272260211295161762 ,  1636719289988127499 }, // x^2^614
+
+{ 1664102028071512887 ,  1120090114167125698 ,  678215188753805783 ,  2138907996085142066 ,  1721907479145867302 ,  1086048121088426493 ,  282597371457515207 ,  2252562071097166874 ,  880501568810700587 ,  1513947683103494818 ,  118765667707656670 ,  133926112036721838 ,  567919878484948061 ,  79410381524546935 ,  1454226739666450321 ,  1554107108291191805 ,  1442417286842488623 }, // x^2^615
+
+{ 1332050352400493680 ,  1321101334747797025 ,  124286076382368784 ,  2016174712781756901 ,  1731820299459346863 ,  545006347537545995 ,  1068874777126444609 ,  767902913608305873 ,  2267594054719694251 ,  1590335123866293983 ,  868156357394677648 ,  2249143243430804218 ,  347528411881171670 ,  1754162230470399965 ,  2070435897818234391 ,  1115080034382184583 ,  727133352156610422 }, // x^2^616
+
+{ 355843900790907311 ,  363191259548386230 ,  1247995738147301243 ,  1637822561418229980 ,  2161312386590710193 ,  827479641142168242 ,  551490280055656820 ,  116433208177525158 ,  219743626260276428 ,  1031342298455165770 ,  956941458946607562 ,  1971095287608216767 ,  1197895453910374577 ,  758722071093994752 ,  777125942555276007 ,  1415023106509593626 ,  318953060159670820 }, // x^2^617
+
+{ 2116009462654498660 ,  816261448355929278 ,  2272633719403001435 ,  573255612065573253 ,  968141095699581750 ,  1763511366639175689 ,  1173137617258455966 ,  481749672060198864 ,  750918552038818940 ,  519539565364235474 ,  1926446971667744959 ,  313336000413603815 ,  515893570456318246 ,  1497339799087374476 ,  478287049521531475 ,  1236625257683288385 ,  359335556903155074 }, // x^2^618
+
+{ 284411256262633501 ,  1789336543617378723 ,  1511246373802383265 ,  937616018328420337 ,  138227234082520224 ,  1667286674451103850 ,  693601919163926025 ,  13505785652584434 ,  720941360686097746 ,  958299499540120680 ,  2165799682030822154 ,  569950597702024056 ,  655621404198221825 ,  1755168971312171278 ,  2043596507506252973 ,  723457099999134940 ,  542348768717201728 }, // x^2^619
+
+{ 14465848372683520 ,  998603764202551247 ,  724938950471219595 ,  325209651844259777 ,  821482747183096187 ,  2233461975545399202 ,  1931346739821426382 ,  2293325179294131554 ,  243302969523090162 ,  1510343779638164774 ,  959965695881709576 ,  1831708054543131034 ,  883847853251423712 ,  1165812877635303160 ,  2007671568816070608 ,  44906786919162811 ,  1277515245984716365 }, // x^2^620
+
+{ 631454168243258384 ,  1518494517880443027 ,  2287322260213126417 ,  73685056254736833 ,  942598569327822559 ,  1532992770144658039 ,  2229403027115926821 ,  842743413258963958 ,  1865437000699847820 ,  369346269649870812 ,  1254468168457467223 ,  559290945795297932 ,  1457535932113599162 ,  1926542806880420086 ,  437581909774872900 ,  1002726378631501533 ,  182139319307728305 }, // x^2^621
+
+{ 1539785164984779758 ,  2045565025200148915 ,  120085588989386291 ,  1851688824536244542 ,  1980194980609785714 ,  1768302648237998273 ,  937420755919973678 ,  2210791368804399387 ,  328708870574712683 ,  1007467584698506438 ,  928652866696652949 ,  2228808519508067530 ,  1642940388306822455 ,  1618552641970232946 ,  108597732587521580 ,  1988412225258659433 ,  1660448156453751001 }, // x^2^622
+
+{ 131526482323749840 ,  503434317381483751 ,  1560765838807596555 ,  1027232124412255852 ,  654779282895217397 ,  2273302731188942421 ,  1551720915692515402 ,  1826828338382462468 ,  604485879763220723 ,  1057737165588422306 ,  106531819221246981 ,  425783162197347937 ,  489984741847314768 ,  523779814301405487 ,  1865105780134618992 ,  1677103130386412459 ,  318671848522267610 }, // x^2^623
+
+{ 1639089103549084783 ,  2260330277551218496 ,  754892153485552977 ,  2047710903557972804 ,  2229320667962732246 ,  307697180850646681 ,  423771572510991367 ,  17754501861697173 ,  1238910344976760345 ,  339140894113319789 ,  1445828176546996833 ,  2042997377381022070 ,  2207437390435628445 ,  551144982338719121 ,  384944764939825909 ,  1474198471229196249 ,  1659242296891127810 }, // x^2^624
+
+{ 795309147716951006 ,  1452711178054134843 ,  1664782861363981358 ,  1006405310566793269 ,  971226122955741921 ,  752466486101956684 ,  1168643572163743942 ,  1369730711422397346 ,  627654132922219702 ,  2104230754718381006 ,  1588968809641801416 ,  440238626377162178 ,  732134947635506053 ,  1160506871469695759 ,  395151165185992129 ,  957703098223865589 ,  1844022625316455575 }, // x^2^625
+
+{ 1964890520513046164 ,  829237203674997277 ,  1045128031232022678 ,  930832386865989645 ,  501678256081635892 ,  917829023495631073 ,  881831635034313098 ,  394230956869842550 ,  71895808688446751 ,  1874849501133161758 ,  645524260474403693 ,  301748792077270066 ,  1063834409487547589 ,  521966335959428443 ,  735272076777014321 ,  753999758410424987 ,  292792781171616060 }, // x^2^626
+
+{ 1976769506118304887 ,  328773719961068411 ,  538585052711106822 ,  2240145575420715174 ,  1291008337036982313 ,  1034350880207773432 ,  1888353090737342761 ,  1978548213503129609 ,  743903982747501766 ,  2240529593386457105 ,  229308641102639039 ,  1678061221588734337 ,  1248567483083113072 ,  1768596733943584627 ,  1268676342298500384 ,  216631313378191648 ,  1830565281699080858 }, // x^2^627
+
+{ 222335339024805767 ,  146026382943095887 ,  1745991091656491682 ,  364572344137903996 ,  491711339269051713 ,  1355802228883764600 ,  2287182613212117809 ,  1628584733201147266 ,  1730760626189599809 ,  1568968802576521126 ,  311144633575093285 ,  1257985237943883059 ,  2071814648333720000 ,  687885153740420654 ,  134464058021476074 ,  1159135545218675174 ,  1046016356313666266 }, // x^2^628
+
+{ 2218238585802550263 ,  1796073488881922547 ,  108147160614248440 ,  2147497677742828723 ,  940377749228338777 ,  166702239352185028 ,  14837186174842935 ,  1881782378259233528 ,  442792017566957734 ,  1550966771279227259 ,  1902483003501989215 ,  660158154476503000 ,  673105302791153020 ,  749237522153030898 ,  904762939332248683 ,  1085180116177134000 ,  32410680807805190 }, // x^2^629
+
+{ 1565357784323540059 ,  2004567400386087250 ,  563556471193393347 ,  2176504294263508416 ,  1354931238262729258 ,  1868302215996805633 ,  1039100426489789297 ,  540592058162156675 ,  2043699722662531658 ,  1444431990121166178 ,  931232743518887886 ,  2055389098129788451 ,  1235540167366576898 ,  1090861860232324736 ,  1732281293507875925 ,  1285969498415451875 ,  1757636524168521125 }, // x^2^630
+
+{ 2235287391127506303 ,  1401636021589024058 ,  2051095569701746947 ,  644542369524563085 ,  1717409066688580033 ,  839878318363782234 ,  964338104498650467 ,  193111125860611415 ,  41226228181058700 ,  659802192172642300 ,  481644033258717967 ,  1649597320109931775 ,  1169882823869469374 ,  2274302558105332483 ,  318951044439334783 ,  1790653470898866515 ,  1273274762300650850 }, // x^2^631
+
+{ 300315997535491308 ,  2172766343039137706 ,  1640387745692897945 ,  1527918068635484048 ,  882433389130238556 ,  519571876740982088 ,  840810181955840317 ,  967721991925117462 ,  757987265005423150 ,  692837321139327525 ,  1347675233763916865 ,  1262263002781039893 ,  366275012213686434 ,  1720518743710286974 ,  788123176573078844 ,  118677909193430055 ,  331168466283501375 }, // x^2^632
+
+{ 2204724021283944336 ,  614667311623152683 ,  1844295995732490949 ,  810904822281747865 ,  1057803411491878257 ,  2192784617269518738 ,  525328525709235908 ,  1625675926696707314 ,  2219520795837515346 ,  1111578739349697071 ,  850264822465802254 ,  1875225445537474991 ,  447934460789202298 ,  1577693216054964532 ,  2069208581718856575 ,  1844971513814599533 ,  1100047974003437056 }, // x^2^633
+
+{ 212063919658939349 ,  1955050766498331130 ,  1016643445611085361 ,  310923000117302398 ,  1108559594383531725 ,  1574096462279058604 ,  2058721019806527059 ,  1495555971494333315 ,  1290148625133033915 ,  1726857770258998357 ,  839176472902748525 ,  1938011636448040623 ,  1509176691584187835 ,  1684996881754410913 ,  339436873350779373 ,  1476161426144378753 ,  1443502235558451143 }, // x^2^634
+
+{ 2154146131717687018 ,  1161290317713326050 ,  362912384144320809 ,  1110332862415808487 ,  716729690755657231 ,  1048283751964241153 ,  1336749180461386759 ,  2280348028206541371 ,  1259918326978354156 ,  1980538993489881425 ,  1333565492277683625 ,  514625871678292821 ,  1173265787958600132 ,  569893472219264300 ,  812175666077334991 ,  1255701238860433930 ,  79993120593823961 }, // x^2^635
+
+{ 409217992666385378 ,  728224420299246128 ,  1959936179576785895 ,  1444295307550571182 ,  998063979246407739 ,  554682441648959407 ,  2073233365208088525 ,  1349861235690939694 ,  1097281043663983707 ,  1649731108960027223 ,  2214398868943408274 ,  1449553717147663103 ,  2196500661345051527 ,  2171871885380914429 ,  491427195052411299 ,  470657045051378516 ,  554494320060727622 }, // x^2^636
+
+{ 2211653713656632985 ,  1473087100459732074 ,  345838890415354412 ,  1445623252958685565 ,  1041419654691395738 ,  2188283266242968331 ,  1417777133392604784 ,  1154814074996216406 ,  1442989098963226134 ,  901383116742510681 ,  497212216058355345 ,  999731586376214828 ,  817831951035275271 ,  159325628653133022 ,  120914548489015393 ,  1449314802372642427 ,  471104110110980884 }, // x^2^637
+
+{ 323038291515073224 ,  265070044955668299 ,  1053676807411822176 ,  2247226908652601286 ,  2083016828474103062 ,  273234004102527575 ,  11857926472424922 ,  2055314572998661967 ,  311921811886787132 ,  1870415531374949840 ,  21265502209057804 ,  678112836970959466 ,  926490438737525589 ,  2027644346888215584 ,  495557988374717777 ,  1817993454480752541 ,  2141242677896973878 }, // x^2^638
+
+{ 1445163042542523231 ,  1767410038768111356 ,  462427149761798838 ,  1007504177155879483 ,  1181224656949316486 ,  357022206754812124 ,  1898271363493710619 ,  1339720449183752193 ,  883758224923786172 ,  36020177946849132 ,  461791639260431974 ,  1411405712850290498 ,  701758173778161093 ,  1222887871945629218 ,  1174213062638935911 ,  326342461405188581 ,  2211032895429848005 } // x^2^639
+
+};
diff --git a/madspace/include/madspace/phasespace/channel_weight_network.hpp b/madspace/include/madspace/phasespace/channel_weight_network.hpp
index b7e0cbb0e..97802381e 100644
--- a/madspace/include/madspace/phasespace/channel_weight_network.hpp
+++ b/madspace/include/madspace/phasespace/channel_weight_network.hpp
@@ -32,7 +32,10 @@ class ChannelWeightNetwork : public FunctionGenerator {
 
     const MLP& mlp() const { return _mlp; }
     const MomentumPreprocessing& preprocessing() const { return _preprocessing; }
-    void initialize_globals(ContextPtr context) const;
+    // nullopt (the default) initializes non-deterministically.
+    void initialize_globals(
+        ContextPtr context, std::optional<std::uint64_t> seed = std::nullopt
+    ) const;
     const std::string& mask_name() const { return _mask_name; }
 
 private:
diff --git a/madspace/include/madspace/phasespace/discrete_flow.hpp b/madspace/include/madspace/phasespace/discrete_flow.hpp
index 5dc05a793..88913bfa5 100644
--- a/madspace/include/madspace/phasespace/discrete_flow.hpp
+++ b/madspace/include/madspace/phasespace/discrete_flow.hpp
@@ -18,7 +18,10 @@ class DiscreteFlow : public Mapping {
     );
     const std::vector<std::size_t>& option_counts() const { return _option_counts; }
     std::size_t condition_dim() const { return _condition_dim; }
-    void initialize_globals(ContextPtr context) const;
+    // nullopt (the default) initializes non-deterministically.
+    void initialize_globals(
+        ContextPtr context, std::optional<std::uint64_t> seed = std::nullopt
+    ) const;
 
 private:
     Result build_forward_impl(
diff --git a/madspace/include/madspace/phasespace/flow.hpp b/madspace/include/madspace/phasespace/flow.hpp
index fcf064084..51149d6aa 100644
--- a/madspace/include/madspace/phasespace/flow.hpp
+++ b/madspace/include/madspace/phasespace/flow.hpp
@@ -1,5 +1,7 @@
 #pragma once
 
+#include <optional>
+
 #include "madspace/phasespace/base.hpp"
 #include "madspace/phasespace/mlp.hpp"
 
@@ -19,8 +21,15 @@ class Flow : public Mapping {
     );
     std::size_t input_dim() const { return _input_dim; }
     std::size_t condition_dim() const { return _condition_dim; }
-    void initialize_globals(ContextPtr context) const;
-    void initialize_from_vegas(ContextPtr context, const std::string& grid_name) const;
+    // nullopt (the default) initializes non-deterministically.
+    void initialize_globals(
+        ContextPtr context, std::optional<std::uint64_t> seed = std::nullopt
+    ) const;
+    void initialize_from_vegas(
+        ContextPtr context,
+        const std::string& grid_name,
+        std::optional<std::uint64_t> seed = std::nullopt
+    ) const;
 
 private:
     Result build_forward_impl(
diff --git a/madspace/include/madspace/phasespace/mlp.hpp b/madspace/include/madspace/phasespace/mlp.hpp
index 45234a8f3..5b6b374b6 100644
--- a/madspace/include/madspace/phasespace/mlp.hpp
+++ b/madspace/include/madspace/phasespace/mlp.hpp
@@ -1,5 +1,8 @@
 #pragma once
 
+#include <cstdint>
+#include <optional>
+
 #include "madspace/driver/context.hpp"
 #include "madspace/phasespace/base.hpp"
 
@@ -17,7 +20,10 @@ class MLP : public FunctionGenerator {
 
     std::size_t input_dim() const { return _input_dim; }
     std::size_t output_dim() const { return _output_dim; }
-    void initialize_globals(ContextPtr context) const;
+    // nullopt (the default) initializes non-deterministically.
+    void initialize_globals(
+        ContextPtr context, std::optional<std::uint64_t> seed = std::nullopt
+    ) const;
     std::string last_layer_bias_name() const {
         return prefixed_name(_prefix, std::format("layer{}.bias", _layers));
     }
diff --git a/madspace/include/madspace/util.hpp b/madspace/include/madspace/util.hpp
index 7b1e43b38..de5d8c9b5 100644
--- a/madspace/include/madspace/util.hpp
+++ b/madspace/include/madspace/util.hpp
@@ -2,6 +2,7 @@
 
 #include <cstdio>
 #include <format>
+#include <initializer_list>
 #include <ranges>
 #include <tuple>
 #include <unordered_map>
diff --git a/madspace/src/compgraphs/function.cpp b/madspace/src/compgraphs/function.cpp
index 6552c0327..cc6ce199d 100644
--- a/madspace/src/compgraphs/function.cpp
+++ b/madspace/src/compgraphs/function.cpp
@@ -361,10 +361,15 @@ FunctionBuilder::instruction(InstructionPtr instruction, const ValueVec& args) {
         }
     }
 
+    // random instructions share one GPU RNG stream, so they always run on the
+    // main stream (0) and are never CSE'd (see Instruction::is_random)
+    bool is_random = instruction->is_random();
+    std::size_t call_stream = is_random ? 0 : _current_stream;
+
     // create local variables for constants
     std::vector<std::size_t> cache_key;
     cache_key.push_back(opcode);
-    cache_key.push_back(_current_stream);
+    cache_key.push_back(call_stream);
     for (auto& arg : params) {
         register_local(arg);
         cache_key.push_back(arg.local_index);
@@ -381,8 +386,8 @@ FunctionBuilder::instruction(InstructionPtr instruction, const ValueVec& args) {
         }
     }
 
-    // check for cached result (for deterministic instructions)
-    if (opcode != opcodes::random && opcode != opcodes::unweight) {
+    // check for cached result
+    if (!is_random) {
         auto find_instr = _instruction_cache.find(cache_key);
         if (find_instr != _instruction_cache.end()) {
             ValueVec call_outputs;
@@ -405,10 +410,12 @@ FunctionBuilder::instruction(InstructionPtr instruction, const ValueVec& args) {
         output_locals.push_back(value.local_index);
     }
     _instructions.push_back(
-        InstructionCall{instruction, params, call_outputs, _current_stream}
+        InstructionCall{instruction, params, call_outputs, call_stream}
     );
     _instruction_use_count.push_back(0);
-    _instruction_cache[cache_key] = output_locals;
+    if (!is_random) {
+        _instruction_cache[cache_key] = output_locals;
+    }
     for (auto& param : params) {
         int param_source = _local_sources.at(param.local_index);
         if (param_source != -1) {
diff --git a/madspace/src/cpu/device.hpp b/madspace/src/cpu/device.hpp
index a29bb9e5b..06184ce67 100644
--- a/madspace/src/cpu/device.hpp
+++ b/madspace/src/cpu/device.hpp
@@ -1,5 +1,9 @@
 #pragma once
 
+#include <cstdint>
+#include <random>
+
+#include "madspace/driver/random.hpp"
 #include "madspace/driver/tensor.hpp"
 #include "madspace/driver/thread_pool.hpp"
 #include "simd.hpp"
@@ -7,6 +11,8 @@
 namespace madspace {
 namespace cpu {
 
+class CpuRuntime;
+
 class CpuDevice : public Device {
 public:
     static constexpr bool is_concurrent = false;
diff --git a/madspace/src/cpu/runtime.cpp b/madspace/src/cpu/runtime.cpp
index 02ed5a5be..99c0eb136 100644
--- a/madspace/src/cpu/runtime.cpp
+++ b/madspace/src/cpu/runtime.cpp
@@ -1,6 +1,7 @@
 #include "runtime.hpp"
 
 #include <algorithm>
+#include <memory>
 #include <random>
 #include <ranges>
 #include <tuple>
@@ -613,12 +614,11 @@ void op_random(
     auto& runtime = instruction.runtime;
     device.foreach (
         flat_view.shape[0],
-        [flat_view, &runtime](std::size_t count, std::size_t offset) mutable {
+        [flat_view, &runtime, &device](std::size_t count, std::size_t offset) mutable {
             auto output_view = TensorView<double, 1>(flat_view);
-            std::uniform_real_distribution<double> dist;
             auto& rand_gen = runtime.rand_gen();
             for (std::size_t i = offset; i < offset + count; ++i) {
-                output_view[i] = dist(rand_gen);
+                output_view[i] = rand_gen.generate_double();
             }
         }
     );
@@ -636,12 +636,13 @@ void op_random_int(
     auto& runtime = instruction.runtime;
     device.foreach (
         flat_view.shape[0],
-        [flat_view, max_val, &runtime](std::size_t count, std::size_t offset) mutable {
+        [flat_view, max_val, &runtime, &device](
+            std::size_t count, std::size_t offset
+        ) mutable {
             auto output_view = TensorView<me_int_t, 1>(flat_view);
-            std::uniform_int_distribution<me_int_t> dist(0, max_val - 1);
             auto& rand_gen = runtime.rand_gen();
             for (std::size_t i = offset; i < offset + count; ++i) {
-                output_view[i] = dist(rand_gen);
+                output_view[i] = rand_gen.generate_int(max_val);
             }
         }
     );
@@ -673,6 +674,7 @@ void op_unweight(
          indices_view_flat,
          uw_weights_view_flat,
          &runtime,
+         &device,
          batch_size,
          &indices,
          &uw_weights,
@@ -682,12 +684,11 @@ void op_unweight(
             TensorView<double, 1> max_weight_view(max_weight_view_flat);
             TensorView<me_int_t, 1> indices_view(indices_view_flat);
             TensorView<double, 1> uw_weights_view(uw_weights_view_flat);
-            std::uniform_real_distribution<double> dist;
             auto& rand_gen = runtime.rand_gen();
             std::size_t count = 0;
             for (std::size_t i = 0; i < batch_size; ++i) {
                 double w = weights_view[i], w_max = max_weight_view[i];
-                if (w != 0. && w > dist(rand_gen) * w_max) {
+                if (w != 0. && w > rand_gen.generate_double() * w_max) {
                     indices_view[count] = i;
                     uw_weights_view[count] = w > w_max ? w : w_max;
                     ++count;
@@ -883,13 +884,7 @@ void op_discrete_histogram(
 CpuRuntime::CpuRuntime(const Function& function, ContextPtr context, bool concurrent) :
     _context(context),
     _input_count(function.inputs().size()),
-    _rand_gens(
-        context->thread_pool(),
-        []() {
-            std::random_device rand_device;
-            return std::mt19937(rand_device());
-        }
-    ),
+    _rand_gens(context->thread_pool(), []() { return MixMaxRandom(); }),
     _concurrent(concurrent) {
     if (context->device()->device_type() != DeviceType::cpu) {
         throw std::runtime_error("Context has incompatible device");
diff --git a/madspace/src/cpu/runtime.hpp b/madspace/src/cpu/runtime.hpp
index 9b70fdbfc..fd18a0462 100644
--- a/madspace/src/cpu/runtime.hpp
+++ b/madspace/src/cpu/runtime.hpp
@@ -1,10 +1,12 @@
 #pragma once
 
+#include <optional>
 #include <random>
 
 #include "madspace/compgraphs/function.hpp"
 #include "madspace/driver/backend.hpp"
 #include "madspace/driver/context.hpp"
+#include "madspace/driver/random.hpp"
 #include "madspace/driver/tensor.hpp"
 
 namespace madspace {
@@ -49,7 +51,8 @@ class CpuRuntime : public Runtime {
     ) override;
 
     Context& context() { return *_context; }
-    std::mt19937& rand_gen() { return _rand_gens.get(); }
+    MixMaxRandom& rand_gen() { return _rand_gens.get(); }
+    void set_seed(DerivedSeed seed) override { rand_gen().set_seed(seed); }
 
 private:
     TensorVec run_single(const TensorVec& inputs) const;
@@ -83,7 +86,7 @@ class CpuRuntime : public Runtime {
     std::vector<Sizes> _grad_global_shapes;
     std::size_t _grad_global_total_size;
     ContextPtr _context;
-    ThreadResource<std::mt19937> _rand_gens;
+    ThreadResource<MixMaxRandom> _rand_gens;
     bool _concurrent;
     SizeVec _ready_instructions_init;
     SizeVec _ready_instructions_backward_init;
diff --git a/madspace/src/cpu/tensor.hpp b/madspace/src/cpu/tensor.hpp
index d4397ef40..cb0bd7afe 100644
--- a/madspace/src/cpu/tensor.hpp
+++ b/madspace/src/cpu/tensor.hpp
@@ -305,6 +305,20 @@ inline void tensor_foreach_impl(
     bool single_job,
     S... scalar_args
 ) {
+    // A broadcast *output* (size(0) != batch_size, e.g. a global's gradient
+    // accumulator during backward) is shared across the whole batch instead of
+    // having one slot per batch element. Splitting the batch across worker threads
+    // would then have multiple threads accumulate into the very same memory
+    // concurrently with no synchronization, so force a single job in that case
+    // (the whole batch handled sequentially by one thread). Broadcast *inputs* are
+    // only ever read, never written, in both forward and backward -- concurrent
+    // reads of the same memory are safe, so they don't need this.
+    for (const Tensor* out : outputs) {
+        if (out->size(0) != batch_size) {
+            single_job = true;
+        }
+    }
+
     // get views to the tensors with the correct types based on the signature of
     // scalar_func
     auto flat_views = std::apply(
diff --git a/madspace/src/driver/channel_generator.cpp b/madspace/src/driver/channel_generator.cpp
index c86da5107..c4e140f97 100644
--- a/madspace/src/driver/channel_generator.cpp
+++ b/madspace/src/driver/channel_generator.cpp
@@ -1,4 +1,6 @@
 #include "madspace/driver/channel_generator.hpp"
+#include "madspace/driver/random.hpp"
+#include "madspace/util.hpp"
 
 using namespace madspace;
 using json = nlohmann::json;
@@ -29,6 +31,44 @@ int particle_extra_flags(
     return flags;
 }
 
+DerivedSeed::SeedType generate_seed_type(bool is_survey, std::size_t survey_pass) {
+    if (!is_survey) {
+        return DerivedSeed::generator_generate;
+    }
+    return survey_pass == 0
+        ? DerivedSeed::first_survey_generate
+        : DerivedSeed::second_survey_generate;
+}
+
+DerivedSeed::SeedType unweight_seed_type(bool is_survey, std::size_t survey_pass) {
+    if (!is_survey) {
+        return DerivedSeed::generator_unweight;
+    }
+    return survey_pass == 0
+        ? DerivedSeed::first_survey_unweight
+        : DerivedSeed::second_survey_unweight;
+}
+
+// Role folded into job_index so integrand_channel and integrand_common get
+// independent streams despite sharing seed_type/channel_index.
+enum RngRole { rng_role_channel = 0, rng_role_common = 1, rng_role_count = 2 };
+
+DerivedSeed generate_phase_seed(
+    std::optional<std::uint64_t> seed,
+    bool is_survey,
+    std::size_t survey_pass,
+    std::size_t job_index,
+    std::size_t channel_index,
+    RngRole role
+) {
+    return DerivedSeed(
+        seed,
+        generate_seed_type(is_survey, survey_pass),
+        job_index * rng_role_count + role,
+        channel_index
+    );
+}
+
 } // namespace
 
 ChannelEventGenerator::ChannelEventGenerator(
@@ -53,7 +93,7 @@ ChannelEventGenerator::ChannelEventGenerator(
         .count_after_cuts = 0,
         .count_after_cuts_opt = 0,
         .count_unweighted = 0.,
-        .count_target = 1.,
+        .count_target = 1,
         .optimized = false,
         .done = false
     },
@@ -177,7 +217,7 @@ ChannelEventGenerator::ChannelEventGenerator(
         .count_after_cuts = 0,
         .count_after_cuts_opt = 0,
         .count_unweighted = 0.,
-        .count_target = 1.,
+        .count_target = 1,
         .optimized = false,
         .done = false
     },
@@ -275,18 +315,17 @@ void ChannelEventGenerator::init_field_indices() {
     _field_indices.rest = _field_indices.random + 1;
 }
 
-void ChannelEventGenerator::unweight_file(std::mt19937& rand_gen) {
+void ChannelEventGenerator::unweight_file(MixMaxRandom& rand_gen) {
     std::size_t buf_size = 1000000;
-    std::uniform_real_distribution<double> rand_dist;
     EventBuffer buffer(0, 0, weight_file_layout);
-    std::size_t accept_count = _unweighted_count;
+    std::size_t accept_count = _unweighted_accept_count;
     for (std::size_t i = _unweighted_count; i < _weight_file.event_count();
          i += buf_size) {
         _weight_file.seek(i);
         _weight_file.read(buffer, buf_size);
         for (std::size_t j = 0; j < buffer.event_count(); ++j) {
             auto weight = buffer.event(j).weight();
-            if (weight / _max_weight < rand_dist(rand_gen)) {
+            if (weight / _max_weight < rand_gen.generate_double()) {
                 weight = 0;
             } else {
                 weight = std::max(weight.value(), _max_weight);
@@ -296,6 +335,8 @@ void ChannelEventGenerator::unweight_file(std::mt19937& rand_gen) {
         _weight_file.seek(i);
         _weight_file.write(buffer);
     }
+    _unweighted_count = _weight_file.event_count();
+    _unweighted_accept_count = accept_count;
     _status.count_unweighted = accept_count;
 }
 
@@ -387,14 +428,35 @@ double ChannelEventGenerator::channel_weight_sum(std::size_t event_count) {
     return weight_sum;
 }
 
+// Assigns job's base seed on the scheduling thread, from its logical identity
+// only (channel, kind, sequence) -- never from context_index -- so it doesn't
+// depend on worker/context assignment.
 void ChannelEventGenerator::start_job(
-    GeneratorBatchJob& job, ResultQueue& result_queue
+    GeneratorBatchJob& job,
+    ResultQueue& result_queue,
+    std::optional<std::uint64_t> seed,
+    bool is_survey,
+    std::size_t survey_pass
 ) {
+    job.rng_seed = seed;
+    job.rng_is_survey = is_survey;
+    job.rng_survey_pass = survey_pass;
+    job.rng_job_index = is_survey ? _survey_rng_seq++ : _generate_rng_seq++;
     _contexts.at(job.context_index)
         ->thread_pool()
         .submit([this, &job, &result_queue]() {
             auto& runtimes = _runtimes.at(job.context_index);
             auto& context = _contexts.at(job.context_index);
+            if (job.rng_seed) {
+                runtimes.integrand_channel->set_seed(generate_phase_seed(
+                    job.rng_seed,
+                    job.rng_is_survey,
+                    job.rng_survey_pass,
+                    job.rng_job_index,
+                    job.channel_index,
+                    rng_role_channel
+                ));
+            }
             std::size_t max_batch_size =
                 context->device()->device_type() == DeviceType::cpu
                 ? _config.cpu_batch_size
@@ -408,6 +470,8 @@ void ChannelEventGenerator::start_job(
             std::size_t total_count = 0, repetitions = 0;
             TensorVec all_ps_points;
             while (true) {
+                // Successive calls draw further from integrand_channel's own
+                // persistent stream (set once above) -- no per-call re-seeding.
                 auto ps_points =
                     runtimes.integrand_channel->run({Tensor({batch_size})});
                 std::size_t acc_count =
@@ -444,9 +508,21 @@ void ChannelEventGenerator::start_job(
                     (target_count - total_count) / cut_eff
                 );
             }
+            if (job.rng_seed) {
+                runtimes.integrand_common->set_seed(generate_phase_seed(
+                    job.rng_seed,
+                    job.rng_is_survey,
+                    job.rng_survey_pass,
+                    job.rng_job_index,
+                    job.channel_index,
+                    rng_role_common
+                ));
+            }
             job.events = runtimes.integrand_common->run(all_ps_points);
 
             job.weights = job.events.at(_field_indices.weight).cpu();
+            // observable_histograms/vegas_histogram/discrete_histogram don't consume
+            // random numbers, so they're never seeded.
             if (runtimes.observable_histograms) {
                 auto hists = runtimes.observable_histograms->run(
                     {job.events.at(_field_indices.weight),
@@ -456,7 +532,7 @@ void ChannelEventGenerator::start_job(
                     job.hists.push_back(item.cpu());
                 }
             }
-            if (job.vegas_batch_size != 0) {
+            if (job.is_vegas_batch) {
                 if (_vegas_optimizer) {
                     auto hist = runtimes.vegas_histogram->run(
                         {job.events.at(_field_indices.random),
@@ -482,15 +558,26 @@ void ChannelEventGenerator::start_job(
         });
 }
 
-void ChannelEventGenerator::start_unweight_job(
+void ChannelEventGenerator::prepare_unweight_job(GeneratorBatchJob& job) const {
+    job.max_weight = _max_weight;
+}
+
+void ChannelEventGenerator::submit_unweight_job(
     GeneratorBatchJob& job, ResultQueue& result_queue
 ) {
-    job.max_weight = _max_weight;
     _contexts.at(job.context_index)
         ->thread_pool()
         .submit([this, &job, &result_queue]() {
             auto& runtimes = _runtimes.at(job.context_index);
             auto& context = _contexts.at(job.context_index);
+            if (job.rng_seed) {
+                runtimes.unweighter->set_seed(DerivedSeed(
+                    job.rng_seed,
+                    unweight_seed_type(job.rng_is_survey, job.rng_survey_pass),
+                    job.rng_job_index,
+                    job.channel_index
+                ));
+            }
             std::vector<Tensor> unweighter_args(
                 job.events.begin(), job.events.begin() + _field_indices.random
             );
@@ -504,6 +591,13 @@ void ChannelEventGenerator::start_unweight_job(
         });
 }
 
+void ChannelEventGenerator::start_unweight_job(
+    GeneratorBatchJob& job, ResultQueue& result_queue
+) {
+    prepare_unweight_job(job);
+    submit_unweight_job(job, result_queue);
+}
+
 std::size_t ChannelEventGenerator::next_vegas_batch_size() {
     std::size_t batch_size = _batch_size;
     _batch_size = std::min(_batch_size * 2, _config.max_batch_size);
@@ -514,6 +608,7 @@ void ChannelEventGenerator::clear_events() {
     _status.count_unweighted = 0;
     _max_weight = 0;
     _unweighted_count = 0;
+    _unweighted_accept_count = 0;
     _status.count_opt = 0;
     _status.count_after_cuts_opt = 0;
     _event_file.clear();
@@ -549,8 +644,9 @@ void ChannelEventGenerator::update_max_weight(Tensor weights) {
 
     double w_sum = 0;
     double max_truncation = _config.max_overweight_truncation *
-        std::min(_status.count_target,
-                 static_cast<double>(_config.freeze_max_weight_after));
+        static_cast<double>(std::min(
+            _status.count_target, _config.freeze_max_weight_after
+        ));
     std::size_t count = 0;
     for (auto w : _large_weights) {
         if (w < _max_weight) {
@@ -563,6 +659,7 @@ void ChannelEventGenerator::update_max_weight(Tensor weights) {
                 _status.count_unweighted *= _max_weight / w;
                 _max_weight = w;
                 _unweighted_count = 0;
+                _unweighted_accept_count = 0;
             }
             break;
         }
diff --git a/madspace/src/driver/event_generator.cpp b/madspace/src/driver/event_generator.cpp
index a5e4760a5..f8625371f 100644
--- a/madspace/src/driver/event_generator.cpp
+++ b/madspace/src/driver/event_generator.cpp
@@ -1,11 +1,15 @@
 #include "madspace/driver/event_generator.hpp"
 
+#include <algorithm>
 #include <cmath>
 #include <filesystem>
 #include <format>
+#include <memory>
+#include <numeric>
 #include <ranges>
 
 #include "madspace/driver/logger.hpp"
+#include "madspace/driver/random.hpp"
 #include "madspace/util.hpp"
 
 using namespace madspace;
@@ -16,7 +20,8 @@ EventGenerator::EventGenerator(
     const std::vector<ContextPtr>& contexts,
     const std::vector<std::shared_ptr<ChannelEventGenerator>>& channels,
     const std::string& status_file,
-    const GeneratorConfig& config
+    const GeneratorConfig& config,
+    std::optional<std::uint64_t> seed
 ) :
     _config(config),
     _status{
@@ -30,7 +35,7 @@ EventGenerator::EventGenerator(
         .count_after_cuts = 0,
         .count_after_cuts_opt = 0,
         .count_unweighted = 0.,
-        .count_target = static_cast<double>(config.target_count),
+        .count_target = config.target_count,
         .optimized = false,
         .done = false
     },
@@ -41,9 +46,24 @@ EventGenerator::EventGenerator(
     _channel_optimizing(channels.size()),
     _channel_integral_fractions(channels.size(), 1.),
     _context_job_counts(contexts.size()),
+    _channel_batch_pending(channels.size()),
+    _channel_ready_gen(channels.size()),
+    _channel_commit_cursor(channels.size()),
+    _channel_cursor_set(channels.size()),
+    _channel_unweight_ready(channels.size()),
+    _channel_unweight_cursor(channels.size()),
+    _context_unweight_queue(contexts.size()),
+    _seed(seed),
+    _deterministic(seed.has_value()),
     _status_file(status_file) {}
 
-void EventGenerator::survey() {
+void EventGenerator::survey(std::size_t survey_pass) {
+    _survey_job = true;
+    _survey_pass = survey_pass;
+    if (_deterministic) {
+        survey_deterministic();
+        return;
+    }
     reset_start_time();
     bool done = false;
     std::size_t min_iters = _config.survey_min_iters;
@@ -80,6 +100,7 @@ void EventGenerator::survey() {
                 .channel_index = i,
                 .unweight = iter >= min_iters - 1,
                 .vegas_batch_size = vegas_batch_size,
+                .is_vegas_batch = true,
             });
             if (iter >= min_iters) {
                 total_event_count += vegas_batch_size;
@@ -136,6 +157,11 @@ void EventGenerator::survey() {
 }
 
 void EventGenerator::generate() {
+    _survey_job = false;
+    if (_deterministic) {
+        generate_deterministic();
+        return;
+    }
     reset_start_time();
     print_gen_init();
 
@@ -168,6 +194,7 @@ void EventGenerator::generate() {
                             .channel_index = channel_index,
                             .unweight = true,
                             .vegas_batch_size = channel->next_vegas_batch_size(),
+                            .is_vegas_batch = true,
                         });
                     }
                 } else {
@@ -227,14 +254,326 @@ void EventGenerator::generate() {
     print_gen_update(true);
 }
 
-bool EventGenerator::start_jobs() {
+// Commits a job's generate stage in ascending job id per channel. Doesn't write
+// events itself -- if job.unweight, that's deferred to commit_unweight_job().
+void EventGenerator::commit_generate_job(GeneratorBatchJob& job) {
+    auto& channel = _channels.at(job.channel_index);
+    auto& channel_job_count = _channel_job_counts.at(job.channel_index);
+    if (job.is_vegas_batch && channel_job_count == job.split_job_count) {
+        channel->clear_events();
+    }
+    channel->integrate(job);
+    update_integral_status();
+    channel->update_max_weight(job.weights);
+    --channel_job_count;
+    --_context_job_counts.at(job.context_index);
+    if (job.unweight) {
+        // Snapshot max_weight here, at this job's fixed commit position.
+        channel->prepare_unweight_job(job);
+        _context_unweight_queue.at(job.context_index).push_back(job.job_id);
+    }
+    finish_channel_job(job);
+    print_gen_update(false);
+}
+
+void EventGenerator::commit_unweight_job(GeneratorBatchJob& job) {
+    auto& channel = _channels.at(job.channel_index);
+    channel->write_events(job.unweighted_events, job.max_weight);
+    update_counts();
+    --_channel_job_counts.at(job.channel_index);
+    --_context_job_counts.at(job.context_index);
+    finish_channel_job(job);
+    print_gen_update(false);
+}
+
+void EventGenerator::finish_channel_job(const GeneratorBatchJob& job) {
+    auto& channel_job_count = _channel_job_counts.at(job.channel_index);
+    if (channel_job_count != 0) {
+        return;
+    }
+    if (job.is_vegas_batch) {
+        _channels.at(job.channel_index)->optimize_vegas(job);
+        _channel_optimizing.at(job.channel_index) = false;
+    } else {
+        _channel_batch_pending.at(job.channel_index) = false;
+    }
+}
+
+void EventGenerator::register_dispatched_ids(std::size_t first_id, std::size_t end_id) {
+    for (std::size_t id = first_id; id < end_id; ++id) {
+        std::size_t channel_index = _running_jobs.at(id).channel_index;
+        if (!_channel_cursor_set.at(channel_index)) {
+            _channel_commit_cursor.at(channel_index) = id;
+            _channel_unweight_cursor.at(channel_index) = id;
+            _channel_cursor_set.at(channel_index) = true;
+        }
+    }
+}
+
+// Deterministic generate path: jobs still run on the pool, but commit strictly
+// in per-channel job-id order, so results are a pure function of the seed.
+void EventGenerator::generate_deterministic() {
+    reset_start_time();
+    print_gen_init();
+
+    while (true) {
+        _abort_check_function();
+
+        for (std::size_t channel_index = 0; channel_index < _channels.size();
+             ++channel_index) {
+            auto& channel = _channels.at(channel_index);
+            double integral_frac = _channel_integral_fractions.at(channel_index);
+            if (integral_frac > 0 &&
+                channel->status().count_unweighted >= channel->status().count_target) {
+                continue;
+            }
+            if (channel->needs_optimization()) {
+                if (!_channel_optimizing.at(channel_index)) {
+                    _channel_optimizing.at(channel_index) = true;
+                    _channel_cursor_set.at(channel_index) = false;
+                    _ready_jobs.push_back({
+                        .channel_index = channel_index,
+                        .unweight = true,
+                        .vegas_batch_size = channel->next_vegas_batch_size(),
+                        .is_vegas_batch = true,
+                    });
+                }
+            } else if (!_channel_batch_pending.at(channel_index)) {
+                _channel_batch_pending.at(channel_index) = true;
+                _channel_cursor_set.at(channel_index) = false;
+                std::size_t batch_jobs = next_batch_job_count(channel_index);
+                // Single ready_job with vegas_batch_size set so start_jobs() dispatches
+                // the whole batch atomically, giving the round a contiguous job-id
+                // range.
+                _ready_jobs.push_back({
+                    .channel_index = channel_index,
+                    .unweight = true,
+                    .vegas_batch_size = batch_jobs * _config.cpu_batch_size,
+                    .is_vegas_batch = false,
+                });
+            }
+        }
+
+        std::size_t job_id_before = _job_id;
+        std::size_t unweight_dispatched = start_jobs();
+        std::size_t round_in_flight = (_job_id - job_id_before) + unweight_dispatched;
+        register_dispatched_ids(job_id_before, _job_id);
+
+        // Wait for the round to fully commit before resyncing cross-channel fractions.
+        while (round_in_flight > 0) {
+            std::size_t job_id = _result_queue.wait();
+            --round_in_flight;
+            auto& job = _running_jobs.at(job_id);
+            if (job.unweighted_events.size() == 0) {
+                auto& ready = _channel_ready_gen.at(job.channel_index);
+                auto& cursor = _channel_commit_cursor.at(job.channel_index);
+                ready.insert(job_id);
+                while (ready.erase(cursor) > 0) {
+                    auto& gen_job = _running_jobs.at(cursor);
+                    commit_generate_job(gen_job);
+                    if (!gen_job.unweight) {
+                        _running_jobs.erase(cursor);
+                    }
+                    ++cursor;
+                }
+            } else {
+                auto& ready = _channel_unweight_ready.at(job.channel_index);
+                auto& cursor = _channel_unweight_cursor.at(job.channel_index);
+                ready.insert(job_id);
+                while (ready.erase(cursor) > 0) {
+                    commit_unweight_job(_running_jobs.at(cursor));
+                    _running_jobs.erase(cursor);
+                    ++cursor;
+                }
+            }
+            std::size_t job_id_refill = _job_id;
+            std::size_t refill_unweight_dispatched = start_jobs();
+            round_in_flight += (_job_id - job_id_refill) + refill_unweight_dispatched;
+            register_dispatched_ids(job_id_refill, _job_id);
+        }
+
+        update_integral_fractions();
+
+        if (_status.done) {
+            unweight_all();
+        }
+        if (_status.done) {
+            break;
+        }
+    }
+    print_gen_update(true);
+}
+
+// Deterministic survey path, same commit-ordering approach as generate_deterministic().
+void EventGenerator::survey_deterministic() {
+    reset_start_time();
+    _commit_cursor = _job_id;
+    _ready_gen.clear();
+    bool done = false;
+    std::size_t min_iters = _config.survey_min_iters;
+    std::size_t max_iters = std::max(min_iters, _config.survey_max_iters);
+    double target_precision = _config.survey_target_precision;
+
+    std::size_t total_event_count = 0;
+    std::size_t done_event_count = 0;
+    for (auto& channel : _channels) {
+        std::size_t chan_batch_size = channel->batch_size();
+        for (std::size_t iter = channel->status().iterations; iter < min_iters;
+             ++iter) {
+            total_event_count += chan_batch_size;
+            chan_batch_size = std::min(chan_batch_size * 2, _config.max_batch_size);
+        }
+    }
+    print_survey_init();
+
+    std::size_t iter = 0;
+    for (; !done && iter < max_iters; ++iter) {
+        for (std::size_t i = 0; auto channel : _channels) {
+            if (channel->status().iterations > iter) {
+                ++i;
+                continue;
+            }
+            if (iter >= min_iters &&
+                channel->cross_section().rel_error() < target_precision) {
+                ++i;
+                continue;
+            }
+            std::size_t vegas_batch_size = channel->next_vegas_batch_size();
+            // Let register_dispatched_ids() re-init the cursor for this new job.
+            _channel_cursor_set.at(i) = false;
+            _ready_jobs.push_back({
+                .channel_index = i,
+                .unweight = iter >= min_iters - 1,
+                .vegas_batch_size = vegas_batch_size,
+                .is_vegas_batch = true,
+            });
+            if (iter >= min_iters) {
+                total_event_count += vegas_batch_size;
+            }
+            ++i;
+        }
+        std::size_t job_id_before = _job_id;
+        std::size_t unweight_dispatched = start_jobs();
+        std::size_t in_flight = (_job_id - job_id_before) + unweight_dispatched;
+        register_dispatched_ids(job_id_before, _job_id);
+        done = true;
+        while (in_flight > 0) {
+            std::size_t job_id = _result_queue.wait();
+            _abort_check_function();
+            --in_flight;
+            auto& job = _running_jobs.at(job_id);
+            if (job.unweighted_events.size() == 0) {
+                // Generate-stage completion, globally ordered via
+                // _ready_gen/_commit_cursor.
+                _ready_gen.insert(job_id);
+                while (_ready_gen.erase(_commit_cursor) > 0) {
+                    auto& gen_job = _running_jobs.at(_commit_cursor);
+                    auto& channel = _channels.at(gen_job.channel_index);
+                    auto& channel_job_count =
+                        _channel_job_counts.at(gen_job.channel_index);
+                    if (channel_job_count == gen_job.split_job_count) {
+                        channel->clear_events();
+                    }
+                    channel->integrate(gen_job);
+                    update_integral();
+                    channel->update_max_weight(gen_job.weights);
+                    --channel_job_count;
+                    --_context_job_counts.at(gen_job.context_index);
+                    if (gen_job.unweight) {
+                        // Snapshot max_weight at this fixed commit position.
+                        channel->prepare_unweight_job(gen_job);
+                        _context_unweight_queue.at(gen_job.context_index)
+                            .push_back(gen_job.job_id);
+                    } else {
+                        done = false;
+                        if (channel_job_count == 0) {
+                            channel->optimize_vegas(gen_job);
+                            done_event_count += gen_job.vegas_batch_size;
+                        }
+                        _running_jobs.erase(_commit_cursor);
+                    }
+                    ++_commit_cursor;
+                }
+            } else {
+                // Unweight-stage completion, ordered per channel only.
+                auto& ready = _channel_unweight_ready.at(job.channel_index);
+                auto& cursor = _channel_unweight_cursor.at(job.channel_index);
+                ready.insert(job_id);
+                while (ready.erase(cursor) > 0) {
+                    auto& uw_job = _running_jobs.at(cursor);
+                    auto& channel = _channels.at(uw_job.channel_index);
+                    channel->write_events(uw_job.unweighted_events, uw_job.max_weight);
+                    update_counts();
+                    auto& channel_job_count =
+                        _channel_job_counts.at(uw_job.channel_index);
+                    --channel_job_count;
+                    --_context_job_counts.at(uw_job.context_index);
+                    if (channel_job_count == 0 &&
+                        channel->cross_section().rel_error() < target_precision) {
+                        done = false;
+                    }
+                    if (channel_job_count == 0) {
+                        channel->optimize_vegas(uw_job);
+                        done_event_count += uw_job.vegas_batch_size;
+                    }
+                    _running_jobs.erase(cursor);
+                    ++cursor;
+                }
+            }
+            std::size_t job_id_refill = _job_id;
+            std::size_t refill_unweight_dispatched = start_jobs();
+            in_flight += (_job_id - job_id_refill) + refill_unweight_dispatched;
+            register_dispatched_ids(job_id_refill, _job_id);
+            print_survey_update(false, done_event_count, total_event_count, iter);
+        }
+    }
+    print_survey_update(true, done_event_count, total_event_count, iter - 1);
+}
+
+// Sizes the next batch from committed data only, so the estimate never depends
+// on jobs still in flight.
+std::size_t EventGenerator::next_batch_job_count(std::size_t channel_index) const {
+    auto& status = _channels.at(channel_index)->status();
+    // No data yet: assume a pessimistic efficiency of 1.0.
+    double efficiency = status.count_opt > 0
+        ? status.count_unweighted / static_cast<double>(status.count_opt)
+        : 1.;
+    // Floor at 1 event so a near-zero observed efficiency can't blow this up.
+    double per_job_estimate =
+        std::max(efficiency * static_cast<double>(_config.cpu_batch_size), 1.);
+    double remaining = std::max(
+        static_cast<double>(status.count_target) - status.count_unweighted, 0.
+    );
+    double estimated = _config.generation_batch_fraction * remaining / per_job_estimate;
+    return std::max(
+        _config.min_batch_jobs, static_cast<std::size_t>(std::ceil(estimated))
+    );
+}
+
+std::size_t EventGenerator::start_jobs() {
     std::size_t ready_index = 0, context_index = 0;
+    std::size_t unweight_dispatched = 0;
     for (auto [context, job_count] : zip(_contexts, _context_job_counts)) {
         // fill the queue to twice the thread count to keep the worker threads busy
         std::size_t target_count = 2 * context->thread_pool().thread_count();
         std::size_t batch_size = context->device()->device_type() == DeviceType::cpu
             ? _config.cpu_batch_size
             : _config.gpu_batch_size;
+
+        // give priority to unweighting jobs to reduce memory usage
+        auto& unweight_queue = _context_unweight_queue.at(context_index);
+        std::size_t unweight_index = 0;
+        for (; job_count < target_count && unweight_index < unweight_queue.size();
+             ++unweight_index) {
+            auto& job = _running_jobs.at(unweight_queue.at(unweight_index));
+            _channels.at(job.channel_index)->submit_unweight_job(job, _result_queue);
+            ++job_count;
+            ++unweight_dispatched;
+        }
+        unweight_queue.erase(
+            unweight_queue.begin(), unweight_queue.begin() + unweight_index
+        );
+
         for (; job_count < target_count && ready_index < _ready_jobs.size();
              ++ready_index) {
             auto ready_job = _ready_jobs.at(ready_index);
@@ -247,22 +586,25 @@ bool EventGenerator::start_jobs() {
                 job.split_job_count = split_job_count * (1 + job.unweight);
                 job.job_id = _job_id;
                 job.context_index = context_index;
-                _channels.at(job.channel_index)->start_job(job, _result_queue);
+                _channels.at(job.channel_index)
+                    ->start_job(job, _result_queue, _seed, _survey_job, _survey_pass);
                 _channel_job_counts.at(job.channel_index) += 1 + job.unweight;
                 ++_job_id;
                 ++job_count;
             }
         }
-        if (ready_index == _ready_jobs.size()) {
-            break;
-        }
         ++context_index;
     }
     _ready_jobs.erase(_ready_jobs.begin(), _ready_jobs.begin() + ready_index);
-    return ready_index > 0;
+    return unweight_dispatched;
 }
 
 void EventGenerator::update_integral() {
+    update_integral_status();
+    update_integral_fractions();
+}
+
+void EventGenerator::update_integral_status() {
     double total_mean = 0., total_var = 0.;
     std::size_t total_count = 0, total_count_opt = 0;
     std::size_t total_count_after_cuts = 0, total_count_after_cuts_opt = 0;
@@ -293,19 +635,48 @@ void EventGenerator::update_integral() {
     _status.count_after_cuts_opt = total_count_after_cuts_opt;
     _status.iterations = iterations;
     _status.optimized = optimized;
+}
+
+void EventGenerator::update_integral_fractions() {
     for (auto [channel, integral_fraction] :
          zip(_channels, _channel_integral_fractions)) {
-        integral_fraction = channel->cross_section().mean() / total_mean;
-        channel->set_target_count(integral_fraction * _config.target_count);
+        integral_fraction = channel->cross_section().mean() / _status.mean;
+    }
+
+    // Distribute events between channels, ensure sum is exactly target_count
+    std::vector<std::size_t> counts(_channels.size());
+    std::vector<double> remainders(_channels.size());
+    std::size_t allocated = 0;
+    for (auto [count, remainder, fraction] :
+         zip(counts, remainders, _channel_integral_fractions)) {
+        double raw = fraction * _config.target_count;
+        double floored = std::floor(std::isfinite(raw) && raw > 0. ? raw : 0.);
+        count = static_cast<std::size_t>(floored);
+        remainder = raw - floored;
+        allocated += count;
+    }
+    std::vector<std::size_t> order(_channels.size());
+    std::iota(order.begin(), order.end(), 0);
+    std::sort(order.begin(), order.end(), [&](std::size_t a, std::size_t b) {
+        double rem_a = remainders.at(a), rem_b = remainders.at(b);
+        return rem_a == rem_b ? a < b : rem_a > rem_b;
+    });
+    std::size_t leftover =
+        allocated <= _config.target_count ? _config.target_count - allocated : 0;
+    leftover = std::min(leftover, order.size());
+    for (std::size_t index : order | std::views::take(leftover)) {
+        ++counts.at(index);
+    }
+    for (auto [channel, count] : zip(_channels, counts)) {
+        channel->set_target_count(count);
     }
 }
 
 void EventGenerator::update_counts() {
     double total_eff_count = 0.;
     bool done = true;
-    for (auto [channel, integral_fraction] :
-         zip(_channels, _channel_integral_fractions)) {
-        double chan_target = integral_fraction * _config.target_count;
+    for (auto& channel : _channels) {
+        std::size_t chan_target = channel->status().count_target;
         if (channel->status().count_unweighted < chan_target) {
             total_eff_count += channel->status().count_unweighted;
             done = false;
@@ -319,6 +690,7 @@ void EventGenerator::update_counts() {
 
 void EventGenerator::combine_to_compact_npy(const std::string& file_name) {
     reset_start_time();
+    MixMaxRandom select_rng(DerivedSeed(_seed, DerivedSeed::combine_select));
     auto [channel_data, particle_count, norm_factor] = init_combine();
     DataLayout layout(
         EventRecord::layout(
@@ -337,7 +709,7 @@ void EventGenerator::combine_to_compact_npy(const std::string& file_name) {
     print_combine_init();
     while (true) {
         _abort_check_function();
-        read_and_combine(channel_data, buffer, norm_factor);
+        read_and_combine(channel_data, buffer, norm_factor, select_rng);
         if (buffer.event_count() == 0) {
             break;
         }
@@ -355,8 +727,9 @@ void EventGenerator::combine_to_lhe_npy(
     const std::string& file_name, LHECompleter& lhe_completer
 ) {
     reset_start_time();
-    std::random_device rand_device;
-    std::mt19937 rand_gen(rand_device());
+    std::optional<std::uint64_t> seed = _seed;
+    MixMaxRandom select_rng(DerivedSeed(seed, DerivedSeed::combine_select));
+    MixMaxRandom rand_gen(DerivedSeed(seed, DerivedSeed::lhe_complete));
     auto [channel_data, particle_count, norm_factor] = init_combine();
     DataLayout in_layout(
         EventRecord::layout(
@@ -388,7 +761,7 @@ void EventGenerator::combine_to_lhe_npy(
     print_combine_init();
     while (true) {
         _abort_check_function();
-        read_and_combine(channel_data, buffer, norm_factor);
+        read_and_combine(channel_data, buffer, norm_factor, select_rng);
         if (buffer.event_count() == 0) {
             break;
         }
@@ -418,11 +791,9 @@ void EventGenerator::combine_to_lhe(
     const std::string& file_name, LHECompleter& lhe_completer, const LHEMeta& meta
 ) {
     reset_start_time();
+    std::optional<std::uint64_t> seed = _seed;
+    MixMaxRandom select_rng(DerivedSeed(seed, DerivedSeed::combine_select));
     ThreadPool pool(_config.combine_thread_count);
-    ThreadResource<std::mt19937> rand_gens(pool, []() {
-        std::random_device rand_device;
-        return std::mt19937(rand_device());
-    });
     auto [channel_data, particle_count, norm_factor] = init_combine();
     std::vector<std::pair<EventBuffer, std::string>> buffers;
     std::vector<std::size_t> idle_buffers;
@@ -444,45 +815,63 @@ void EventGenerator::combine_to_lhe(
     std::size_t event_count = 0;
     std::size_t last_update_count = 0;
     bool done = false;
+    // Batches are seeded by submission order and written in that same order
+    // (buffered in ready_slots until their turn), independent of completion timing.
+    std::size_t next_batch_seq = 0;
+    std::size_t write_cursor = 0;
+    std::vector<std::size_t> slot_batch_seq(buffers.size());
+    std::unordered_map<std::size_t, std::size_t> ready_slots;
     print_combine_init();
     while (true) {
         _abort_check_function();
         while (idle_buffers.size() > 0 && !done) {
-            std::size_t job_id = idle_buffers.back();
-            auto& [in_buffer, out_buffer] = buffers.at(job_id);
-            read_and_combine(channel_data, in_buffer, norm_factor);
+            std::size_t slot = idle_buffers.back();
+            auto& [in_buffer, out_buffer] = buffers.at(slot);
+            read_and_combine(channel_data, in_buffer, norm_factor, select_rng);
             if (in_buffer.event_count() == 0) {
                 done = true;
                 break;
             }
             idle_buffers.pop_back();
+            std::size_t batch_seq = next_batch_seq++;
+            slot_batch_seq.at(slot) = batch_seq;
             pool.submit(
-                [job_id, this, &in_buffer, &out_buffer, &lhe_completer, &rand_gens] {
+                [slot, batch_seq, seed, this, &in_buffer, &out_buffer, &lhe_completer] {
+                    MixMaxRandom rand_gen(
+                        DerivedSeed(seed, DerivedSeed::lhe_complete, batch_seq)
+                    );
                     LHEEvent lhe_event;
                     out_buffer.clear();
                     for (std::size_t i = 0; i < in_buffer.event_count(); ++i) {
                         fill_lhe_event(
-                            lhe_completer, lhe_event, in_buffer, i, rand_gens.get()
+                            lhe_completer, lhe_event, in_buffer, i, rand_gen
                         );
                         lhe_event.format_to(out_buffer);
                     }
-                    return job_id;
+                    return slot;
                 }
             );
         }
 
-        auto done_jobs = pool.wait_multiple();
-        for (std::size_t job_id : done_jobs) {
-            auto& [in_buffer, out_buffer] = buffers.at(job_id);
-            idle_buffers.push_back(job_id);
+        auto done_slots = pool.wait_multiple();
+        for (std::size_t slot : done_slots) {
+            ready_slots.emplace(slot_batch_seq.at(slot), slot);
+        }
+        for (auto it = ready_slots.find(write_cursor); it != ready_slots.end();
+             it = ready_slots.find(write_cursor)) {
+            std::size_t slot = it->second;
+            auto& [in_buffer, out_buffer] = buffers.at(slot);
             event_file.write_string(out_buffer);
             event_count += in_buffer.event_count();
+            idle_buffers.push_back(slot);
+            ready_slots.erase(it);
+            ++write_cursor;
             if (event_count - last_update_count > 10000) {
                 print_combine_update(event_count);
                 last_update_count = event_count;
             }
         }
-        if (done_jobs.size() == 0 && done) {
+        if (done_slots.size() == 0 && done) {
             break;
         }
     }
@@ -506,15 +895,17 @@ void EventGenerator::add_timing_data(const std::string& key) {
 }
 
 void EventGenerator::unweight_all() {
-    std::random_device rand_device;
-    std::mt19937 rand_gen(rand_device());
+    std::size_t call_index = _unweight_call_index++;
     bool done = true;
     double total_eff_count = 0.;
-    for (auto [channel, integral_fraction] :
-         zip(_channels, _channel_integral_fractions)) {
+    for (std::size_t channel_index = 0; auto& channel : _channels) {
+        MixMaxRandom rand_gen(
+            DerivedSeed(_seed, DerivedSeed::unweight_pass, call_index, channel_index)
+        );
         channel->unweight_file(rand_gen);
+        ++channel_index;
 
-        double chan_target = integral_fraction * _config.target_count;
+        std::size_t chan_target = channel->status().count_target;
         if (channel->status().count_unweighted < chan_target) {
             total_eff_count += channel->status().count_unweighted;
             done = false;
@@ -572,11 +963,12 @@ EventGenerator::init_combine() {
     std::size_t count_sum = 0;
     std::size_t particle_count = 0;
     double weight_sum = 0.;
-    for (auto [channel, integral_fraction] :
-         zip(_channels, _channel_integral_fractions)) {
+    for (auto& channel : _channels) {
         particle_count =
             std::max(particle_count, channel->event_file().particle_count());
-        std::size_t count = std::round(integral_fraction * _config.target_count);
+        // Exact apportioned target (see update_integral_fractions()), so counts
+        // sum to exactly _config.target_count.
+        std::size_t count = channel->status().count_target;
         count_sum += count;
         channel->event_file().seek(0);
         weight_sum += channel->channel_weight_sum(count);
@@ -596,7 +988,8 @@ EventGenerator::init_combine() {
 void EventGenerator::read_and_combine(
     std::vector<EventGenerator::CombineChannelData>& channel_data,
     EventBuffer& buffer,
-    double norm_factor
+    double norm_factor,
+    MixMaxRandom& rand_gen
 ) {
     std::size_t batch_size = 1000;
     std::size_t event_count = std::min(batch_size, channel_data.back().cum_count);
@@ -607,11 +1000,8 @@ void EventGenerator::read_and_combine(
     bool has_partial =
         _channels.at(0)->event_layout_extra_flags() & EventRecord::f_partial_weights;
 
-    std::random_device rand_device;
-    std::mt19937 rand_gen(rand_device());
     for (std::size_t event_index = 0; event_index < event_count; ++event_index) {
-        std::size_t random_index = std::uniform_int_distribution<
-            std::size_t>(0, channel_data.back().cum_count - 1)(rand_gen);
+        std::size_t random_index = rand_gen.generate_int(channel_data.back().cum_count);
         auto sampled_chan = std::lower_bound(
             channel_data.begin(),
             channel_data.end(),
@@ -688,7 +1078,7 @@ void EventGenerator::fill_lhe_event(
     LHEEvent& lhe_event,
     EventBuffer& buffer,
     std::size_t event_index,
-    std::mt19937& rand_gen
+    MixMaxRandom& rand_gen
 ) {
     EventRecord event_in = buffer.event(event_index);
     lhe_event.weight = event_in.weight();
@@ -739,14 +1129,14 @@ void EventGenerator::write_status(const std::string& status, bool force_write) {
     std::ofstream f(status_tmp_file);
     nlohmann::json j{
         {"status", status},
+        {"seed", _seed.value_or(0)},
         {"process", _status},
         {"channels", channel_status()},
         {"run_times", _timing_data},
         {"histograms", histograms()},
     };
     f << j.dump();
-    // rename atomically deletes the old file and replaces it with the new one
-    // such that the status file exists at all times
+    // Atomic rename keeps the status file present at all times.
     std::filesystem::rename(status_tmp_file, _status_file);
 }
 
diff --git a/madspace/src/driver/lhe_output.cpp b/madspace/src/driver/lhe_output.cpp
index 3f96379b7..f99ea8f73 100644
--- a/madspace/src/driver/lhe_output.cpp
+++ b/madspace/src/driver/lhe_output.cpp
@@ -394,7 +394,7 @@ void LHECompleter::complete_event_data(
     int color_index,
     int flavor_index,
     int helicity_index,
-    std::mt19937& rand_gen
+    MixMaxRandom& rand_gen
 ) {
     auto& subproc_data = _subproc_data.at(subprocess_index);
     if (event.particles.size() != subproc_data.particle_count) {
@@ -423,8 +423,7 @@ void LHECompleter::complete_event_data(
 
     auto [pdg_index, pdg_count] =
         _pdg_id_and_count.at(subproc_data.pdg_id_offset + flavor_index);
-    std::uniform_int_distribution<std::size_t> dist(0, pdg_count - 1);
-    std::size_t pdg_random = dist(rand_gen);
+    std::size_t pdg_random = rand_gen.generate_int(pdg_count);
     std::size_t pdg_offset = pdg_index + subproc_data.particle_count * pdg_random;
 
     for (std::size_t particle_index = 0; auto& particle : event.particles) {
diff --git a/madspace/src/driver/madnis_training.cpp b/madspace/src/driver/madnis_training.cpp
index c94d84fb9..665080be8 100644
--- a/madspace/src/driver/madnis_training.cpp
+++ b/madspace/src/driver/madnis_training.cpp
@@ -1,5 +1,6 @@
 #include "madspace/driver/madnis_training.hpp"
 
+#include "madspace/driver/random.hpp"
 #include "madspace/phasespace/batch_sampler.hpp"
 
 using namespace madspace;
@@ -9,13 +10,17 @@ MadnisTraining::MadnisTraining(
     ContextPtr optimizer_context,
     const Config& config,
     const std::vector<std::shared_ptr<Integrand>>& integrands,
-    const std::optional<ChannelWeightNetwork>& cwnet
+    const std::optional<ChannelWeightNetwork>& cwnet,
+    std::optional<std::uint64_t> seed,
+    std::size_t channel_index_offset
 ) :
     _generator_context(generator_context),
     _optimizer_context(optimizer_context),
     _config(config),
     _channels(integrands.size()),
-    _cwnet(cwnet) {
+    _cwnet(cwnet),
+    _seed(seed),
+    _channel_index_offset(channel_index_offset) {
     for (std::size_t index = 0;
          auto [integrand, channel] : zip(integrands, _channels)) {
         channel.index = index;
@@ -54,7 +59,7 @@ void MadnisTraining::train_step(std::size_t batch_index) {
         _config.buffer_capacity > 0 && batch_index % (_config.buffered_steps + 1) != 0;
     TensorVec training_batch;
     while (true) {
-        start_generator_jobs(channel_sizes);
+        maybe_start_generator_jobs(channel_sizes, !try_buffered);
         if (try_buffered) {
             if (check_buffered_training_batch(channel_sizes)) {
                 training_batch = build_buffered_training_batch(channel_sizes);
@@ -71,22 +76,26 @@ void MadnisTraining::train_step(std::size_t batch_index) {
     update_history(results, channel_sizes);
     if (_channels.size() > 0 && _cwnet &&
         (batch_index + 1) % _config.channel_dropping_interval == 0) {
-        std::vector<std::size_t> job_ids;
-        while ((job_ids = gen_thread_pool.wait_multiple()).size() != 0) {
-            process_job_results(job_ids);
-        }
+        process_all_jobs();
         drop_channels();
     }
     if (batch_index ==
         static_cast<std::size_t>(
             (1 - _config.fixed_cwnet_fraction) * _config.batches
         )) {
-        std::vector<std::size_t> job_ids;
-        while ((job_ids = gen_thread_pool.wait_multiple()).size() != 0) {
-            process_job_results(job_ids);
-        }
+        process_all_jobs();
         freeze_cwnet();
     }
+    if (batch_index + 1 == _config.batches) {
+        process_all_jobs();
+    }
+}
+
+void MadnisTraining::process_all_jobs() {
+    std::vector<std::size_t> job_ids;
+    while ((job_ids = _generator_context->thread_pool().wait_multiple()).size() != 0) {
+        process_job_results(job_ids);
+    }
 }
 
 std::vector<std::size_t> MadnisTraining::active_channels() const {
@@ -262,9 +271,18 @@ void MadnisTraining::start_generator_jobs(const std::vector<std::size_t>& counts
     if (_running_jobs.size() > 0) {
         return;
     }
+    bool is_gpu = _generator_context->device()->device_type() != DeviceType::cpu;
+    if (_config.reproducible) {
+        // flush buffer samples staged since the last round (see process_job_results)
+        for (auto& channel : _channels) {
+            for (auto& pending : channel.pending_buffer_samples) {
+                buffer_store(channel, pending);
+            }
+            channel.pending_buffer_samples.clear();
+        }
+    }
     _generator_params.copy_from(_optimizer->parameters());
     std::size_t chan_count = counts.size();
-    bool is_gpu = _generator_context->device()->device_type() != DeviceType::cpu;
     std::size_t batch_size = is_gpu
         ? _config.gpu_generator_batch_granularity
         : _config.cpu_generator_batch_size;
@@ -280,14 +298,29 @@ void MadnisTraining::start_generator_jobs(const std::vector<std::size_t>& counts
             ? (target_count - channel.sample_count + batch_size - 1) / batch_size
             : 0;
     }
-    std::size_t available_jobs = _generator_context->thread_pool().thread_count();
-    std::vector<std::size_t> channel_sizes;
     std::size_t gpu_subbatches =
         (_config.gpu_generator_batch_size + _config.gpu_generator_batch_granularity -
          1) /
         _config.gpu_generator_batch_granularity;
+    // reproducible mode: dispatch exactly what this round needs, uncapped by
+    // thread count, so round contents don't depend on thread count
+    std::size_t available_jobs;
+    if (_config.reproducible) {
+        available_jobs = 0;
+        for (auto count : missing_batch_counts) {
+            available_jobs += count;
+        }
+        for (auto count : target_batch_counts) {
+            available_jobs += count;
+        }
+    } else {
+        available_jobs = _generator_context->thread_pool().thread_count();
+        if (is_gpu) {
+            available_jobs *= gpu_subbatches;
+        }
+    }
+    std::vector<std::size_t> channel_sizes;
     if (is_gpu) {
-        available_jobs *= gpu_subbatches;
         channel_sizes.resize(chan_count, 0);
     }
 
@@ -349,6 +382,32 @@ void MadnisTraining::start_generator_jobs(const std::vector<std::size_t>& counts
     }
 }
 
+// Reproducible mode (CPU only) additionally requires an online attempt and the
+// online cache to be exhausted or about to be, both deterministic conditions.
+void MadnisTraining::maybe_start_generator_jobs(
+    const std::vector<std::size_t>& counts, bool is_online_attempt
+) {
+    if (_running_jobs.size() > 0) {
+        return;
+    }
+    if (_config.reproducible) {
+        if (!is_online_attempt) {
+            return;
+        }
+        bool depletion_imminent = false;
+        for (auto [channel, count] : zip(_channels, counts)) {
+            if (count >= channel.sample_count) {
+                depletion_imminent = true;
+                break;
+            }
+        }
+        if (!depletion_imminent) {
+            return;
+        }
+    }
+    start_generator_jobs(counts);
+}
+
 TensorVec MadnisTraining::permute_tensors(const TensorVec& tensors) const {
     TensorVec ret;
     ret.reserve(tensors.size());
@@ -358,20 +417,42 @@ TensorVec MadnisTraining::permute_tensors(const TensorVec& tensors) const {
     return ret;
 }
 
+// Seed depends only on channel index and channel-local sequence, not on which
+// worker thread runs the job.
 void MadnisTraining::start_single_job(
     std::size_t channel_index, std::size_t batch_size
 ) {
     std::size_t job_id = _job_id;
     ++_job_id;
+    auto& channel = _channels.at(channel_index);
+    std::size_t channel_seq = channel.next_dispatch_seq++;
+    std::size_t global_channel_index = _channel_index_offset + channel_index;
     auto& job = std::get<0>(_running_jobs.emplace(job_id, SampleJob{}))->second;
+    job.dispatch_seq = channel_seq;
     _generator_context->thread_pool().submit(
-        [this, channel_index, batch_size, job_id, &job]() {
+        [this, channel_index, global_channel_index, batch_size, job_id, &job]() {
             auto& channel = _channels.at(channel_index);
+            if (_seed) {
+                channel.generator_runtime->set_seed(DerivedSeed(
+                    _seed,
+                    DerivedSeed::madnis_generate,
+                    job.dispatch_seq,
+                    global_channel_index
+                ));
+            }
             auto samples = channel.generator_runtime->run({Tensor({batch_size})});
             job.samples.tensors = permute_tensors(samples);
             job.samples.size = samples.at(0).size(0);
             job.samples.channel_index = channel_index;
             if (channel.unweighter_runtime) {
+                if (_seed) {
+                    channel.unweighter_runtime->set_seed(DerivedSeed(
+                        _seed,
+                        DerivedSeed::madnis_unweight,
+                        job.dispatch_seq,
+                        global_channel_index
+                    ));
+                }
                 auto unw_samples = channel.unweighter_runtime->run(samples);
                 job.unweighted_samples.tensors = permute_tensors(unw_samples);
                 job.unweighted_samples.size = unw_samples.at(0).size(0);
@@ -382,15 +463,34 @@ void MadnisTraining::start_single_job(
     );
 }
 
+// Dispatch sequence is global (not per-channel), since one job spans all channels.
 void MadnisTraining::start_multi_job(const std::vector<std::size_t> batch_sizes) {
     std::size_t job_id = _job_id;
     ++_job_id;
+    std::size_t dispatch_seq = _multi_job_next_dispatch_seq++;
     auto& job = std::get<0>(_running_jobs.emplace(job_id, SampleJob{}))->second;
+    job.dispatch_seq = dispatch_seq;
     _generator_context->thread_pool().submit([this, batch_sizes, job_id, &job]() {
+        if (_seed) {
+            _multi_channel_generator->set_seed(DerivedSeed(
+                _seed,
+                DerivedSeed::madnis_generate,
+                job.dispatch_seq,
+                _channel_index_offset
+            ));
+        }
         auto samples = _multi_channel_generator->run({Tensor(batch_sizes)});
         job.samples.tensors = permute_tensors(samples);
         job.samples.channel_sizes = samples.back().batch_sizes();
         if (_multi_channel_unweighter) {
+            if (_seed) {
+                _multi_channel_unweighter->set_seed(DerivedSeed(
+                    _seed,
+                    DerivedSeed::madnis_unweight,
+                    job.dispatch_seq,
+                    _channel_index_offset
+                ));
+            }
             auto unw_samples = _multi_channel_unweighter->run(samples);
             job.unweighted_samples.tensors = permute_tensors(unw_samples);
             job.unweighted_samples.channel_sizes = unw_samples.back().batch_sizes();
@@ -502,54 +602,94 @@ MadnisTraining::build_buffered_training_batch(const std::vector<size_t>& counts)
         }
     }
     args.emplace_back(counts);
+    std::size_t batch_seq = _buffered_batch_seq++;
+    if (_seed) {
+        _multi_channel_sampler->set_seed(DerivedSeed(
+            _seed, DerivedSeed::madnis_sample_buffer, batch_seq, _channel_index_offset
+        ));
+    }
     return _multi_channel_sampler->run(args);
 }
 
 void MadnisTraining::process_job_results(const std::vector<std::size_t>& job_ids) {
     for (auto job_id : job_ids) {
-        auto job = std::move(_running_jobs.extract(job_id).mapped());
+        auto& job = _running_jobs.at(job_id);
+        // mark ready, committed below strictly in dispatch order
         if (job.samples.channel_sizes.size() == 0) {
-            auto& channel = _channels.at(job.samples.channel_index);
-            channel.sample_count += job.samples.size;
-            channel.sample_batches.push_back(std::move(job.samples));
-            if (job.unweighted_samples.size > 0) {
-                buffer_store(channel, job.unweighted_samples);
-            }
+            _channels.at(job.samples.channel_index)
+                .ready_job_ids.emplace(job.dispatch_seq, job_id);
         } else {
-            std::size_t offset = 0, unw_offset = 0, chan_index = 0;
-            SampleBatch chan_unweighted_samples;
-            for (auto [channel, chan_size] :
-                 zip(_channels, job.samples.channel_sizes)) {
-                if (chan_size == 0) {
-                    ++chan_index;
-                    continue;
+            _multi_job_ready_job_ids.emplace(job.dispatch_seq, job_id);
+        }
+    }
+
+    // commit each channel's ready single-channel jobs strictly in dispatch order
+    for (auto& channel : _channels) {
+        for (auto it = channel.ready_job_ids.find(channel.commit_cursor);
+             it != channel.ready_job_ids.end();
+             it = channel.ready_job_ids.find(channel.commit_cursor)) {
+            auto committed_job = std::move(_running_jobs.extract(it->second).mapped());
+            channel.ready_job_ids.erase(it);
+            ++channel.commit_cursor;
+            channel.sample_count += committed_job.samples.size;
+            channel.sample_batches.push_back(std::move(committed_job.samples));
+            if (committed_job.unweighted_samples.size > 0) {
+                if (_config.reproducible) {
+                    // flushed into buffer at the start of the next round instead
+                    channel.pending_buffer_samples.push_back(
+                        std::move(committed_job.unweighted_samples)
+                    );
+                } else {
+                    buffer_store(channel, committed_job.unweighted_samples);
                 }
-                channel.sample_count += chan_size;
-                channel.sample_batches.emplace_back();
-                auto& batch = channel.sample_batches.back();
-                batch.tensors.reserve(job.samples.tensors.size());
-                for (auto& tensor : job.samples.tensors) {
-                    batch.tensors.push_back(
-                        tensor.slice(0, offset, offset + chan_size)
+            }
+        }
+    }
+
+    // commit ready multi-channel (GPU) jobs strictly in dispatch order
+    for (auto it = _multi_job_ready_job_ids.find(_multi_job_commit_cursor);
+         it != _multi_job_ready_job_ids.end();
+         it = _multi_job_ready_job_ids.find(_multi_job_commit_cursor)) {
+        auto multi_job = std::move(_running_jobs.extract(it->second).mapped());
+        _multi_job_ready_job_ids.erase(it);
+        ++_multi_job_commit_cursor;
+        std::size_t offset = 0, unw_offset = 0, chan_index = 0;
+        SampleBatch chan_unweighted_samples;
+        for (auto [channel, chan_size] :
+             zip(_channels, multi_job.samples.channel_sizes)) {
+            if (chan_size == 0) {
+                ++chan_index;
+                continue;
+            }
+            channel.sample_count += chan_size;
+            channel.sample_batches.emplace_back();
+            auto& batch = channel.sample_batches.back();
+            batch.tensors.reserve(multi_job.samples.tensors.size());
+            for (auto& tensor : multi_job.samples.tensors) {
+                batch.tensors.push_back(tensor.slice(0, offset, offset + chan_size));
+            }
+            if (multi_job.unweighted_samples.channel_sizes.size() > 0) {
+                std::size_t unw_chan_size =
+                    multi_job.unweighted_samples.channel_sizes.at(chan_index);
+                chan_unweighted_samples.tensors.clear();
+                chan_unweighted_samples.size = unw_chan_size;
+                for (auto& tensor : multi_job.unweighted_samples.tensors) {
+                    chan_unweighted_samples.tensors.push_back(
+                        tensor.slice(0, unw_offset, unw_offset + unw_chan_size)
                     );
                 }
-                if (job.unweighted_samples.channel_sizes.size() > 0) {
-                    std::size_t unw_chan_size =
-                        job.unweighted_samples.channel_sizes.at(chan_index);
-                    chan_unweighted_samples.tensors.clear();
-                    chan_unweighted_samples.size = unw_chan_size;
-                    for (auto& tensor : job.unweighted_samples.tensors) {
-                        chan_unweighted_samples.tensors.push_back(
-                            tensor.slice(0, unw_offset, unw_offset + unw_chan_size)
-                        );
-                    }
+                if (_config.reproducible) {
+                    channel.pending_buffer_samples.push_back(
+                        std::move(chan_unweighted_samples)
+                    );
+                } else {
                     buffer_store(channel, chan_unweighted_samples);
-                    unw_offset += unw_chan_size;
                 }
-                batch.size = chan_size;
-                offset += chan_size;
-                ++chan_index;
+                unw_offset += unw_chan_size;
             }
+            batch.size = chan_size;
+            offset += chan_size;
+            ++chan_index;
         }
     }
 }
@@ -642,7 +782,8 @@ void MadnisTraining::drop_channels() {
     }
     std::vector<std::size_t> indices(_channels.size());
     std::iota(indices.begin(), indices.end(), 0);
-    std::sort(indices.begin(), indices.end(), [&](auto i, auto j) {
+    // stable_sort: exact ties (e.g. charge-conjugate channels) must break consistently
+    std::stable_sort(indices.begin(), indices.end(), [&](auto i, auto j) {
         return abs_means.at(i) < abs_means.at(j);
     });
 
@@ -668,8 +809,7 @@ void MadnisTraining::drop_channels() {
                         return _active_flavors_count.at(flav_index) == 0;
                     }
                 )) {
-                // cannot drop this channel because one of its flavors is not
-                // available in any other channel
+                // a flavor of this channel has no other channel to fall back to
                 continue;
             }
             for (std::size_t flav_index : active_flavors) {
@@ -719,14 +859,24 @@ MultiMadnisTraining::MultiMadnisTraining(
     ContextPtr optimizer_context,
     const MadnisTraining::Config& config,
     const nested_vector2<std::shared_ptr<Integrand>>& integrands,
-    const std::vector<std::optional<ChannelWeightNetwork>>& cwnets
+    const std::vector<std::optional<ChannelWeightNetwork>>& cwnets,
+    std::optional<std::uint64_t> seed
 ) :
     _config(config) {
     _subprocesses.reserve(integrands.size());
+    // each subprocess gets its own channel_index slice so their streams don't collide
+    std::size_t channel_index_offset = 0;
     for (auto [integ, cwnet] : zip(integrands, cwnets)) {
         _subprocesses.emplace_back(
-            generator_context, optimizer_context, config, integ, cwnet
+            generator_context,
+            optimizer_context,
+            config,
+            integ,
+            cwnet,
+            seed,
+            channel_index_offset
         );
+        channel_index_offset += integ.size();
     }
 }
 
diff --git a/madspace/src/driver/random.cpp b/madspace/src/driver/random.cpp
new file mode 100644
index 000000000..d27c60e8c
--- /dev/null
+++ b/madspace/src/driver/random.cpp
@@ -0,0 +1,50 @@
+#include "madspace/driver/random.hpp"
+
+#include <format>
+#include <random>
+
+namespace madspace {
+
+DerivedSeed::DerivedSeed(
+    const std::optional<std::uint64_t>& seed,
+    SeedType seed_type,
+    std::size_t job_index,
+    std::size_t channel_index,
+    std::size_t stream_index
+) {
+    std::uint64_t main_seed;
+    if (seed) {
+        main_seed = seed.value();
+    } else {
+        std::random_device rand_device;
+        main_seed = rand_device();
+    }
+    if (channel_index >= max_channel_count) {
+        throw std::invalid_argument(
+            std::format(
+                "channel index {} exceeded {}", channel_index, max_channel_count - 1
+            )
+        );
+    }
+    if (job_index >= max_job_count) {
+        throw std::invalid_argument(
+            std::format("job index {} exceeded {}", job_index, max_job_count - 1)
+        );
+    }
+    if (stream_index >= max_stream_count) {
+        throw std::invalid_argument(
+            std::format(
+                "stream index {} exceeded {}", stream_index, max_stream_count - 1
+            )
+        );
+    }
+    seed_parts[0] = static_cast<std::uint32_t>(main_seed >> 32);
+    seed_parts[1] = static_cast<std::uint32_t>(main_seed & 0xFFFFFFFFULL);
+    seed_parts[2] = job_index;
+    seed_parts[3] =
+        ((static_cast<std::uint32_t>(seed_type) << 28) |
+         (static_cast<std::uint32_t>(channel_index) << 16) |
+         static_cast<std::uint32_t>(stream_index));
+}
+
+} // namespace madspace
diff --git a/madspace/src/gpu/runtime.cu b/madspace/src/gpu/runtime.cu
index 87ff0b0be..0e01bd1e2 100644
--- a/madspace/src/gpu/runtime.cu
+++ b/madspace/src/gpu/runtime.cu
@@ -1610,11 +1610,14 @@ GpuRuntime::GpuRuntime(const Function& function_arg, ContextPtr context) :
     );
 }
 
-TensorVec GpuRuntime::run(const TensorVec& inputs) {
+TensorVec GpuRuntime::run(const TensorVec& inputs, std::optional<std::uint64_t> seed) {
     auto& gpu_device = *static_cast<const GpuDevice*>(_context->device());
     auto& streams = _streams.get();
     auto& events = _events.get();
     gpu_device.activate();
+    if (seed) {
+        check_error(gpurandSetPseudoRandomGeneratorSeed(gpurand_generator(), *seed));
+    }
     auto locals = _locals_init;
     std::copy(inputs.begin(), inputs.end(), locals.begin());
     gpuStream_t main_stream = streams.at(0);
@@ -1650,12 +1653,17 @@ TensorVec GpuRuntime::run(const TensorVec& inputs) {
 }
 
 std::tuple<TensorVec, TensorVec, std::vector<bool>> GpuRuntime::run_with_grad(
-    const TensorVec& inputs, const std::vector<bool>& input_requires_grad
+    const TensorVec& inputs,
+    const std::vector<bool>& input_requires_grad,
+    std::optional<std::uint64_t> seed
 ) {
     auto& gpu_device = *static_cast<const GpuDevice*>(_context->device());
     auto& streams = _streams.get();
     auto& events = _events.get();
     gpu_device.activate();
+    if (seed) {
+        check_error(gpurandSetPseudoRandomGeneratorSeed(gpurand_generator(), *seed));
+    }
     auto locals = _locals_init;
     auto requires_grad = _requires_grad_init;
     std::vector<bool> store_local(locals.size());
diff --git a/madspace/src/gpu/runtime.hpp b/madspace/src/gpu/runtime.hpp
index 30930de26..567b28b1b 100644
--- a/madspace/src/gpu/runtime.hpp
+++ b/madspace/src/gpu/runtime.hpp
@@ -6,6 +6,7 @@
 #include "madspace/driver/tensor.hpp"
 
 #include <memory>
+#include <optional>
 
 namespace madspace {
 namespace gpu {
@@ -31,9 +32,13 @@ class GpuRuntime : public Runtime {
     };
 
     GpuRuntime(const Function& function, ContextPtr context);
-    TensorVec run(const TensorVec& inputs) override;
+    TensorVec
+    run(const TensorVec& inputs,
+        std::optional<std::uint64_t> seed = std::nullopt) override;
     std::tuple<TensorVec, TensorVec, std::vector<bool>> run_with_grad(
-        const TensorVec& inputs, const std::vector<bool>& input_requires_grad
+        const TensorVec& inputs,
+        const std::vector<bool>& input_requires_grad,
+        std::optional<std::uint64_t> seed = std::nullopt
     ) override;
     std::pair<TensorVec, TensorVec> run_backward(
         const TensorVec& output_grads,
diff --git a/madspace/src/phasespace/channel_weight_network.cpp b/madspace/src/phasespace/channel_weight_network.cpp
index d84a7a389..4a71ad747 100644
--- a/madspace/src/phasespace/channel_weight_network.cpp
+++ b/madspace/src/phasespace/channel_weight_network.cpp
@@ -67,8 +67,10 @@ NamedVector<Value> ChannelWeightNetwork::build_function_impl(
     };
 }
 
-void ChannelWeightNetwork::initialize_globals(ContextPtr context) const {
-    _mlp.initialize_globals(context);
+void ChannelWeightNetwork::initialize_globals(
+    ContextPtr context, std::optional<std::uint64_t> seed
+) const {
+    _mlp.initialize_globals(context, seed);
 
     context->define_global(_mask_name, DataType::dt_float, {_channel_count});
     bool is_cpu = context->device() == cpu_device();
diff --git a/madspace/src/phasespace/discrete_flow.cpp b/madspace/src/phasespace/discrete_flow.cpp
index e7119e2e1..5f283fcb1 100644
--- a/madspace/src/phasespace/discrete_flow.cpp
+++ b/madspace/src/phasespace/discrete_flow.cpp
@@ -139,13 +139,15 @@ Mapping::Result DiscreteFlow::build_transform(
     };
 }
 
-void DiscreteFlow::initialize_globals(ContextPtr context) const {
+void DiscreteFlow::initialize_globals(
+    ContextPtr context, std::optional<std::uint64_t> seed
+) const {
     if (_first_prob_name) {
         initialize_uniform_probs(
             context, _first_prob_name.value(), _option_counts.at(0)
         );
     }
     for (auto& subnet : _subnets) {
-        subnet.initialize_globals(context);
+        subnet.initialize_globals(context, seed);
     }
 }
diff --git a/madspace/src/phasespace/flow.cpp b/madspace/src/phasespace/flow.cpp
index d3eefd906..1f353fcf1 100644
--- a/madspace/src/phasespace/flow.cpp
+++ b/madspace/src/phasespace/flow.cpp
@@ -211,17 +211,19 @@ Flow::Flow(
     }
 }
 
-void Flow::initialize_globals(ContextPtr context) const {
+void Flow::initialize_globals(
+    ContextPtr context, std::optional<std::uint64_t> seed
+) const {
     for (auto& block : _coupling_blocks) {
-        block.subnet1.initialize_globals(context);
-        block.subnet2.initialize_globals(context);
+        block.subnet1.initialize_globals(context, seed);
+        block.subnet2.initialize_globals(context, seed);
     }
 }
 
 void Flow::initialize_from_vegas(
-    ContextPtr context, const std::string& grid_name
+    ContextPtr context, const std::string& grid_name, std::optional<std::uint64_t> seed
 ) const {
-    initialize_globals(context);
+    initialize_globals(context, seed);
     auto& last_block = _coupling_blocks.at(_coupling_blocks.size() - 1);
     vegas_init(
         context,
diff --git a/madspace/src/phasespace/mlp.cpp b/madspace/src/phasespace/mlp.cpp
index 920f36594..4de6a6537 100644
--- a/madspace/src/phasespace/mlp.cpp
+++ b/madspace/src/phasespace/mlp.cpp
@@ -1,7 +1,8 @@
 #include "madspace/phasespace/mlp.hpp"
 
 #include <format>
-#include <random>
+
+#include "madspace/driver/random.hpp"
 
 using namespace madspace;
 
@@ -51,11 +52,9 @@ void initialize_layer(
     std::size_t output_dim,
     const std::string& prefix,
     int layer_index,
-    std::mt19937& rand_gen,
+    std::optional<std::uint64_t> seed,
     bool zeros
 ) {
-    double bound = 1 / std::sqrt(input_dim);
-    std::uniform_real_distribution<double> rand_dist(-bound, bound);
     auto weight_name =
         prefixed_name(prefix, std::format("layer{}.weight", layer_index));
     auto bias_name = prefixed_name(prefix, std::format("layer{}.bias", layer_index));
@@ -75,14 +74,35 @@ void initialize_layer(
     }
 
     auto weight_view = weight_tensor.view<double, 3>()[0];
-    for (std::size_t i = 0; i < output_dim; ++i) {
-        for (std::size_t j = 0; j < input_dim; ++j) {
-            weight_view[i][j] = zeros ? 0. : rand_dist(rand_gen);
-        }
-    }
     auto bias_view = bias_tensor.view<double, 2>()[0];
-    for (std::size_t i = 0; i < output_dim; ++i) {
-        bias_view[i] = zeros ? 0. : rand_dist(rand_gen);
+    if (zeros) {
+        for (std::size_t i = 0; i < output_dim; ++i) {
+            for (std::size_t j = 0; j < input_dim; ++j) {
+                weight_view[i][j] = 0.;
+            }
+            bias_view[i] = 0.;
+        }
+    } else {
+        double bound = 1 / std::sqrt(input_dim);
+        auto uniform = [&](MixMaxRandom& gen) {
+            return bound * (2. * gen.generate_double() - 1.);
+        };
+        // Independently seeded per tensor, via a fresh unique_seed_index() rather
+        // than a hash of the tensor's name.
+        MixMaxRandom weight_rand_gen(
+            DerivedSeed(seed, DerivedSeed::global_init, context->unique_seed_index())
+        );
+        for (std::size_t i = 0; i < output_dim; ++i) {
+            for (std::size_t j = 0; j < input_dim; ++j) {
+                weight_view[i][j] = uniform(weight_rand_gen);
+            }
+        }
+        MixMaxRandom bias_rand_gen(
+            DerivedSeed(seed, DerivedSeed::global_init, context->unique_seed_index())
+        );
+        for (std::size_t i = 0; i < output_dim; ++i) {
+            bias_view[i] = uniform(bias_rand_gen);
+        }
     }
 
     if (!is_cpu) {
@@ -133,15 +153,15 @@ MLP::build_function_impl(FunctionBuilder& fb, const NamedVector<Value>& args) co
     };
 }
 
-void MLP::initialize_globals(ContextPtr context) const {
-    std::random_device rand_device;
-    std::mt19937 rand_gen(rand_device());
+void MLP::initialize_globals(
+    ContextPtr context, std::optional<std::uint64_t> seed
+) const {
     std::size_t dim = _input_dim;
     for (std::size_t i = 1; i < _layers; ++i) {
-        initialize_layer(context, dim, _hidden_dim, _prefix, i, rand_gen, false);
+        initialize_layer(context, dim, _hidden_dim, _prefix, i, seed, false);
         dim = _hidden_dim;
     }
-    initialize_layer(context, dim, _output_dim, _prefix, _layers, rand_gen, true);
+    initialize_layer(context, dim, _output_dim, _prefix, _layers, seed, true);
 }
 
 std::vector<std::string> MLP::global_names() const {
diff --git a/madspace/src/python/madspace.cpp b/madspace/src/python/madspace.cpp
index 72c029e7c..0beb10945 100644
--- a/madspace/src/python/madspace.cpp
+++ b/madspace/src/python/madspace.cpp
@@ -842,7 +842,12 @@ PYBIND11_MODULE(_madspace_py, m) {
     )
         .def("input_dim", &MLP::input_dim)
         .def("output_dim", &MLP::output_dim)
-        .def("initialize_globals", &MLP::initialize_globals, py::arg("context"));
+        .def(
+            "initialize_globals",
+            &MLP::initialize_globals,
+            py::arg("context"),
+            py::arg("seed") = std::nullopt
+        );
 
     py::classh<Flow, Mapping>(m, "Flow")
         .def(
@@ -866,12 +871,18 @@ PYBIND11_MODULE(_madspace_py, m) {
         )
         .def("input_dim", &Flow::input_dim)
         .def("condition_dim", &Flow::condition_dim)
-        .def("initialize_globals", &Flow::initialize_globals, py::arg("context"))
+        .def(
+            "initialize_globals",
+            &Flow::initialize_globals,
+            py::arg("context"),
+            py::arg("seed") = std::nullopt
+        )
         .def(
             "initialize_from_vegas",
             &Flow::initialize_from_vegas,
             py::arg("context"),
-            py::arg("grid_name")
+            py::arg("grid_name"),
+            py::arg("seed") = std::nullopt
         );
 
     py::classh<PropagatorChannelWeights, FunctionGenerator>(
@@ -927,7 +938,8 @@ PYBIND11_MODULE(_madspace_py, m) {
         .def(
             "initialize_globals",
             &ChannelWeightNetwork::initialize_globals,
-            py::arg("context")
+            py::arg("context"),
+            py::arg("seed") = std::nullopt
         );
 
     py::classh<DiscreteHistogram, FunctionGenerator>(m, "DiscreteHistogram")
@@ -972,7 +984,10 @@ PYBIND11_MODULE(_madspace_py, m) {
         .def("option_counts", &DiscreteFlow::option_counts)
         .def("condition_dim", &DiscreteFlow::condition_dim)
         .def(
-            "initialize_globals", &DiscreteFlow::initialize_globals, py::arg("context")
+            "initialize_globals",
+            &DiscreteFlow::initialize_globals,
+            py::arg("context"),
+            py::arg("seed") = std::nullopt
         );
 
     py::classh<VegasGridOptimizer>(m, "VegasGridOptimizer")
@@ -1363,7 +1378,8 @@ PYBIND11_MODULE(_madspace_py, m) {
         )
         .def_readwrite(
             "softclip_threshold", &MadnisTraining::Config::softclip_threshold
-        );
+        )
+        .def_readwrite("reproducible", &MadnisTraining::Config::reproducible);
 
     py::classh<MadnisTraining>(m, "MadnisTraining")
         .def(
@@ -1372,12 +1388,14 @@ PYBIND11_MODULE(_madspace_py, m) {
                 ContextPtr,
                 const MadnisTraining::Config&,
                 const std::vector<std::shared_ptr<Integrand>>&,
-                const std::optional<ChannelWeightNetwork>&>(),
+                const std::optional<ChannelWeightNetwork>&,
+                std::optional<std::uint64_t>>(),
             py::arg("generator_context"),
             py::arg("optimizer_context"),
             py::arg("config"),
             py::arg("integrands"),
-            py::arg("cwnet")
+            py::arg("cwnet"),
+            py::arg("seed") = std::nullopt
         )
         .def("train_step", &MadnisTraining::train_step, py::arg("batch_index"))
         .def("active_channels", &MadnisTraining::active_channels)
@@ -1390,12 +1408,14 @@ PYBIND11_MODULE(_madspace_py, m) {
                 ContextPtr,
                 const MadnisTraining::Config&,
                 const nested_vector2<std::shared_ptr<Integrand>>&,
-                const std::vector<std::optional<ChannelWeightNetwork>>&>(),
+                const std::vector<std::optional<ChannelWeightNetwork>>&,
+                std::optional<std::uint64_t>>(),
             py::arg("generator_context"),
             py::arg("optimizer_context"),
             py::arg("config"),
             py::arg("integrands"),
-            py::arg("cwnets")
+            py::arg("cwnets"),
+            py::arg("seed") = std::nullopt
         )
         .def("train", &MultiMadnisTraining::train)
         .def("active_channels", &MultiMadnisTraining::active_channels);
@@ -1429,7 +1449,11 @@ PYBIND11_MODULE(_madspace_py, m) {
         .def_readwrite(
             "cut_efficiency_threshold", &GeneratorConfig::cut_efficiency_threshold
         )
-        .def_readwrite("max_cut_repetitions", &GeneratorConfig::max_cut_repetitions);
+        .def_readwrite("max_cut_repetitions", &GeneratorConfig::max_cut_repetitions)
+        .def_readwrite(
+            "generation_batch_fraction", &GeneratorConfig::generation_batch_fraction
+        )
+        .def_readwrite("min_batch_jobs", &GeneratorConfig::min_batch_jobs);
 
     py::classh<GeneratorStatus>(m, "GeneratorStatus")
         .def(py::init<>())
@@ -1611,9 +1635,9 @@ PYBIND11_MODULE(_madspace_py, m) {
         .def_readwrite("pdg_color_types", &LHECompleter::SubprocArgs::pdg_color_types)
         .def_readwrite("helicities", &LHECompleter::SubprocArgs::helicities)
         .def_readwrite("pdg_ids", &LHECompleter::SubprocArgs::pdg_ids);
-    py::classh<std::mt19937>(m, "RandGen")
+    py::classh<MixMaxRandom>(m, "MixMaxRandom")
         .def(py::init<>())
-        .def(py::init<std::mt19937::result_type>(), py::arg("seed"));
+        .def(py::init<std::uint64_t>(), py::arg("seed"));
     py::classh<LHECompleter>(m, "LHECompleter")
         .def(
             py::init<const std::vector<LHECompleter::SubprocArgs>&, double>(),
@@ -1712,7 +1736,8 @@ PYBIND11_MODULE(_madspace_py, m) {
                 const std::vector<ContextPtr>&,
                 const std::vector<std::shared_ptr<ChannelEventGenerator>>&,
                 const std::string&,
-                const GeneratorConfig&>(),
+                const GeneratorConfig&,
+                std::optional<std::uint64_t>>(),
             py::arg("contexts"),
             py::arg("channels"),
             py::arg("status_file") = "",
@@ -1720,9 +1745,10 @@ PYBIND11_MODULE(_madspace_py, m) {
                 "config",
                 EventGenerator::default_config,
                 "EventGenerator.default_config"
-            )
+            ),
+            py::arg("seed") = std::nullopt
         )
-        .def("survey", &EventGenerator::survey)
+        .def("survey", &EventGenerator::survey, py::arg("survey_pass") = 0)
         .def("generate", &EventGenerator::generate)
         .def(
             "combine_to_compact_npy",
diff --git a/madspace/tests/test_lhe.py b/madspace/tests/test_lhe.py
index 46d3d0d1a..f5338c2a2 100644
--- a/madspace/tests/test_lhe.py
+++ b/madspace/tests/test_lhe.py
@@ -148,7 +148,7 @@ def external_particles(event):
 @pytest.fixture(scope="module")
 def events(lhe_completer, mapping):
     p_ext = sample_external_momenta(mapping, 300, seed=1234)
-    rand_gen = ms.RandGen(2024)
+    rand_gen = ms.MixMaxRandom(2024)
     result = []
     for row in p_ext:
         event = build_event(row)
@@ -309,7 +309,7 @@ def test_external_color_flows_match_input_across_subprocesses(
     _, _, subproc_args = topology_and_args
     lhe_completer = ms.LHECompleter([subproc_args, subproc_args], bw_cutoff=BW_CUTOFF)
     p_ext = sample_external_momenta(mapping, 20, seed=99)
-    rand_gen = ms.RandGen(7)
+    rand_gen = ms.MixMaxRandom(7)
     color_flows = subproc_meta["color_flows"]
 
     for subprocess_index in (0, 1):
@@ -329,7 +329,7 @@ def test_external_color_flows_match_input_across_subprocesses(
 def test_external_spins_match_helicity_table(lhe_completer, mapping, subproc_meta):
     helicities = subproc_meta["helicities"]
     p_ext = sample_external_momenta(mapping, 5, seed=55)
-    rand_gen = ms.RandGen(3)
+    rand_gen = ms.MixMaxRandom(3)
     for helicity_index in [0, 3, 10]:
         for row in p_ext:
             event = build_event(row)
@@ -343,7 +343,7 @@ def test_external_spins_match_helicity_table(lhe_completer, mapping, subproc_met
 
 def test_external_flavors_are_valid_options(lhe_completer, mapping, subproc_meta):
     p_ext = sample_external_momenta(mapping, 30, seed=77)
-    rand_gen = ms.RandGen(9)
+    rand_gen = ms.MixMaxRandom(9)
     for flavor_index, flavor in enumerate(subproc_meta["flavors"]):
         options = [tuple(option) for option in flavor["options"]]
         for row in p_ext:
@@ -363,8 +363,10 @@ def test_save_load_roundtrip(lhe_completer, mapping):
         p_ext = sample_external_momenta(mapping, 20, seed=321)
         for row in p_ext:
             event_a, event_b = build_event(row), build_event(row)
-            lhe_completer.complete_event_data(event_a, 0, 0, 0, 0, 0, ms.RandGen(11))
-            loaded.complete_event_data(event_b, 0, 0, 0, 0, 0, ms.RandGen(11))
+            lhe_completer.complete_event_data(
+                event_a, 0, 0, 0, 0, 0, ms.MixMaxRandom(11)
+            )
+            loaded.complete_event_data(event_b, 0, 0, 0, 0, 0, ms.MixMaxRandom(11))
             assert len(event_a.particles) == len(event_b.particles)
             for pa, pb in zip(event_a.particles, event_b.particles):
                 assert pa.pdg_id == pb.pdg_id
@@ -379,11 +381,11 @@ def test_wrong_particle_count_raises(lhe_completer):
     event = ms.LHEEvent()
     event.particles = [ms.LHEParticle()] * 3
     with pytest.raises(RuntimeError):
-        lhe_completer.complete_event_data(event, 0, 0, 0, 0, 0, ms.RandGen(1))
+        lhe_completer.complete_event_data(event, 0, 0, 0, 0, 0, ms.MixMaxRandom(1))
 
 
 def test_invalid_color_index_raises(lhe_completer, mapping):
     p_ext = sample_external_momenta(mapping, 1, seed=1)
     event = build_event(p_ext[0])
     with pytest.raises(RuntimeError):
-        lhe_completer.complete_event_data(event, 0, 0, 999, 0, 0, ms.RandGen(1))
+        lhe_completer.complete_event_data(event, 0, 0, 999, 0, 0, ms.MixMaxRandom(1))
diff --git a/tests/acceptance_tests/test_mg7_reproducibility.py b/tests/acceptance_tests/test_mg7_reproducibility.py
new file mode 100644
index 000000000..6f9edd8e3
--- /dev/null
+++ b/tests/acceptance_tests/test_mg7_reproducibility.py
@@ -0,0 +1,379 @@
+################################################################################
+#
+# Copyright (c) 2009 The MadGraph5_aMC@NLO Development team and Contributors
+#
+# This file is a part of the MadGraph5_aMC@NLO project, an application which
+# automatically generates Feynman diagrams and matrix elements for arbitrary
+# high-energy processes in the Standard Model and beyond.
+#
+# It is subject to the MadGraph5_aMC@NLO license which should accompany this
+# distribution.
+#
+# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
+#
+################################################################################
+"""MG7 (madspace) reproducibility tests.
+
+Guards the seeded event-generation path: given the same run_card ``seed``,
+independent mg7 runs must produce byte-identical LHE output (hashed with
+sha256); a different seed must change the result. Covers:
+
+  * ``test_vegas_reproducibility_mg7`` -- a plain survey (VEGAS-optimized)
+    + generate run, through ``bin/generate_events``.
+  * ``test_madnis_reproducible_mode_mg7`` -- a full survey + madnis training +
+    generate run, through ``bin/generate_events`` (no gridpack export), with
+    ``madnis.reproducible = True``: training is required to be
+    thread-count-independent and byte-identical for the same seed, both with
+    online-only training and with buffered (off-policy replay) training
+    enabled. GPU multi-channel batches are not covered by this mode.
+  * ``test_gridpack_reproducibility_mg7`` -- a gridpack trained with madnis,
+    then standalone event generation from that gridpack (its own
+    ``bin/generate_events --seed``), which is the normal way a gridpack is
+    used and does not re-run survey/training.
+  * ``test_gridpack_reproducibility_vegas_mg7`` -- same as above, but for a
+    plain VEGAS-optimized gridpack (no madnis training).
+
+Process: ``p p > t t~`` (hadronic, needs the NNPDF23_lo_as_0130_qed PDF grid;
+self-skips if the mg7 runtime stack / PDF is unavailable, same as
+test_check_xsec_processes_mg7.py).
+
+Run locally with e.g.::
+
+    ./tests/test_manager.py test_.*reproducib.*_mg7 -pA -t0 -l INFO
+
+One knob is read from the environment so the CI can dial it without touching
+the code:
+
+  * ``MG7_REPRO_EVENTS`` -- events per run (default 10000).
+"""
+
+from __future__ import absolute_import
+from __future__ import division
+
+import glob
+import hashlib
+import os
+import shutil
+import subprocess
+import sys
+import tempfile
+import unittest
+
+import madgraph.interface.master_interface as MGCmd
+from madgraph.various.banner import RunCardMG7
+
+pjoin = os.path.join
+
+_PROCESS = 'p p > t t~'
+_EVENTS = int(os.environ.get('MG7_REPRO_EVENTS', 10000))
+_SEED_A = 424242
+_SEED_B = 909090
+
+
+def _mg7_datadir_or_skip(test):
+    """Return an LHAPDF data dir that contains the NNPDF23_lo_as_0130_qed set,
+    or ``skipTest`` (on *test*) when the mg7 runtime stack (madspace + LHAPDF +
+    the run_card.toml default PDF) is unavailable."""
+    try:
+        import madspace
+        has_mg7 = hasattr(madspace, 'ChannelEventGenerator')
+    except ImportError:
+        has_mg7 = False
+    if not has_mg7:
+        test.skipTest('mg7 runtime stack (madspace) unavailable')
+
+    candidates = []
+    if os.environ.get('LHAPDF_DATA_PATH'):
+        candidates.extend(os.environ['LHAPDF_DATA_PATH'].split(os.pathsep))
+    try:
+        out = subprocess.check_output(['lhapdf-config', '--datadir'],
+                                      stderr=subprocess.DEVNULL).decode().strip()
+        if out:
+            candidates.append(out)
+    except Exception:
+        pass
+    for d in candidates:
+        if d and os.path.isdir(d) and glob.glob(pjoin(d, 'NNPDF23_lo_as_0130_qed*')):
+            return d
+    test.skipTest('NNPDF23_lo_as_0130_qed LHAPDF data not found '
+                  '(set $LHAPDF_DATA_PATH)')
+
+
+def _lhe_hash(path):
+    """sha256 of an LHE file's physics content: everything from ``<init>``
+    onward, i.e. excluding the ``<header>`` block. The header embeds the
+    run_card.toml/param_card/proc_card verbatim, so it legitimately differs
+    between runs that vary settings unrelated to the physics content (e.g.
+    the seed itself, or the cpu thread pool size) -- only the generated
+    events (and the <init> cross-section info, itself a deterministic
+    function of the seed) should drive this hash."""
+    with open(path, 'rb') as f:
+        content = f.read()
+    marker = b'</header>\n'
+    index = content.find(marker)
+    if index == -1:
+        raise AssertionError('no </header> found in %s' % path)
+    return hashlib.sha256(content[index + len(marker):]).hexdigest()
+
+
+def _find_lhe(run_path):
+    """Locate the LHE file produced under Events/<run>/, or fail."""
+    matches = sorted(glob.glob(pjoin(run_path, 'Events', '*', 'events.lhe')))
+    if not matches:
+        raise AssertionError('no events.lhe produced under %s' % run_path)
+    return matches[-1]
+
+
+def _run(cmd, cwd, log_path, env, what):
+    """Run *cmd*, raising with the log tail on a non-zero exit."""
+    with open(log_path, 'w') as logfh:
+        ret = subprocess.call(cmd, cwd=cwd, env=env, stdout=logfh,
+                              stderr=subprocess.STDOUT)
+    if ret != 0:
+        with open(log_path) as f:
+            tail = ''.join(f.readlines()[-60:])
+        raise AssertionError('%s failed (exit %d, see %s)\n\n%s'
+                             % (what, ret, log_path, tail))
+
+
+class MG7ReproducibilityTest(unittest.TestCase):
+    """LHE-hash reproducibility of mg7 (madspace) event generation for a
+    fixed run_card seed, and non-reproducibility across distinct seeds."""
+
+    def setUp(self):
+        self.path = tempfile.mkdtemp(prefix='mg7_repro_')
+
+    def tearDown(self):
+        shutil.rmtree(self.path, ignore_errors=True)
+
+    def _output_process(self, run_name):
+        run_dir = pjoin(self.path, run_name)
+        mg = MGCmd.MasterCmd()
+        mg.no_notification()
+        mg.exec_cmd('set automatic_html_opening False --no_save')
+        mg.exec_cmd('generate %s' % _PROCESS)
+        mg.exec_cmd('output mg7 %s' % run_dir)
+        return run_dir
+
+    def _set_run_card(self, toml_path, **settings):
+        """settings keys are 'section.key' RunCardMG7 names."""
+        rc = RunCardMG7(toml_path)
+        for key, value in settings.items():
+            rc.set(key, value, user=True)
+        rc.write(toml_path)
+
+    def _generate_and_hash(self, run_dir, datadir, seed, tag, **run_card_settings):
+        """Set run_card seed (plus any extra 'section.key' settings), run
+        bin/generate_events -f from scratch, and return the sha256 of the
+        resulting LHE file."""
+        shutil.rmtree(pjoin(run_dir, 'Events'), ignore_errors=True)
+        toml = pjoin(run_dir, 'Cards', 'run_card.toml')
+        self._set_run_card(toml, **{'run.seed': seed}, **run_card_settings)
+        env = dict(os.environ)
+        env['LHAPDF_DATA_PATH'] = datadir
+        _run([sys.executable, pjoin(run_dir, 'bin', 'generate_events'), '-f'],
+             run_dir, pjoin(run_dir, 'gen_%s.log' % tag), env,
+             'mg7 generate_events')
+        return _lhe_hash(_find_lhe(run_dir))
+
+    def _generate_and_hash_gridpack(self, gridpack_dir, datadir, seed, tag):
+        """Run the gridpack's own bin/generate_events --seed from scratch, and
+        return the sha256 of the resulting LHE file."""
+        shutil.rmtree(pjoin(gridpack_dir, 'Events'), ignore_errors=True)
+        env = dict(os.environ)
+        env['LHAPDF_DATA_PATH'] = datadir
+        _run([sys.executable, pjoin(gridpack_dir, 'bin', 'generate_events'),
+              '--seed', str(seed), '--events', str(_EVENTS),
+              '--output_format', 'lhe'],
+             gridpack_dir, pjoin(gridpack_dir, 'gen_%s.log' % tag), env,
+             'gridpack generate_events')
+        return _lhe_hash(_find_lhe(gridpack_dir))
+
+    def test_vegas_reproducibility_mg7(self):
+        """Two VEGAS-optimized runs with the same seed produce byte-identical
+        LHE files, including when the cpu thread pool is shrunk to 1; a
+        different seed changes the result."""
+        datadir = _mg7_datadir_or_skip(self)
+        run_dir = self._output_process('vegas')
+        toml = pjoin(run_dir, 'Cards', 'run_card.toml')
+        self._set_run_card(
+            toml,
+            **{
+                'generation.events': _EVENTS,
+                # Keep the test scoped to madspace's own seeding: LHE-level
+                # post-processing (systematics) is orthogonal and would only
+                # add runtime here.
+                'postprocessing.systematics': False,
+            }
+        )
+
+        # Every call sets 'run.cpu_thread_pool_size' explicitly so it never
+        # leaks from a previous call via the run_card left on disk.
+        hash_a1 = self._generate_and_hash(
+            run_dir, datadir, _SEED_A, 'a1', **{'run.cpu_thread_pool_size': -1})
+        hash_a2 = self._generate_and_hash(
+            run_dir, datadir, _SEED_A, 'a2', **{'run.cpu_thread_pool_size': -1})
+        self.assertEqual(hash_a1, hash_a2,
+                         'same seed produced different LHE content')
+
+        hash_a_single = self._generate_and_hash(
+            run_dir, datadir, _SEED_A, 'a_single',
+            **{'run.cpu_thread_pool_size': 1})
+        self.assertEqual(hash_a1, hash_a_single,
+                         'same seed produced different LHE content with a '
+                         'single-threaded cpu thread pool')
+
+        hash_b = self._generate_and_hash(
+            run_dir, datadir, _SEED_B, 'b', **{'run.cpu_thread_pool_size': -1})
+        self.assertNotEqual(hash_a1, hash_b,
+                            'different seeds produced identical LHE content')
+
+    def test_madnis_reproducible_mode_mg7(self):
+        """With madnis.reproducible = True (CPU codepath), madnis training
+        itself -- not just event generation -- is byte-identical for the same
+        seed independent of the cpu thread pool size, both with online-only
+        training and with buffered (off-policy replay) training enabled.
+        Covers the scheduling fix in MadnisTraining::maybe_start_generator_jobs
+        / start_generator_jobs (deterministic round dispatch + deferred buffer
+        flush) and the seeding of BufferUnweighter/BatchSampler (see
+        salt::madnis_train_unweight / salt::buffered_batch_sample in
+        random.hpp) -- both required for this test to pass reliably."""
+        datadir = _mg7_datadir_or_skip(self)
+        run_dir = self._output_process('madnis_reproducible')
+        toml = pjoin(run_dir, 'Cards', 'run_card.toml')
+        base_settings = {
+            'generation.events': _EVENTS,
+            'postprocessing.systematics': False,
+            'madnis.enable': True,
+            'madnis.train_batches': 60,
+            'madnis.reproducible': True,
+        }
+
+        # Online-only training (buffering disabled).
+        self._set_run_card(toml, **base_settings, **{'madnis.buffer_capacity': 0})
+        hash_online_multi = self._generate_and_hash(
+            run_dir, datadir, _SEED_A, 'online_multi',
+            **{'run.cpu_thread_pool_size': -1})
+        hash_online_single = self._generate_and_hash(
+            run_dir, datadir, _SEED_A, 'online_single',
+            **{'run.cpu_thread_pool_size': 1})
+        self.assertEqual(
+            hash_online_multi, hash_online_single,
+            'madnis.reproducible training (buffering disabled) produced '
+            'different LHE content with a single-threaded cpu thread pool')
+
+        # Buffered (off-policy replay) training enabled.
+        self._set_run_card(
+            toml,
+            **base_settings,
+            **{
+                'madnis.buffer_capacity': 3000,
+                'madnis.buffered_steps': 3,
+                'madnis.minimum_buffer_size': 500,
+            }
+        )
+        hash_buffered_multi = self._generate_and_hash(
+            run_dir, datadir, _SEED_A, 'buffered_multi',
+            **{'run.cpu_thread_pool_size': -1})
+        hash_buffered_single = self._generate_and_hash(
+            run_dir, datadir, _SEED_A, 'buffered_single',
+            **{'run.cpu_thread_pool_size': 1})
+        self.assertEqual(
+            hash_buffered_multi, hash_buffered_single,
+            'madnis.reproducible training (buffering enabled) produced '
+            'different LHE content with a single-threaded cpu thread pool')
+
+        hash_buffered_b = self._generate_and_hash(
+            run_dir, datadir, _SEED_B, 'buffered_b',
+            **{'run.cpu_thread_pool_size': -1})
+        self.assertNotEqual(
+            hash_buffered_multi, hash_buffered_b,
+            'different seeds produced identical LHE content')
+
+    def test_gridpack_reproducibility_mg7(self):
+        """A gridpack trained with madnis: two standalone event-generation
+        runs from that gridpack with the same seed produce byte-identical LHE
+        files; a different seed changes the result. Training itself is not
+        required to be deterministic -- only the trained gridpack's event
+        generation is under test here."""
+        datadir = _mg7_datadir_or_skip(self)
+        run_dir = self._output_process('gridpack')
+        toml = pjoin(run_dir, 'Cards', 'run_card.toml')
+        self._set_run_card(
+            toml,
+            **{
+                'generation.events': _EVENTS,
+                'postprocessing.systematics': False,
+                'gridpack.save_gridpack': True,
+                'madnis.enable': True,
+                'madnis.train_batches': 100,
+            }
+        )
+
+        env = dict(os.environ)
+        env['LHAPDF_DATA_PATH'] = datadir
+        _run([sys.executable, pjoin(run_dir, 'bin', 'generate_events'), '-f'],
+             run_dir, pjoin(run_dir, 'train.log'), env,
+             'mg7 madnis training run')
+
+        gridpack_dir = pjoin(run_dir, 'Events', 'run_01', 'gridpack')
+        self.assertTrue(os.path.isdir(gridpack_dir),
+                        'gridpack was not produced under %s' % run_dir)
+
+        hash_a1 = self._generate_and_hash_gridpack(
+            gridpack_dir, datadir, _SEED_A, 'a1')
+        hash_a2 = self._generate_and_hash_gridpack(
+            gridpack_dir, datadir, _SEED_A, 'a2')
+        self.assertEqual(hash_a1, hash_a2,
+                         'same seed produced different LHE content')
+
+        hash_b = self._generate_and_hash_gridpack(
+            gridpack_dir, datadir, _SEED_B, 'b')
+        self.assertNotEqual(hash_a1, hash_b,
+                            'different seeds produced identical LHE content')
+
+    def test_gridpack_reproducibility_vegas_mg7(self):
+        """A plain VEGAS-optimized gridpack (no madnis training): two
+        standalone event-generation runs from that gridpack with the same
+        seed produce byte-identical LHE files; a different seed changes the
+        result. Companion to test_gridpack_reproducibility_mg7, which covers
+        the madnis-trained case; the survey/optimization itself is not
+        required to be deterministic here either -- only the gridpack's own
+        event generation is under test."""
+        datadir = _mg7_datadir_or_skip(self)
+        run_dir = self._output_process('gridpack_vegas')
+        toml = pjoin(run_dir, 'Cards', 'run_card.toml')
+        self._set_run_card(
+            toml,
+            **{
+                'generation.events': _EVENTS,
+                'postprocessing.systematics': False,
+                'gridpack.save_gridpack': True,
+            }
+        )
+
+        env = dict(os.environ)
+        env['LHAPDF_DATA_PATH'] = datadir
+        _run([sys.executable, pjoin(run_dir, 'bin', 'generate_events'), '-f'],
+             run_dir, pjoin(run_dir, 'survey.log'), env,
+             'mg7 vegas survey run')
+
+        gridpack_dir = pjoin(run_dir, 'Events', 'run_01', 'gridpack')
+        self.assertTrue(os.path.isdir(gridpack_dir),
+                        'gridpack was not produced under %s' % run_dir)
+
+        hash_a1 = self._generate_and_hash_gridpack(
+            gridpack_dir, datadir, _SEED_A, 'a1')
+        hash_a2 = self._generate_and_hash_gridpack(
+            gridpack_dir, datadir, _SEED_A, 'a2')
+        self.assertEqual(hash_a1, hash_a2,
+                         'same seed produced different LHE content')
+
+        hash_b = self._generate_and_hash_gridpack(
+            gridpack_dir, datadir, _SEED_B, 'b')
+        self.assertNotEqual(hash_a1, hash_b,
+                            'different seeds produced identical LHE content')
+
+
+if __name__ == '__main__':
+    unittest.main()