From 7af1fc7600d8b8838012b551b48bf9e29d170c24 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 12 Aug 2026 03:57:16 -0700
Subject: [PATCH 01/30] Enable free-threaded CPython support

---
 mkl_random/mklrand.pyx | 1 +
 pyproject.toml         | 3 ++-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index edbcd400..215ae25d 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -24,6 +24,7 @@
 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 # cython: language_level=3
+# cython: freethreading_compatible=True
 
 cdef extern from "Python.h":
     void* PyMem_Malloc(size_t n)
diff --git a/pyproject.toml b/pyproject.toml
index c841c9a3..983c3601 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -28,7 +28,7 @@ build-backend = "mesonpy"
 requires = [
   "meson-python>=0.16.0",
   "ninja",
-  "Cython",
+  "Cython>=3.1.0",
   "numpy",
   "mkl-devel",
   "cmake"
@@ -50,6 +50,7 @@ classifiers = [
   "Programming Language :: Python :: 3.12",
   "Programming Language :: Python :: 3.13",
   "Programming Language :: Python :: 3.14",
+  "Programming Language :: Python :: Free Threading :: 2 - Beta",
   "Programming Language :: Python :: Implementation :: CPython",
   "Topic :: Software Development",
   "Topic :: Scientific/Engineering",

From ac81826eebc16006d03a7f64f8c0fa858b90d8dc Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 12 Aug 2026 03:58:02 -0700
Subject: [PATCH 02/30] Update conda recipes for free-threaded builds

---
 conda-recipe-cf/meta.yaml | 4 +---
 conda-recipe/meta.yaml    | 4 +---
 2 files changed, 2 insertions(+), 6 deletions(-)

diff --git a/conda-recipe-cf/meta.yaml b/conda-recipe-cf/meta.yaml
index c1dcd498..96b019f5 100644
--- a/conda-recipe-cf/meta.yaml
+++ b/conda-recipe-cf/meta.yaml
@@ -21,15 +21,13 @@ requirements:
       - cmake
       - ninja
       - python
-      - python-gil     # [py>=314]
       - mkl-devel
-      - cython
+      - cython >=3.1.0
       - numpy
       - pip
       - wheel  >=0.41.3
     run:
       - python
-      - python-gil     # [py>=314]
       - {{ pin_compatible('numpy', min_pin="x.x", max_pin="x") }}
       - {{ pin_compatible('mkl', min_pin="x.x", max_pin="x") }}
 
diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml
index dbdedbbe..03d03ca0 100644
--- a/conda-recipe/meta.yaml
+++ b/conda-recipe/meta.yaml
@@ -24,15 +24,13 @@ requirements:
       - cmake
       - ninja
       - python
-      - python-gil     # [py>=314]
       - mkl-devel
-      - cython
+      - cython >=3.1.0
       - numpy
       - pip
       - wheel  >=0.41.3
     run:
       - python
-      - python-gil     # [py>=314]
       - {{ pin_compatible('mkl', min_pin="x.x", max_pin="x") }}
       - numpy >=1.26.4
 

From dd3c9c99756bfb25b435cc1a043ee0ff5e3a8ee1 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 12 Aug 2026 04:13:05 -0700
Subject: [PATCH 03/30] Extend GH workflow matrices with free-threaded Python
 3.14t

---
 .github/workflows/build-with-clang.yml        |  4 +-
 .../workflows/build-with-standard-clang.yml   |  4 +-
 .github/workflows/build_pip.yml               | 19 +++-
 .github/workflows/conda-package-cf.yml        | 90 +++++++++++++------
 .github/workflows/conda-package.yml           | 74 ++++++++++-----
 5 files changed, 137 insertions(+), 54 deletions(-)

diff --git a/.github/workflows/build-with-clang.yml b/.github/workflows/build-with-clang.yml
index 45137e15..c3a58aca 100644
--- a/.github/workflows/build-with-clang.yml
+++ b/.github/workflows/build-with-clang.yml
@@ -12,7 +12,7 @@ jobs:
 
     strategy:
       matrix:
-        python: ["3.10", "3.11", "3.12", "3.13", "3.14"]
+        python: ["3.10", "3.11", "3.12", "3.13", "3.14", "3.14t"]
         numpy_version: ["'numpy>=2'"]
     env:
       ONEAPI_ROOT: /opt/intel/oneapi
@@ -53,7 +53,7 @@ jobs:
 
       - name: Install mkl_random dependencies
         run: |
-          pip install meson-python ninja cython cmake
+          pip install meson-python ninja "cython>=3.1.0" cmake
           pip install ${{ matrix.numpy_version }}
 
       - name: List oneAPI folder content
diff --git a/.github/workflows/build-with-standard-clang.yml b/.github/workflows/build-with-standard-clang.yml
index 97c1d48b..39b59501 100644
--- a/.github/workflows/build-with-standard-clang.yml
+++ b/.github/workflows/build-with-standard-clang.yml
@@ -13,7 +13,7 @@ jobs:
 
     strategy:
       matrix:
-        python: ["3.10", "3.11", "3.12", "3.13", "3.14"]
+        python: ["3.10", "3.11", "3.12", "3.13", "3.14", "3.14t"]
         numpy_version: ["'numpy>=2'"]
 
     env:
@@ -47,7 +47,7 @@ jobs:
 
       - name: Install mkl_random dependencies
         run: |
-          pip install meson-python ninja cmake cython mkl-devel
+          pip install meson-python ninja cmake "cython>=3.1.0" mkl-devel
           pip install ${{ matrix.numpy_version }}
 
       - name: Build mkl_random
diff --git a/.github/workflows/build_pip.yml b/.github/workflows/build_pip.yml
index 798b1ff7..43062e4d 100644
--- a/.github/workflows/build_pip.yml
+++ b/.github/workflows/build_pip.yml
@@ -17,8 +17,21 @@ jobs:
 
     strategy:
       matrix:
-        python: ["3.10", "3.11", "3.12", "3.13", "3.14"]
+        python: ["3.10", "3.11", "3.12", "3.13"]
         use_pre: ["", "--pre"]
+        include:
+          - python: "3.14"
+            python_spec: "3.14.* *_cp314"
+            use_pre: ""
+          - python: "3.14"
+            python_spec: "3.14.* *_cp314"
+            use_pre: "--pre"
+          - python: "3.14"
+            python_spec: "3.14.* *_cp314t"
+            use_pre: ""
+          - python: "3.14"
+            python_spec: "3.14.* *_cp314t"
+            use_pre: "--pre"
 
     steps:
       - uses: actions/checkout@de0fac2e4500dabe0009e67214ff5f5447ce83dd # v6.0.2
@@ -34,11 +47,11 @@ jobs:
 
       - name: Install MKL
         run: |
-          conda install mkl-devel
+          conda install mkl-devel "python=${{ matrix.python_spec || matrix.python }}"
 
       - name: Build with pip
         run: |
-          pip install --no-cache-dir meson-python ninja cmake cython
+          pip install --no-cache-dir meson-python ninja cmake "cython>=3.1.0"
           pip install --no-cache-dir numpy ${{ matrix.use_pre }}
           pip install -e . --no-build-isolation --no-deps --verbose
           pip install --no-cache-dir pytest
diff --git a/.github/workflows/conda-package-cf.yml b/.github/workflows/conda-package-cf.yml
index 0fb666eb..516e1b09 100644
--- a/.github/workflows/conda-package-cf.yml
+++ b/.github/workflows/conda-package-cf.yml
@@ -30,6 +30,12 @@ jobs:
             numpy: "2.3"
           - python: "3.14"
             numpy: "2.3"
+            python_spec: "3.14.* *_cp314"
+            python_tag: "3.14"
+          - python: "3.14"
+            numpy: "2.3"
+            python_spec: "3.14.* *_cp314t"
+            python_tag: "3.14t"
 
     steps:
       - name: Cancel Previous Runs
@@ -52,9 +58,9 @@ jobs:
         with:
           path: ~/.conda/pkgs
           key:
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-${{hashFiles('**/meta.yaml') }}
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-${{hashFiles('**/meta.yaml') }}
           restore-keys: |
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-
             ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-
 
       - name: Add conda to system path
@@ -66,7 +72,7 @@ jobs:
       - name: Build conda package
         run: |
           CHANNELS=(-c conda-forge --override-channels)
-          VERSIONS=(--python "${{ matrix.python }}" --numpy "${{ matrix.numpy }}")
+          VERSIONS=(--python "${{ matrix.python_spec || matrix.python }}" --numpy "${{ matrix.numpy }}")
           TEST=(--no-test)
 
           conda build \
@@ -78,7 +84,7 @@ jobs:
       - name: Upload artifact
         uses: actions/upload-artifact@043fb46d1a93c77aae656e7c1c64a875d1fc6a0a # v7.0.1
         with:
-          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python }}
+          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python_tag || matrix.python }}
           path: /usr/share/miniconda/conda-bld/linux-64/${{ env.PACKAGE_NAME }}-*.conda
 
   build_windows:
@@ -99,6 +105,12 @@ jobs:
             numpy: "2.3"
           - python: "3.14"
             numpy: "2.3"
+            python_spec: "3.14.* *_cp314"
+            python_tag: "3.14"
+          - python: "3.14"
+            numpy: "2.3"
+            python_spec: "3.14.* *_cp314t"
+            python_tag: "3.14t"
 
     env:
       conda-bld: C:\Miniconda\conda-bld\win-64\
@@ -127,9 +139,9 @@ jobs:
         with:
           path: /home/runner/conda_pkgs_dir
           key:
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-${{hashFiles('**/meta.yaml') }}
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-${{hashFiles('**/meta.yaml') }}
           restore-keys: |
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-
             ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-
 
       - name: Store conda paths as envs
@@ -143,12 +155,12 @@ jobs:
           conda list -n base
 
       - name: Build conda package
-        run: conda build --no-test --python ${{ matrix.python }} --numpy ${{ matrix.numpy }} -c conda-forge --override-channels conda-recipe-cf
+        run: conda build --no-test --python "${{ matrix.python_spec || matrix.python }}" --numpy ${{ matrix.numpy }} -c conda-forge --override-channels conda-recipe-cf
 
       - name: Upload artifact
         uses: actions/upload-artifact@043fb46d1a93c77aae656e7c1c64a875d1fc6a0a # v7.0.1
         with:
-          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python }}
+          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python_tag || matrix.python }}
           path: ${{ env.CONDA_BLD }}${{ env.PACKAGE_NAME }}-*.conda
 
   test_linux:
@@ -168,6 +180,12 @@ jobs:
             numpy: "2.3"
           - python: "3.14"
             numpy: "2.3"
+            python_spec: "3.14.* *_cp314"
+            python_tag: "3.14"
+          - python: "3.14"
+            numpy: "2.3"
+            python_spec: "3.14.* *_cp314t"
+            python_tag: "3.14t"
 
     env:
       CHANNELS: -c conda-forge --override-channels
@@ -176,7 +194,7 @@ jobs:
       - name: Download artifact
         uses: actions/download-artifact@3e5f45b2cfb9172054b4087a40e8e0b5a5461e7c # v8.0.1
         with:
-          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python }}
+          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python_tag || matrix.python }}
 
       - name: Add conda to system path
         run: echo "$CONDA/bin" >> "$GITHUB_PATH"
@@ -199,7 +217,7 @@ jobs:
           . "$CONDA/etc/profile.d/conda.sh"
           PACKAGE_VERSION="$(python -c "${VER_SCRIPT1} ${VER_SCRIPT2}")"
           export PACKAGE_VERSION
-          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" "python=${{ matrix.python }}" "numpy=${{ matrix.numpy }}" -c "$GITHUB_WORKSPACE/channel" -c conda-forge --override-channels --only-deps --dry-run > lockfile
+          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" "python=${{ matrix.python_spec || matrix.python }}" "numpy=${{ matrix.numpy }}" -c "$GITHUB_WORKSPACE/channel" -c conda-forge --override-channels --only-deps --dry-run > lockfile
           cat lockfile
 
       - name: Set pkgs_dirs
@@ -213,9 +231,9 @@ jobs:
         with:
           path: ~/.conda/pkgs
           key:
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-${{hashFiles('lockfile') }}
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-${{hashFiles('lockfile') }}
           restore-keys: |
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-
             ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-
 
       - name: Install mkl_random
@@ -223,7 +241,7 @@ jobs:
           . "$CONDA/etc/profile.d/conda.sh"
           PACKAGE_VERSION="$(python -c "${VER_SCRIPT1} ${VER_SCRIPT2}")"
           export PACKAGE_VERSION
-          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" pytest "python=${{ matrix.python }}" "numpy=${{ matrix.numpy }}" -c "$GITHUB_WORKSPACE/channel" -c conda-forge --override-channels
+          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" pytest "python=${{ matrix.python_spec || matrix.python }}" "numpy=${{ matrix.numpy }}" -c "$GITHUB_WORKSPACE/channel" -c conda-forge --override-channels
           # Test installed packages
           conda list
 
@@ -252,6 +270,12 @@ jobs:
             numpy: "2.3"
           - python: "3.14"
             numpy: "2.3"
+            python_spec: "3.14.* *_cp314"
+            python_tag: "3.14"
+          - python: "3.14"
+            numpy: "2.3"
+            python_spec: "3.14.* *_cp314t"
+            python_tag: "3.14t"
 
     env:
       CHANNELS: -c conda-forge --override-channels
@@ -260,7 +284,7 @@ jobs:
       - name: Download artifact
         uses: actions/download-artifact@3e5f45b2cfb9172054b4087a40e8e0b5a5461e7c # v8.0.1
         with:
-          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python }}
+          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python_tag || matrix.python }}
 
       - uses: conda-incubator/setup-miniconda@8ee1f361103df19b6f8c8655fd3967a8ecb162d5 # v4.0.1
         with:
@@ -292,7 +316,7 @@ jobs:
           FOR /F "tokens=* USEBACKQ" %%F IN (`python -c "%SCRIPT%"`) DO (
              SET PACKAGE_VERSION=%%F
           )
-          conda install -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% python=${{ matrix.python }} numpy=${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }} --only-deps --dry-run > lockfile
+          conda install -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% "python=${{ matrix.python_spec || matrix.python }}" numpy=${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }} --only-deps --dry-run > lockfile
 
       - name: Cache conda packages
         uses: actions/cache@55cc8345863c7cc4c66a329aec7e433d2d1c52a9 # v6.1.0
@@ -301,9 +325,9 @@ jobs:
         with:
           path: /home/runner/conda_pkgs_dir
           key:
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-${{hashFiles('lockfile') }}
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-${{hashFiles('lockfile') }}
           restore-keys: |
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-
             ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-
 
       - name: Install mkl_random
@@ -315,7 +339,7 @@ jobs:
           FOR /F "tokens=* USEBACKQ" %%F IN (`python -c "%SCRIPT%"`) DO (
              SET PACKAGE_VERSION=%%F
           )
-          conda create -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% pytest python=${{ matrix.python }} numpy=${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }}
+          conda create -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% pytest "python=${{ matrix.python_spec || matrix.python }}" numpy=${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }}
           # Test installed packages
           conda list
 
@@ -339,6 +363,12 @@ jobs:
             numpy: "2.3"
           - python: "3.14"
             numpy: "2.3"
+            python_spec: "3.14.* *_cp314"
+            python_tag: "3.14"
+          - python: "3.14"
+            numpy: "2.3"
+            python_spec: "3.14.* *_cp314t"
+            python_tag: "3.14t"
 
     steps:
       - name: Cancel Previous Runs
@@ -364,9 +394,9 @@ jobs:
         with:
           path: /Users/runner/conda_pkgs_dir
           key:
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-${{hashFiles('**/meta.yaml') }}
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-${{hashFiles('**/meta.yaml') }}
           restore-keys: |
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-
             ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-
 
       - name: Install conda-build
@@ -384,7 +414,7 @@ jobs:
         shell: bash -el {0}
         run: |
           CHANNELS=(-c conda-forge --override-channels)
-          VERSIONS=(--python "${{ matrix.python }}" --numpy "${{ matrix.numpy }}")
+          VERSIONS=(--python "${{ matrix.python_spec || matrix.python }}" --numpy "${{ matrix.numpy }}")
           TEST=(--no-test)
 
           conda build \
@@ -396,7 +426,7 @@ jobs:
       - name: Upload artifact
         uses: actions/upload-artifact@043fb46d1a93c77aae656e7c1c64a875d1fc6a0a # v7.0.1
         with:
-          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python }}
+          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python_tag || matrix.python }}
           path: ${{ env.CONDA_BLD }}${{ env.PACKAGE_NAME }}-*.conda
 
   test_osx:
@@ -415,12 +445,18 @@ jobs:
             numpy: "2.3"
           - python: "3.14"
             numpy: "2.3"
+            python_spec: "3.14.* *_cp314"
+            python_tag: "3.14"
+          - python: "3.14"
+            numpy: "2.3"
+            python_spec: "3.14.* *_cp314t"
+            python_tag: "3.14t"
 
     steps:
       - name: Download artifact
         uses: actions/download-artifact@3e5f45b2cfb9172054b4087a40e8e0b5a5461e7c # v8.0.1
         with:
-          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python }}
+          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python_tag || matrix.python }}
 
       - uses: conda-incubator/setup-miniconda@8ee1f361103df19b6f8c8655fd3967a8ecb162d5 # v4.0.1
         with:
@@ -448,7 +484,7 @@ jobs:
         run: |
           PACKAGE_VERSION="$(python -c "${VER_SCRIPT1} ${VER_SCRIPT2}")"
           export PACKAGE_VERSION
-          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" "python=${{ matrix.python }}" "numpy=${{ matrix.numpy }}" -c "$GITHUB_WORKSPACE/channel" -c conda-forge --override-channels --only-deps --dry-run > lockfile
+          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" "python=${{ matrix.python_spec || matrix.python }}" "numpy=${{ matrix.numpy }}" -c "$GITHUB_WORKSPACE/channel" -c conda-forge --override-channels --only-deps --dry-run > lockfile
           cat lockfile
 
       - name: Cache conda packages
@@ -458,9 +494,9 @@ jobs:
         with:
           path: /Users/runner/conda_pkgs_dir
           key:
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-${{hashFiles('lockfile') }}
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-${{hashFiles('lockfile') }}
           restore-keys: |
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-
             ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-
 
       - name: Install mkl_random
@@ -468,7 +504,7 @@ jobs:
         run: |
           PACKAGE_VERSION="$(python -c "${VER_SCRIPT1} ${VER_SCRIPT2}")"
           export PACKAGE_VERSION
-          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" pytest "python=${{ matrix.python }}" "numpy=${{ matrix.numpy }}" -c "$GITHUB_WORKSPACE/channel" -c conda-forge --override-channels
+          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" pytest "python=${{ matrix.python_spec || matrix.python }}" "numpy=${{ matrix.numpy }}" -c "$GITHUB_WORKSPACE/channel" -c conda-forge --override-channels
           # Test installed packages
           conda list -n "${{ env.TEST_ENV_NAME }}"
 
diff --git a/.github/workflows/conda-package.yml b/.github/workflows/conda-package.yml
index a9333f6c..cef4892a 100644
--- a/.github/workflows/conda-package.yml
+++ b/.github/workflows/conda-package.yml
@@ -32,6 +32,12 @@ jobs:
             numpy: "2.3"
           - python: "3.14"
             numpy: "2.3"
+            python_spec: "3.14.* *_cp314"
+            python_tag: "3.14"
+          - python: "3.14"
+            numpy: "2.3"
+            python_spec: "3.14.* *_cp314t"
+            python_tag: "3.14t"
 
     steps:
       - name: Cancel Previous Runs
@@ -54,9 +60,9 @@ jobs:
         with:
           path: ~/.conda/pkgs
           key:
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-${{hashFiles('**/meta.yaml') }}
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-${{hashFiles('**/meta.yaml') }}
           restore-keys: |
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-
             ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-
 
       - name: Add conda to system path
@@ -85,7 +91,7 @@ jobs:
       - name: Build conda package
         run: |
           CHANNELS=(-c https://software.repos.intel.com/python/conda -c conda-forge --override-channels)
-          VERSIONS=(--python "${{ matrix.python }}" --numpy "${{ matrix.numpy }}")
+          VERSIONS=(--python "${{ matrix.python_spec || matrix.python }}" --numpy "${{ matrix.numpy }}")
           TEST=(--no-test)
 
           conda build \
@@ -97,7 +103,7 @@ jobs:
       - name: Upload artifact
         uses: actions/upload-artifact@043fb46d1a93c77aae656e7c1c64a875d1fc6a0a # v7.0.1
         with:
-          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python }}
+          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python_tag || matrix.python }}
           path: ${{ env.CONDA_BLD }}${{ env.PACKAGE_NAME }}-*.conda
 
   build_windows:
@@ -118,6 +124,12 @@ jobs:
             numpy: "2.3"
           - python: "3.14"
             numpy: "2.3"
+            python_spec: "3.14.* *_cp314"
+            python_tag: "3.14"
+          - python: "3.14"
+            numpy: "2.3"
+            python_spec: "3.14.* *_cp314t"
+            python_tag: "3.14t"
 
     steps:
       - name: Cancel Previous Runs
@@ -145,9 +157,9 @@ jobs:
         with:
           path: /home/runner/conda_pkgs_dir
           key:
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-${{hashFiles('**/meta.yaml') }}
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-${{hashFiles('**/meta.yaml') }}
           restore-keys: |
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-
             ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-
 
       - name: Store conda paths as envs
@@ -164,12 +176,12 @@ jobs:
           conda list -n base
 
       - name: Build conda package
-        run: conda build --no-test --python "${{ matrix.python }}" --numpy "${{ matrix.numpy }}" -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels conda-recipe
+        run: conda build --no-test --python "${{ matrix.python_spec || matrix.python }}" --numpy "${{ matrix.numpy }}" -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels conda-recipe
 
       - name: Upload artifact
         uses: actions/upload-artifact@043fb46d1a93c77aae656e7c1c64a875d1fc6a0a # v7.0.1
         with:
-          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python }}
+          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python_tag || matrix.python }}
           path: ${{ env.CONDA_BLD }}${{ env.PACKAGE_NAME }}-*.conda
 
   test_linux:
@@ -178,9 +190,20 @@ jobs:
 
     strategy:
       matrix:
-        python: ["3.10", "3.11", "3.12", "3.13", "3.14"]
+        python: ["3.10", "3.11", "3.12", "3.13"]
         numpy: ['numpy">=2"']
         experimental: [false]
+        include:
+          - python: "3.14"
+            python_spec: "3.14.* *_cp314"
+            python_tag: "3.14"
+            numpy: 'numpy">=2"'
+            experimental: false
+          - python: "3.14"
+            python_spec: "3.14.* *_cp314t"
+            python_tag: "3.14t"
+            numpy: 'numpy">=2"'
+            experimental: false
 
     env:
       CHANNELS: -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels
@@ -189,7 +212,7 @@ jobs:
       - name: Download artifact
         uses: actions/download-artifact@3e5f45b2cfb9172054b4087a40e8e0b5a5461e7c # v8.0.1
         with:
-          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python }}
+          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python_tag || matrix.python }}
 
       - name: Add conda to system path
         run: echo "$CONDA/bin" >> "$GITHUB_PATH"
@@ -225,7 +248,7 @@ jobs:
           . "$CONDA/etc/profile.d/conda.sh"
           PACKAGE_VERSION="$(python -c "${VER_SCRIPT1} ${VER_SCRIPT2}")"
           export PACKAGE_VERSION
-          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" "python=${{ matrix.python }}" ${{ matrix.numpy }} -c "$GITHUB_WORKSPACE/channel" -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels --only-deps --dry-run > lockfile
+          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" "python=${{ matrix.python_spec || matrix.python }}" ${{ matrix.numpy }} -c "$GITHUB_WORKSPACE/channel" -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels --only-deps --dry-run > lockfile
           cat lockfile
 
       - name: Set pkgs_dirs
@@ -239,9 +262,9 @@ jobs:
         with:
           path: ~/.conda/pkgs
           key:
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-${{hashFiles('lockfile') }}
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-${{hashFiles('lockfile') }}
           restore-keys: |
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-
             ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-
 
       - name: Install mkl_random
@@ -249,7 +272,7 @@ jobs:
           . "$CONDA/etc/profile.d/conda.sh"
           PACKAGE_VERSION="$(python -c "${VER_SCRIPT1} ${VER_SCRIPT2}")"
           export PACKAGE_VERSION
-          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" pytest "python=${{ matrix.python }}" ${{ matrix.numpy }} -c "$GITHUB_WORKSPACE/channel" -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels
+          conda create -n "${{ env.TEST_ENV_NAME }}" "${PACKAGE_NAME}=${PACKAGE_VERSION}" pytest "python=${{ matrix.python_spec || matrix.python }}" ${{ matrix.numpy }} -c "$GITHUB_WORKSPACE/channel" -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels
           conda activate "${{ env.TEST_ENV_NAME }}"
 
           # Test installed packages
@@ -267,9 +290,20 @@ jobs:
 
     strategy:
       matrix:
-        python: ["3.10", "3.11", "3.12", "3.13", "3.14"]
+        python: ["3.10", "3.11", "3.12", "3.13"]
         numpy: ['numpy">=2"']
         experimental: [false]
+        include:
+          - python: "3.14"
+            python_spec: "3.14.* *_cp314"
+            python_tag: "3.14"
+            numpy: 'numpy">=2"'
+            experimental: false
+          - python: "3.14"
+            python_spec: "3.14.* *_cp314t"
+            python_tag: "3.14t"
+            numpy: 'numpy">=2"'
+            experimental: false
 
     env:
       CHANNELS: -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels
@@ -278,7 +312,7 @@ jobs:
       - name: Download artifact
         uses: actions/download-artifact@3e5f45b2cfb9172054b4087a40e8e0b5a5461e7c # v8.0.1
         with:
-          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python }}
+          name: ${{ env.PACKAGE_NAME }} ${{ runner.os }} Python ${{ matrix.python_tag || matrix.python }}
 
       - uses: conda-incubator/setup-miniconda@8ee1f361103df19b6f8c8655fd3967a8ecb162d5 # v4.0.1
         with:
@@ -330,7 +364,7 @@ jobs:
           FOR /F "tokens=* USEBACKQ" %%F IN (`python -c "%SCRIPT%"`) DO (
              SET PACKAGE_VERSION=%%F
           )
-          conda install -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% python=${{ matrix.python }} ${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }} --only-deps --dry-run > lockfile
+          conda install -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% "python=${{ matrix.python_spec || matrix.python }}" ${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }} --only-deps --dry-run > lockfile
 
       - name: Cache conda packages
         uses: actions/cache@55cc8345863c7cc4c66a329aec7e433d2d1c52a9 # v6.1.0
@@ -339,9 +373,9 @@ jobs:
         with:
           path: /home/runner/conda_pkgs_dir
           key:
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-${{hashFiles('lockfile') }}
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-${{hashFiles('lockfile') }}
           restore-keys: |
-            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python }}-
+            ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-python-${{ matrix.python_tag || matrix.python }}-
             ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-
 
         # add intel-openmp as an explicit dependency
@@ -356,7 +390,7 @@ jobs:
              SET PACKAGE_VERSION=%%F
           )
           SET "WORKAROUND_DEPENDENCIES=intel-openmp"
-          conda create -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% %WORKAROUND_DEPENDENCIES% pytest python=${{ matrix.python }} ${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }}
+          conda create -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% %WORKAROUND_DEPENDENCIES% pytest "python=${{ matrix.python_spec || matrix.python }}" ${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }}
           conda activate ${{ env.TEST_ENV_NAME }}
 
           # Test installed packages

From 9b41e87a9e0c5f4a7f442be952a3b10a7801b4a9 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 12 Aug 2026 04:20:41 -0700
Subject: [PATCH 04/30] Add free-threading test for mkl_random

---
 mkl_random/tests/test_freethreading.py | 129 +++++++++++++++++++++++++
 1 file changed, 129 insertions(+)
 create mode 100644 mkl_random/tests/test_freethreading.py

diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
new file mode 100644
index 00000000..d5c4a5ba
--- /dev/null
+++ b/mkl_random/tests/test_freethreading.py
@@ -0,0 +1,129 @@
+# Copyright (c) 2026, Intel Corporation
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+#     * Redistributions of source code must retain the above copyright notice,
+#       this list of conditions and the following disclaimer.
+#     * Redistributions in binary form must reproduce the above copyright
+#       notice, this list of conditions and the following disclaimer in the
+#       documentation and/or other materials provided with the distribution.
+#     * Neither the name of Intel Corporation nor the names of its contributors
+#       may be used to endorse or promote products derived from this software
+#       without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+import os
+import sys
+import sysconfig
+import threading
+
+import numpy as np
+import pytest
+
+# Oversubscription: MKL spawns its own thread pool per calling thread, so
+# generating from many Python threads concurrently can spawn far more OS
+# threads than cores. Cap it before mkl_random/MKL initialize.
+os.environ.setdefault("MKL_NUM_THREADS", "1")
+
+import mkl_random  # noqa: E402
+
+FREE_THREADED = bool(sysconfig.get_config_var("Py_GIL_DISABLED"))
+
+
+def test_concurrent_sampling_per_instance():
+    # Each thread owns a private MKLRandomState seeded identically, so the
+    # per-instance lock + `nogil` sampling must reproduce the single-threaded
+    # result exactly regardless of concurrency.
+    n_threads = 4
+    size = 10**5 + 1  # large enough that per-thread nogil sampling overlaps
+    seed = 1234
+
+    expected = mkl_random.MKLRandomState(seed).normal(size=size)
+
+    results = [None] * n_threads
+    errors = []
+
+    def worker(i):
+        try:
+            rs = mkl_random.MKLRandomState(seed)
+            results[i] = rs.normal(size=size)
+        except Exception as exc:  # pylint: disable=broad-except
+            errors.append(exc)
+
+    threads = [
+        threading.Thread(target=worker, args=(i,)) for i in range(n_threads)
+    ]
+    for t in threads:
+        t.start()
+    for t in threads:
+        t.join()
+
+    assert not errors
+
+    for i in range(n_threads):
+        np.testing.assert_array_equal(results[i], expected)
+
+
+def test_concurrent_shared_singleton():
+    # Module-level functions share a single lock-guarded RandomState. Hammering
+    # it from many threads must not corrupt state, crash, or return garbage.
+    n_threads = 8
+    size = 10**5 + 1
+    results = [None] * n_threads
+    errors = []
+
+    def worker(i):
+        try:
+            results[i] = mkl_random.uniform(size=size)
+        except Exception as exc:  # pylint: disable=broad-except
+            errors.append(exc)
+
+    threads = [
+        threading.Thread(target=worker, args=(i,)) for i in range(n_threads)
+    ]
+    for t in threads:
+        t.start()
+    for t in threads:
+        t.join()
+
+    assert not errors
+
+    for i in range(n_threads):
+        assert results[i].shape == (size,)
+        assert np.all(np.isfinite(results[i]))
+
+
+def test_concurrent_patch_restore():
+    n_threads = 8
+    n_iters = 20
+
+    def worker():
+        for _ in range(n_iters):
+            mkl_random.patch_numpy_random()
+            mkl_random.restore_numpy_random()
+
+    threads = [threading.Thread(target=worker) for _ in range(n_threads)]
+    for t in threads:
+        t.start()
+    for t in threads:
+        t.join()
+
+    assert not mkl_random.is_patched()
+
+
+@pytest.mark.skipif(
+    not FREE_THREADED, reason="requires a free-threaded CPython build"
+)
+def test_gil_not_reenabled_on_import():
+    assert not sys._is_gil_enabled()  # pylint: disable=no-member

From 4de50ebc1e12fbb77c73ab0d9289569a7141c7b4 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 12 Aug 2026 04:21:46 -0700
Subject: [PATCH 05/30] Update README build deps to cython>=3.1.0

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 7e5f06a8..ce6d04ba 100644
--- a/README.md
+++ b/README.md
@@ -151,7 +151,7 @@ If you already have `mkl` and `numpy` installed (from your system or a conda env
 and want to reuse them instead of pulling fresh copies into an isolated build, first
 install the build dependencies:
 ```sh
-pip install meson-python cmake ninja cython numpy mkl-devel
+pip install meson-python cmake ninja "cython>=3.1.0" numpy mkl-devel
 ```
 
 then build against the existing installation with:

From 8e5a8e5885b2bb3e514c927fa3f9a084c590ad9e Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 12 Aug 2026 04:57:37 -0700
Subject: [PATCH 06/30] Add a missing continue-on-error to conda-package.yml

---
 .github/workflows/conda-package.yml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.github/workflows/conda-package.yml b/.github/workflows/conda-package.yml
index cef4892a..2334fb27 100644
--- a/.github/workflows/conda-package.yml
+++ b/.github/workflows/conda-package.yml
@@ -204,6 +204,7 @@ jobs:
             python_tag: "3.14t"
             numpy: 'numpy">=2"'
             experimental: false
+    continue-on-error: ${{ matrix.experimental }}
 
     env:
       CHANNELS: -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels
@@ -304,6 +305,7 @@ jobs:
             python_tag: "3.14t"
             numpy: 'numpy">=2"'
             experimental: false
+    continue-on-error: ${{ matrix.experimental }}
 
     env:
       CHANNELS: -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels

From 4b9cfe7111f36eb11de4887ff4ad1f404802504f Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 12 Aug 2026 05:04:59 -0700
Subject: [PATCH 07/30] Add gh-159 to CHANGELOG.md

---
 CHANGELOG.md | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index eff67b71..52eb980b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,8 +7,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 # [dev] (MM/DD/YYYY)
 
 ### Added
+* Added support for free-threaded (GIL-disabled) CPython builds: the Cython extension is compiled with `freethreading_compatible=True`, so importing `mkl_random` no longer re-enables the GIL [gh-159](https://github.com/IntelPython/mkl_random/pull/159)
 
 ### Changed
+* Raised the minimum build-time `Cython` requirement to `3.1.0`, the first release providing the `freethreading_compatible` directive [gh-159](https://github.com/IntelPython/mkl_random/pull/159)
+
+### Removed
+* Removed the `python-gil` constraint from the conda recipes, which pinned `mkl_random` to GIL-enabled Python 3.14 builds [gh-159](https://github.com/IntelPython/mkl_random/pull/159)
 
 ### Fixed
 

From be5c348130ec8ac5784ef42d9e3f3ae082a9f73f Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 12 Aug 2026 13:51:52 -0700
Subject: [PATCH 08/30] Create fresh Windows test env for cp314t

---
 .github/workflows/conda-package-cf.yml | 5 ++---
 .github/workflows/conda-package.yml    | 6 ++----
 2 files changed, 4 insertions(+), 7 deletions(-)

diff --git a/.github/workflows/conda-package-cf.yml b/.github/workflows/conda-package-cf.yml
index 516e1b09..4590abcc 100644
--- a/.github/workflows/conda-package-cf.yml
+++ b/.github/workflows/conda-package-cf.yml
@@ -290,8 +290,7 @@ jobs:
         with:
           miniforge-version: latest
           channels: conda-forge
-          activate-environment: ${{ env.TEST_ENV_NAME }}
-          python-version: ${{ matrix.python }}
+          activate-environment: base
 
       - name: Install conda-index
         run: |
@@ -316,7 +315,7 @@ jobs:
           FOR /F "tokens=* USEBACKQ" %%F IN (`python -c "%SCRIPT%"`) DO (
              SET PACKAGE_VERSION=%%F
           )
-          conda install -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% "python=${{ matrix.python_spec || matrix.python }}" numpy=${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }} --only-deps --dry-run > lockfile
+          conda create -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% "python=${{ matrix.python_spec || matrix.python }}" numpy=${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }} --only-deps --dry-run > lockfile
 
       - name: Cache conda packages
         uses: actions/cache@55cc8345863c7cc4c66a329aec7e433d2d1c52a9 # v6.1.0
diff --git a/.github/workflows/conda-package.yml b/.github/workflows/conda-package.yml
index 2334fb27..946324ba 100644
--- a/.github/workflows/conda-package.yml
+++ b/.github/workflows/conda-package.yml
@@ -318,11 +318,9 @@ jobs:
 
       - uses: conda-incubator/setup-miniconda@8ee1f361103df19b6f8c8655fd3967a8ecb162d5 # v4.0.1
         with:
-          auto-update-conda: true
           miniforge-version: latest
           channels: conda-forge
-          activate-environment: ${{ env.TEST_ENV_NAME }}
-          python-version: ${{ matrix.python }}
+          activate-environment: base
 
       - name: Install conda-index
         run: |
@@ -366,7 +364,7 @@ jobs:
           FOR /F "tokens=* USEBACKQ" %%F IN (`python -c "%SCRIPT%"`) DO (
              SET PACKAGE_VERSION=%%F
           )
-          conda install -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% "python=${{ matrix.python_spec || matrix.python }}" ${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }} --only-deps --dry-run > lockfile
+          conda create -n ${{ env.TEST_ENV_NAME }} ${{ env.PACKAGE_NAME }}=%PACKAGE_VERSION% "python=${{ matrix.python_spec || matrix.python }}" ${{ matrix.numpy }} -c ${{ env.GITHUB_WORKSPACE }}/channel ${{ env.CHANNELS }} --only-deps --dry-run > lockfile
 
       - name: Cache conda packages
         uses: actions/cache@55cc8345863c7cc4c66a329aec7e433d2d1c52a9 # v6.1.0

From 33e696b0365dbea7a050e8ea31567d72e3193837 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Mon, 7 Sep 2026 04:09:44 -0700
Subject: [PATCH 09/30] Fix data races in unguarded RNG methods under
 free-threading

---
 mkl_random/mklrand.pyx | 92 +++++++++++++++++++++++-------------------
 1 file changed, 50 insertions(+), 42 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index f0a2a5e9..c4491c1a 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -1632,13 +1632,13 @@ cdef class _MKLRandomState:
         cdef unsigned int stream_id
         cdef cnp.ndarray obj "arrayObject_obj"
 
-        if (brng):
-            brng_token, stream_id = _parse_brng_argument(brng)
-        else:
-            brng_token = <irk_brng_t> irk_get_brng_and_stream_mkl(
-                self.internal_state, &stream_id
-            )
         with self.lock:
+            if (brng):
+                brng_token, stream_id = _parse_brng_argument(brng)
+            else:
+                brng_token = <irk_brng_t> irk_get_brng_and_stream_mkl(
+                    self.internal_state, &stream_id
+                )
             try:
                 if seed is None:
                     _errcode = irk_randomseed_mkl(
@@ -6599,7 +6599,8 @@ cdef class _MKLRandomState:
             raise ValueError("n < 0")
         # numpy#20483: Avoids divide by 0
         niter = sz // d if d else 0
-        irk_multinomial_vec(self.internal_state, niter, mnix, n, d, pix)
+        with self.lock, nogil:
+            irk_multinomial_vec(self.internal_state, niter, mnix, n, d, pix)
 
         return multin
 
@@ -6976,7 +6977,8 @@ cdef class MKLRandomState(_MKLRandomState):
         """
         cdef int err, brng_id
 
-        err = irk_leapfrog_stream_mkl(self.internal_state, k, nstreams)
+        with self.lock:
+            err = irk_leapfrog_stream_mkl(self.internal_state, k, nstreams)
 
         if err == -1:
             raise ValueError("The stream state buffer is corrupted")
@@ -6997,7 +6999,8 @@ cdef class MKLRandomState(_MKLRandomState):
         """
         cdef int err, brng_id
 
-        err = irk_skipahead_stream_mkl(self.internal_state, nskips)
+        with self.lock:
+            err = irk_skipahead_stream_mkl(self.internal_state, nskips)
 
         if err == -1:
             raise ValueError("The stream state buffer is corrupted")
@@ -7080,9 +7083,10 @@ cdef class MKLRandomState(_MKLRandomState):
 
         if ((<int> lo) == lo) and ((<int>hi) == hi):
             if size is None:
-                irk_discrete_uniform_vec(
-                    self.internal_state, 1, &rv_int, <int>lo, <int>hi
-                )
+                with self.lock, nogil:
+                    irk_discrete_uniform_vec(
+                        self.internal_state, 1, &rv_int, <int>lo, <int>hi
+                    )
                 return rv_int
             else:
                 array = <cnp.ndarray>np.empty(size, np.int32)
@@ -7099,9 +7103,10 @@ cdef class MKLRandomState(_MKLRandomState):
                 return array
         else:
             if size is None:
-                irk_discrete_uniform_long_vec(
-                    self.internal_state, 1, &rv_long, lo, hi
-                )
+                with self.lock, nogil:
+                    irk_discrete_uniform_long_vec(
+                        self.internal_state, 1, &rv_long, lo, hi
+                    )
                 return rv_long
             else:
                 array = <cnp.ndarray>np.empty(size, int)
@@ -7284,35 +7289,38 @@ cdef class MKLRandomState(_MKLRandomState):
             method, [ICDF, BOXMULLER2, BOXMULLER], _method_alias_dict_gaussian
         )
         if (method is ICDF):
-            irk_multinormal_vec_ICDF(
-                self.internal_state,
-                n,
-                res_data,
-                dim,
-                mean_data,
-                t_data,
-                storage_mode
-            )
+            with self.lock, nogil:
+                irk_multinormal_vec_ICDF(
+                    self.internal_state,
+                    n,
+                    res_data,
+                    dim,
+                    mean_data,
+                    t_data,
+                    storage_mode
+                )
         elif (method is BOXMULLER2):
-            irk_multinormal_vec_BM2(
-                self.internal_state,
-                n,
-                res_data,
-                dim,
-                mean_data,
-                t_data,
-                storage_mode
-            )
+            with self.lock, nogil:
+                irk_multinormal_vec_BM2(
+                    self.internal_state,
+                    n,
+                    res_data,
+                    dim,
+                    mean_data,
+                    t_data,
+                    storage_mode
+                )
         else:
-            irk_multinormal_vec_BM1(
-                self.internal_state,
-                n,
-                res_data,
-                dim,
-                mean_data,
-                t_data,
-                storage_mode
-            )
+            with self.lock, nogil:
+                irk_multinormal_vec_BM1(
+                    self.internal_state,
+                    n,
+                    res_data,
+                    dim,
+                    mean_data,
+                    t_data,
+                    storage_mode
+                )
 
         return resarr
 

From 5d864420119b0606e33e664e647fcb191627fdfe Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Mon, 7 Sep 2026 04:30:18 -0700
Subject: [PATCH 10/30] Add shared-stream race regression test for
 free-threading

---
 mkl_random/tests/test_freethreading.py | 93 +++++++++++++++-----------
 1 file changed, 53 insertions(+), 40 deletions(-)

diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index d5c4a5ba..a2aaa3d1 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -41,28 +41,18 @@
 FREE_THREADED = bool(sysconfig.get_config_var("Py_GIL_DISABLED"))
 
 
-def test_concurrent_sampling_per_instance():
-    # Each thread owns a private MKLRandomState seeded identically, so the
-    # per-instance lock + `nogil` sampling must reproduce the single-threaded
-    # result exactly regardless of concurrency.
-    n_threads = 4
-    size = 10**5 + 1  # large enough that per-thread nogil sampling overlaps
-    seed = 1234
-
-    expected = mkl_random.MKLRandomState(seed).normal(size=size)
-
-    results = [None] * n_threads
+def _run_on_threads(worker, n_threads):
+    # Run worker(i) on n_threads and fail if any thread raised.
     errors = []
 
-    def worker(i):
+    def wrapped(i):
         try:
-            rs = mkl_random.MKLRandomState(seed)
-            results[i] = rs.normal(size=size)
+            worker(i)
         except Exception as exc:  # pylint: disable=broad-except
             errors.append(exc)
 
     threads = [
-        threading.Thread(target=worker, args=(i,)) for i in range(n_threads)
+        threading.Thread(target=wrapped, args=(i,)) for i in range(n_threads)
     ]
     for t in threads:
         t.start()
@@ -71,8 +61,24 @@ def worker(i):
 
     assert not errors
 
-    for i in range(n_threads):
-        np.testing.assert_array_equal(results[i], expected)
+
+def test_concurrent_sampling_per_instance():
+    # Each thread owns a private MKLRandomState seeded identically, so the
+    # per-instance lock + `nogil` sampling must reproduce the single-threaded
+    # result exactly regardless of concurrency.
+    n_threads = 4
+    size = 10**5 + 1  # large enough that per-thread nogil sampling overlaps
+    seed = 1234
+    expected = mkl_random.MKLRandomState(seed).normal(size=size)
+    results = [None] * n_threads
+
+    def worker(i):
+        results[i] = mkl_random.MKLRandomState(seed).normal(size=size)
+
+    _run_on_threads(worker, n_threads)
+
+    for r in results:
+        np.testing.assert_array_equal(r, expected)
 
 
 def test_concurrent_shared_singleton():
@@ -81,47 +87,54 @@ def test_concurrent_shared_singleton():
     n_threads = 8
     size = 10**5 + 1
     results = [None] * n_threads
-    errors = []
 
     def worker(i):
-        try:
-            results[i] = mkl_random.uniform(size=size)
-        except Exception as exc:  # pylint: disable=broad-except
-            errors.append(exc)
+        results[i] = mkl_random.uniform(size=size)
 
-    threads = [
-        threading.Thread(target=worker, args=(i,)) for i in range(n_threads)
-    ]
-    for t in threads:
-        t.start()
-    for t in threads:
-        t.join()
+    _run_on_threads(worker, n_threads)
 
-    assert not errors
-
-    for i in range(n_threads):
-        assert results[i].shape == (size,)
-        assert np.all(np.isfinite(results[i]))
+    for r in results:
+        assert r.shape == (size,)
+        assert np.all(np.isfinite(r))
 
 
 def test_concurrent_patch_restore():
     n_threads = 8
     n_iters = 20
 
-    def worker():
+    def worker(_i):
         for _ in range(n_iters):
             mkl_random.patch_numpy_random()
             mkl_random.restore_numpy_random()
 
-    threads = [threading.Thread(target=worker) for _ in range(n_threads)]
-    for t in threads:
-        t.start()
-    for t in threads:
-        t.join()
+    _run_on_threads(worker, n_threads)
 
     assert not mkl_random.is_patched()
 
 
+def test_concurrent_multinormal_cholesky_shared():
+    # A data race on the shared stream reuses values; assert few duplicates.
+    n_threads = 8
+    mean = np.zeros(1)
+    ch = np.eye(1)
+    rng = mkl_random.MKLRandomState(12345)
+    chunks = [None] * n_threads
+
+    def worker(i):
+        parts = [
+            rng.multinormal_cholesky(mean, ch, size=4000).ravel()
+            for _ in range(25)
+        ]
+        chunks[i] = np.concatenate(parts)
+
+    _run_on_threads(worker, n_threads)
+
+    allvals = np.concatenate(chunks)
+    assert np.all(np.isfinite(allvals))
+    dup_frac = 1.0 - np.unique(allvals).size / allvals.size
+    assert dup_frac < 0.01, f"shared stream corrupted: {dup_frac:.3%} dups"
+
+
 @pytest.mark.skipif(
     not FREE_THREADED, reason="requires a free-threaded CPython build"
 )

From a6d8a29fa7cc6b7ddce3b670ec89f5408fcb44f8 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 9 Sep 2026 03:38:50 -0700
Subject: [PATCH 11/30] Fix re-entrancy deadlocks in shuffle and _seed_impl

---
 mkl_random/mklrand.pyx                 | 54 +++++++++++++-------------
 mkl_random/tests/test_freethreading.py | 19 +++++++++
 2 files changed, 47 insertions(+), 26 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index c4491c1a..de1460b0 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -1632,10 +1632,13 @@ cdef class _MKLRandomState:
         cdef unsigned int stream_id
         cdef cnp.ndarray obj "arrayObject_obj"
 
+        if (brng):
+            # Parse before the lock to avoid warn
+            brng_token, stream_id = _parse_brng_argument(brng)
+
         with self.lock:
-            if (brng):
-                brng_token, stream_id = _parse_brng_argument(brng)
-            else:
+            if not brng:
+                # Reads state->stream, which a concurrent seed can free.
                 brng_token = <irk_brng_t> irk_get_brng_and_stream_mkl(
                     self.internal_state, &stream_id
                 )
@@ -6775,6 +6778,8 @@ cdef class _MKLRandomState:
 
         u = <cnp.ndarray>self.random_sample(n - 1)
         u_data = <double*>cnp.PyArray_DATA(u)
+        # Indices are already drawn under the lock; the swaps touch no stream
+        # state and run unlocked (locking across the callback would deadlock).
 
         if type(x) is np.ndarray and x.ndim == 1 and x.size:
             # Fast, statically typed path: shuffle the underlying buffer.
@@ -6789,18 +6794,17 @@ cdef class _MKLRandomState:
             # when the function exits.
             buf = np.empty(itemsize, dtype=np.int8)  # GC'd at function exit
             buf_ptr = cnp.PyArray_BYTES(buf)
-            with self.lock:
-                # We trick gcc into providing a specialized implementation for
-                # the most common case, yielding a ~33% performance improvement.
-                # Note that apparently, only one branch can ever be specialized.
-                if itemsize == sizeof(cnp.npy_intp):
-                    self._shuffle_raw(
-                        n, sizeof(cnp.npy_intp), stride, x_ptr, buf_ptr, u_data
-                    )
-                else:
-                    self._shuffle_raw(
-                        n, itemsize, stride, x_ptr, buf_ptr, u_data
-                    )
+            # We trick gcc into providing a specialized implementation for
+            # the most common case, yielding a ~33% performance improvement.
+            # Note that apparently, only one branch can ever be specialized.
+            if itemsize == sizeof(cnp.npy_intp):
+                self._shuffle_raw(
+                    n, sizeof(cnp.npy_intp), stride, x_ptr, buf_ptr, u_data
+                )
+            else:
+                self._shuffle_raw(
+                    n, itemsize, stride, x_ptr, buf_ptr, u_data
+                )
         elif isinstance(x, np.ndarray):
             if x.size == 0:
                 # shuffling is a no-op
@@ -6814,13 +6818,12 @@ cdef class _MKLRandomState:
                         UserWarning, stacklevel=1)  # Cython adds no stacklevel
 
             buf = np.empty_like(x[0, ...])
-            with self.lock:
-                for i in reversed(range(1, n)):
-                    j = <cnp.npy_intp>floor((i + 1) * u_data[i - 1])
-                    if (j < i):
-                        buf[...] = x[j]
-                        x[j] = x[i]
-                        x[i] = buf
+            for i in reversed(range(1, n)):
+                j = <cnp.npy_intp>floor((i + 1) * u_data[i - 1])
+                if (j < i):
+                    buf[...] = x[j]
+                    x[j] = x[i]
+                    x[i] = buf
         else:
             # Untyped path.
             if not isinstance(x, Sequence):
@@ -6831,10 +6834,9 @@ cdef class _MKLRandomState:
                     "E.g., non-numpy array/tensor objects with view semantics "
                     "may contain duplicates after shuffling.",
                     UserWarning, stacklevel=1)  # Cython does not add a level
-            with self.lock:
-                for i in reversed(range(1, n)):
-                    j = <cnp.npy_intp>floor((i + 1) * u_data[i - 1])
-                    x[i], x[j] = x[j], x[i]
+            for i in reversed(range(1, n)):
+                j = <cnp.npy_intp>floor((i + 1) * u_data[i - 1])
+                x[i], x[j] = x[j], x[i]
 
     cdef inline _shuffle_raw(
         self,
diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index a2aaa3d1..d8dfd937 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -135,6 +135,25 @@ def worker(i):
     assert dup_frac < 0.01, f"shared stream corrupted: {dup_frac:.3%} dups"
 
 
+def test_shuffle_reentrancy():
+    # shuffle must not hold the lock across a user callback.
+    rs = mkl_random.MKLRandomState(1)
+
+    class ReentrantList(list):
+        def __setitem__(self, i, v):
+            rs.uniform(size=1)
+            super().__setitem__(i, v)
+
+    done = threading.Event()
+
+    def run():
+        rs.shuffle(ReentrantList(range(8)))
+        done.set()
+
+    threading.Thread(target=run, daemon=True).start()
+    assert done.wait(timeout=30), "shuffle deadlocked on re-entrant callback"
+
+
 @pytest.mark.skipif(
     not FREE_THREADED, reason="requires a free-threaded CPython build"
 )

From 09e4cff0f2a3c22dfd074ebb79b6666632526629 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 9 Sep 2026 03:59:42 -0700
Subject: [PATCH 12/30] Lock size=None sampling paths and cover them with a
 multiset test

---
 mkl_random/mklrand.pyx                 | 99 +++++++++++++++-----------
 mkl_random/tests/test_freethreading.py | 60 +++++++++++-----
 2 files changed, 99 insertions(+), 60 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index de1460b0..2a0fcfac 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -563,7 +563,8 @@ cdef object vec_cont0_array(
     cdef cnp.npy_intp length
 
     if size is None:
-        func(state, 1, &res)
+        with lock, nogil:
+            func(state, 1, &res)
         return res
     else:
         array = <cnp.ndarray>np.empty(size, np.float64)
@@ -583,7 +584,8 @@ cdef object vec_cont1_array_sc(
     cdef cnp.npy_intp length
 
     if size is None:
-        func(state, 1, &res, a)
+        with lock, nogil:
+            func(state, 1, &res, a)
         return res
     else:
         array = <cnp.ndarray>np.empty(size, np.float64)
@@ -679,7 +681,8 @@ cdef object vec_cont2_array_sc(
     cdef cnp.npy_intp length
 
     if size is None:
-        func(state, 1, &res, a, b)
+        with lock, nogil:
+            func(state, 1, &res, a, b)
         return res
     else:
         array = <cnp.ndarray>np.empty(size, np.float64)
@@ -778,7 +781,8 @@ cdef object vec_cont3_array_sc(
     cdef cnp.npy_intp length
 
     if size is None:
-        func(state, 1, &res, a, b, c)
+        with lock, nogil:
+            func(state, 1, &res, a, b, c)
         return res
     else:
         array = <cnp.ndarray>np.empty(size, np.float64)
@@ -885,7 +889,8 @@ cdef object vec_long_disc0_array(
     cdef cnp.npy_intp length
 
     if size is None:
-        func(state, 1, &res)
+        with lock, nogil:
+            func(state, 1, &res)
         return res
     array = <cnp.ndarray>np.empty(size, np.dtype("long"))
     length = cnp.PyArray_SIZE(array)
@@ -910,7 +915,8 @@ cdef object vec_discnp_array_sc(
     cdef cnp.npy_intp length
 
     if size is None:
-        func(state, 1, &res, n, p)
+        with lock, nogil:
+            func(state, 1, &res, n, p)
         return res
     else:
         array = <cnp.ndarray>np.empty(size, np.intc)
@@ -1009,7 +1015,8 @@ cdef object vec_discdd_array_sc(
     cdef cnp.npy_intp length
 
     if size is None:
-        func(state, 1, &res, n, p)
+        with lock, nogil:
+            func(state, 1, &res, n, p)
         return res
     else:
         array = <cnp.ndarray>np.empty(size, np.intc)
@@ -1109,7 +1116,8 @@ cdef object vec_discnmN_array_sc(
     cdef cnp.npy_intp length
 
     if size is None:
-        func(state, 1, &res, n, m, N)
+        with lock, nogil:
+            func(state, 1, &res, n, m, N)
         return res
     else:
         array = <cnp.ndarray>np.empty(size, np.intc)
@@ -1220,7 +1228,8 @@ cdef object vec_discd_array_sc(
     cdef cnp.npy_intp length
 
     if size is None:
-        func(state, 1, &res, a)
+        with lock, nogil:
+            func(state, 1, &res, a)
         return res
     else:
         array = <cnp.ndarray>np.empty(size, np.intc)
@@ -1244,7 +1253,8 @@ cdef object vec_long_discd_array_sc(
     cdef cnp.npy_intp length
 
     if size is None:
-        func(state, 1, &res, a)
+        with lock, nogil:
+            func(state, 1, &res, a)
         return res
     else:
         array = <cnp.ndarray>np.empty(size, np.dtype("long"))
@@ -1976,13 +1986,14 @@ cdef class _MKLRandomState:
         cdef cnp.npy_intp cnt
 
         if size is None:
-            irk_rand_bool_vec(self.internal_state, 1, &buf, low, high)
+            with self.lock, nogil:
+                irk_rand_bool_vec(self.internal_state, 1, &buf, low, high)
             return np.bool_(buf)
         else:
             array = <cnp.ndarray>np.empty(size, np.bool_)
             cnt = cnp.PyArray_SIZE(array)
             out = <cnp.npy_bool *>cnp.PyArray_DATA(array)
-            with nogil:
+            with self.lock, nogil:
                 irk_rand_bool_vec(self.internal_state, cnt, out, low, high)
             return array
 
@@ -1999,13 +2010,14 @@ cdef class _MKLRandomState:
         cdef cnp.npy_intp cnt
 
         if size is None:
-            irk_rand_int8_vec(self.internal_state, 1, &buf, low, high)
+            with self.lock, nogil:
+                irk_rand_int8_vec(self.internal_state, 1, &buf, low, high)
             return np.int8(<cnp.npy_int8>buf)
         else:
             array = <cnp.ndarray>np.empty(size, np.int8)
             cnt = cnp.PyArray_SIZE(array)
             out = <cnp.npy_int8 *>cnp.PyArray_DATA(array)
-            with nogil:
+            with self.lock, nogil:
                 irk_rand_int8_vec(self.internal_state, cnt, out, low, high)
             return array
 
@@ -2022,13 +2034,14 @@ cdef class _MKLRandomState:
         cdef cnp.npy_intp cnt
 
         if size is None:
-            irk_rand_int16_vec(self.internal_state, 1, &buf, low, high)
+            with self.lock, nogil:
+                irk_rand_int16_vec(self.internal_state, 1, &buf, low, high)
             return np.int16(<cnp.npy_int16>buf)
         else:
             array = <cnp.ndarray>np.empty(size, np.int16)
             cnt = cnp.PyArray_SIZE(array)
             out = <cnp.npy_int16 *>cnp.PyArray_DATA(array)
-            with nogil:
+            with self.lock, nogil:
                 irk_rand_int16_vec(self.internal_state, cnt, out, low, high)
             return array
 
@@ -2066,13 +2079,14 @@ cdef class _MKLRandomState:
         cdef cnp.npy_intp cnt
 
         if size is None:
-            irk_rand_int32_vec(self.internal_state, 1, &buf, low, high)
+            with self.lock, nogil:
+                irk_rand_int32_vec(self.internal_state, 1, &buf, low, high)
             return np.int32(buf)
         else:
             array = <cnp.ndarray>np.empty(size, np.int32)
             cnt = cnp.PyArray_SIZE(array)
             out = <cnp.npy_int32 *>cnp.PyArray_DATA(array)
-            with nogil:
+            with self.lock, nogil:
                 irk_rand_int32_vec(self.internal_state, cnt, out, low, high)
             return array
 
@@ -2089,13 +2103,14 @@ cdef class _MKLRandomState:
         cdef cnp.npy_intp cnt
 
         if size is None:
-            irk_rand_int64_vec(self.internal_state, 1, &buf, low, high)
+            with self.lock, nogil:
+                irk_rand_int64_vec(self.internal_state, 1, &buf, low, high)
             return np.int64(buf)
         else:
             array = <cnp.ndarray>np.empty(size, np.int64)
             cnt = cnp.PyArray_SIZE(array)
             out = <cnp.npy_int64 *>cnp.PyArray_DATA(array)
-            with nogil:
+            with self.lock, nogil:
                 irk_rand_int64_vec(self.internal_state, cnt, out, low, high)
             return array
 
@@ -2112,13 +2127,14 @@ cdef class _MKLRandomState:
         cdef cnp.npy_intp cnt
 
         if size is None:
-            irk_rand_uint8_vec(self.internal_state, 1, &buf, low, high)
+            with self.lock, nogil:
+                irk_rand_uint8_vec(self.internal_state, 1, &buf, low, high)
             return np.uint8(buf)
         else:
             array = <cnp.ndarray>np.empty(size, np.uint8)
             cnt = cnp.PyArray_SIZE(array)
             out = <cnp.npy_uint8 *>cnp.PyArray_DATA(array)
-            with nogil:
+            with self.lock, nogil:
                 irk_rand_uint8_vec(self.internal_state, cnt, out, low, high)
             return array
 
@@ -2135,13 +2151,14 @@ cdef class _MKLRandomState:
         cdef cnp.npy_intp cnt
 
         if size is None:
-            irk_rand_uint16_vec(self.internal_state, 1, &buf, low, high)
+            with self.lock, nogil:
+                irk_rand_uint16_vec(self.internal_state, 1, &buf, low, high)
             return np.uint16(buf)
         else:
             array = <cnp.ndarray>np.empty(size, np.uint16)
             cnt = cnp.PyArray_SIZE(array)
             out = <cnp.npy_uint16 *>cnp.PyArray_DATA(array)
-            with nogil:
+            with self.lock, nogil:
                 irk_rand_uint16_vec(self.internal_state, cnt, out, low, high)
             return array
 
@@ -2158,13 +2175,14 @@ cdef class _MKLRandomState:
         cdef cnp.npy_intp cnt
 
         if size is None:
-            irk_rand_uint32_vec(self.internal_state, 1, &buf, low, high)
+            with self.lock, nogil:
+                irk_rand_uint32_vec(self.internal_state, 1, &buf, low, high)
             return np.uint32(buf)
         else:
             array = <cnp.ndarray>np.empty(size, np.uint32)
             cnt = cnp.PyArray_SIZE(array)
             out = <cnp.npy_uint32 *>cnp.PyArray_DATA(array)
-            with nogil:
+            with self.lock, nogil:
                 irk_rand_uint32_vec(self.internal_state, cnt, out, low, high)
             return array
 
@@ -2181,13 +2199,14 @@ cdef class _MKLRandomState:
         cdef cnp.npy_intp cnt
 
         if size is None:
-            irk_rand_uint64_vec(self.internal_state, 1, &buf, low, high)
+            with self.lock, nogil:
+                irk_rand_uint64_vec(self.internal_state, 1, &buf, low, high)
             return np.uint64(buf)
         else:
             array = <cnp.ndarray>np.empty(size, np.uint64)
             cnt = cnp.PyArray_SIZE(array)
             out = <cnp.npy_uint64 *>cnp.PyArray_DATA(array)
-            with nogil:
+            with self.lock, nogil:
                 irk_rand_uint64_vec(self.internal_state, cnt, out, low, high)
             return array
 
@@ -2198,7 +2217,7 @@ cdef class _MKLRandomState:
         cdef cnp.npy_bool *out_p = <cnp.npy_bool *>cnp.PyArray_DATA(out)
         cdef cnp.npy_bool *low_p = <cnp.npy_bool *>cnp.PyArray_DATA(low)
         cdef cnp.npy_bool *high_p = <cnp.npy_bool *>cnp.PyArray_DATA(high)
-        with nogil:
+        with self.lock, nogil:
             irk_rand_bool_broadcast(
                 self.internal_state, cnt, out_p, low_p, high_p
             )
@@ -2209,7 +2228,7 @@ cdef class _MKLRandomState:
         cdef cnp.npy_int8 *out_p = <cnp.npy_int8 *>cnp.PyArray_DATA(out)
         cdef cnp.npy_int8 *low_p = <cnp.npy_int8 *>cnp.PyArray_DATA(low)
         cdef cnp.npy_int8 *high_p = <cnp.npy_int8 *>cnp.PyArray_DATA(high)
-        with nogil:
+        with self.lock, nogil:
             irk_rand_int8_broadcast(
                 self.internal_state, cnt, out_p, low_p, high_p
             )
@@ -2220,7 +2239,7 @@ cdef class _MKLRandomState:
         cdef cnp.npy_int16 *out_p = <cnp.npy_int16 *>cnp.PyArray_DATA(out)
         cdef cnp.npy_int16 *low_p = <cnp.npy_int16 *>cnp.PyArray_DATA(low)
         cdef cnp.npy_int16 *high_p = <cnp.npy_int16 *>cnp.PyArray_DATA(high)
-        with nogil:
+        with self.lock, nogil:
             irk_rand_int16_broadcast(
                 self.internal_state, cnt, out_p, low_p, high_p
             )
@@ -2231,7 +2250,7 @@ cdef class _MKLRandomState:
         cdef cnp.npy_int32 *out_p = <cnp.npy_int32 *>cnp.PyArray_DATA(out)
         cdef cnp.npy_int32 *low_p = <cnp.npy_int32 *>cnp.PyArray_DATA(low)
         cdef cnp.npy_int32 *high_p = <cnp.npy_int32 *>cnp.PyArray_DATA(high)
-        with nogil:
+        with self.lock, nogil:
             irk_rand_int32_broadcast(
                 self.internal_state, cnt, out_p, low_p, high_p
             )
@@ -2242,7 +2261,7 @@ cdef class _MKLRandomState:
         cdef cnp.npy_int64 *out_p = <cnp.npy_int64 *>cnp.PyArray_DATA(out)
         cdef cnp.npy_int64 *low_p = <cnp.npy_int64 *>cnp.PyArray_DATA(low)
         cdef cnp.npy_int64 *high_p = <cnp.npy_int64 *>cnp.PyArray_DATA(high)
-        with nogil:
+        with self.lock, nogil:
             irk_rand_int64_broadcast(
                 self.internal_state, cnt, out_p, low_p, high_p
             )
@@ -2253,7 +2272,7 @@ cdef class _MKLRandomState:
         cdef cnp.npy_uint8 *out_p = <cnp.npy_uint8 *>cnp.PyArray_DATA(out)
         cdef cnp.npy_uint8 *low_p = <cnp.npy_uint8 *>cnp.PyArray_DATA(low)
         cdef cnp.npy_uint8 *high_p = <cnp.npy_uint8 *>cnp.PyArray_DATA(high)
-        with nogil:
+        with self.lock, nogil:
             irk_rand_uint8_broadcast(
                 self.internal_state, cnt, out_p, low_p, high_p
             )
@@ -2264,7 +2283,7 @@ cdef class _MKLRandomState:
         cdef cnp.npy_uint16 *out_p = <cnp.npy_uint16 *>cnp.PyArray_DATA(out)
         cdef cnp.npy_uint16 *low_p = <cnp.npy_uint16 *>cnp.PyArray_DATA(low)
         cdef cnp.npy_uint16 *high_p = <cnp.npy_uint16 *>cnp.PyArray_DATA(high)
-        with nogil:
+        with self.lock, nogil:
             irk_rand_uint16_broadcast(
                 self.internal_state, cnt, out_p, low_p, high_p
             )
@@ -2275,7 +2294,7 @@ cdef class _MKLRandomState:
         cdef cnp.npy_uint32 *out_p = <cnp.npy_uint32 *>cnp.PyArray_DATA(out)
         cdef cnp.npy_uint32 *low_p = <cnp.npy_uint32 *>cnp.PyArray_DATA(low)
         cdef cnp.npy_uint32 *high_p = <cnp.npy_uint32 *>cnp.PyArray_DATA(high)
-        with nogil:
+        with self.lock, nogil:
             irk_rand_uint32_broadcast(
                 self.internal_state, cnt, out_p, low_p, high_p
             )
@@ -2286,7 +2305,7 @@ cdef class _MKLRandomState:
         cdef cnp.npy_uint64 *out_p = <cnp.npy_uint64 *>cnp.PyArray_DATA(out)
         cdef cnp.npy_uint64 *low_p = <cnp.npy_uint64 *>cnp.PyArray_DATA(low)
         cdef cnp.npy_uint64 *high_p = <cnp.npy_uint64 *>cnp.PyArray_DATA(high)
-        with nogil:
+        with self.lock, nogil:
             irk_rand_uint64_broadcast(
                 self.internal_state, cnt, out_p, low_p, high_p
             )
@@ -2330,8 +2349,7 @@ cdef class _MKLRandomState:
         high_c = np.ascontiguousarray(high_incl, dtype=_dtype)
         out = np.empty(out_shape, dtype=_dtype)
 
-        with self.lock:
-            broadcast_func(low_c, high_c, out)
+        broadcast_func(low_c, high_c, out)
 
         return out
 
@@ -2450,8 +2468,7 @@ cdef class _MKLRandomState:
             if low >= high:
                 raise ValueError("low >= high")
 
-            with self.lock:
-                ret = randfunc(low, high - 1, size)
+            ret = randfunc(low, high - 1, size)
 
             if size is None and dtype in (bool, int):
                 return dtype(ret)
diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index d8dfd937..3559bfea 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -27,6 +27,7 @@
 import sys
 import sysconfig
 import threading
+from collections import Counter
 
 import numpy as np
 import pytest
@@ -62,6 +63,23 @@ def wrapped(i):
     assert not errors
 
 
+def _draw_concurrently(rs, call, k):
+    # k threads each draw once, released together by a barrier.
+    out = [None] * k
+    barrier = threading.Barrier(k)
+
+    def body(i):
+        barrier.wait()
+        out[i] = repr(call(rs))
+
+    threads = [threading.Thread(target=body, args=(i,)) for i in range(k)]
+    for t in threads:
+        t.start()
+    for t in threads:
+        t.join()
+    return Counter(out)
+
+
 def test_concurrent_sampling_per_instance():
     # Each thread owns a private MKLRandomState seeded identically, so the
     # per-instance lock + `nogil` sampling must reproduce the single-threaded
@@ -112,27 +130,31 @@ def worker(_i):
     assert not mkl_random.is_patched()
 
 
-def test_concurrent_multinormal_cholesky_shared():
-    # A data race on the shared stream reuses values; assert few duplicates.
-    n_threads = 8
-    mean = np.zeros(1)
-    ch = np.eye(1)
-    rng = mkl_random.MKLRandomState(12345)
-    chunks = [None] * n_threads
-
-    def worker(i):
-        parts = [
-            rng.multinormal_cholesky(mean, ch, size=4000).ravel()
-            for _ in range(25)
-        ]
-        chunks[i] = np.concatenate(parts)
+_MULTISET_CALLS = {
+    "normal": lambda rs: rs.normal(),
+    "poisson": lambda rs: rs.poisson(3.0),
+    "randint": lambda rs: rs.randint(0, 2**30),
+    "_rand_int32": lambda rs: rs._rand_int32(0, 2**30, None),
+    "multinomial": lambda rs: rs.multinomial(8, [0.25] * 4),
+    "mvn_cholesky": lambda rs: rs.multinormal_cholesky(np.zeros(3), np.eye(3)),
+}
 
-    _run_on_threads(worker, n_threads)
 
-    allvals = np.concatenate(chunks)
-    assert np.all(np.isfinite(allvals))
-    dup_frac = 1.0 - np.unique(allvals).size / allvals.size
-    assert dup_frac < 0.01, f"shared stream corrupted: {dup_frac:.3%} dups"
+@pytest.mark.skipif(
+    not FREE_THREADED, reason="race only manifests without the GIL"
+)
+@pytest.mark.parametrize(
+    "call", _MULTISET_CALLS.values(), ids=list(_MULTISET_CALLS)
+)
+def test_shared_stream_multiset_invariant(call):
+    # Concurrent draws must match the serial multiset; a mismatch = race.
+    k, rounds, seed = 32, 20, 777
+    rs = mkl_random.MKLRandomState(seed)
+    rs.seed(seed)
+    ref = Counter(repr(call(rs)) for _ in range(k))
+    for _ in range(rounds):
+        rs.seed(seed)
+        assert _draw_concurrently(rs, call, k) == ref
 
 
 def test_shuffle_reentrancy():

From 8b746d4e33bf2cbccae053bbb3b98197a5184801 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 9 Sep 2026 04:14:39 -0700
Subject: [PATCH 13/30] Make get_state atomic against concurrent BRNG changes

---
 mkl_random/mklrand.pyx                 | 19 +++++++-----
 mkl_random/tests/test_freethreading.py | 40 ++++++++++++++++++++++++++
 2 files changed, 52 insertions(+), 7 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index 2a0fcfac..548268a1 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -1757,16 +1757,21 @@ cdef class _MKLRandomState:
         MKL Documentation: https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html  # no-cython-lint
 
         """
-        cdef int state_buffer_size
+        cdef int state_buffer_size, cur_size
         cdef int brng_id
         cdef void *bytesPtr
 
-        with self.lock:
-            state_buffer_size = irk_get_stream_size(self.internal_state)
-        bytestring = empty_py_bytes(state_buffer_size, &bytesPtr)
-        with self.lock:
-            brng_id = irk_get_brng_mkl(self.internal_state)
-            irk_get_state_mkl(self.internal_state, <char *>bytesPtr)
+        # The BRNG size can change between sizing and saving, so retry.
+        while True:
+            with self.lock:
+                state_buffer_size = irk_get_stream_size(self.internal_state)
+            bytestring = empty_py_bytes(state_buffer_size, &bytesPtr)
+            with self.lock:
+                cur_size = irk_get_stream_size(self.internal_state)
+                if cur_size == state_buffer_size:
+                    brng_id = irk_get_brng_mkl(self.internal_state)
+                    irk_get_state_mkl(self.internal_state, <char *>bytesPtr)
+                    break
 
         brng_name = _brng_id_to_name(brng_id)
         if legacy:
diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index 3559bfea..ed86c394 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -24,6 +24,7 @@
 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 import os
+import subprocess
 import sys
 import sysconfig
 import threading
@@ -157,6 +158,45 @@ def test_shared_stream_multiset_invariant(call):
         assert _draw_concurrently(rs, call, k) == ref
 
 
+_GET_STATE_RACE = """
+import threading
+import mkl_random
+rs = mkl_random.MKLRandomState(1, brng="MRG32K3A")
+
+def flip():
+    for _ in range(20000):
+        rs.seed(1, brng="MRG32K3A")
+        rs.seed(1, brng="SFMT19937")
+
+def grab():
+    for _ in range(20000):
+        rs.get_state()
+
+ts = [threading.Thread(target=grab) for _ in range(3)]
+ts.append(threading.Thread(target=flip))
+for t in ts:
+    t.start()
+for t in ts:
+    t.join()
+"""
+
+
+@pytest.mark.skipif(
+    not FREE_THREADED, reason="race only manifests without the GIL"
+)
+def test_get_state_race_no_heap_overflow():
+    # get_state racing a BRNG change must not overflow the buffer (a crash).
+    env = dict(os.environ, MKL_NUM_THREADS="1", PYTHONMALLOC="debug")
+    proc = subprocess.run(
+        [sys.executable, "-c", _GET_STATE_RACE],
+        env=env,
+        timeout=120,
+        capture_output=True,
+        text=True,
+    )
+    assert proc.returncode == 0, proc.stderr[-2000:]
+
+
 def test_shuffle_reentrancy():
     # shuffle must not hold the lock across a user callback.
     rs = mkl_random.MKLRandomState(1)

From 4a0fb519ac6185b4d247ac240893333cad7b9c7e Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 9 Sep 2026 04:20:18 -0700
Subject: [PATCH 14/30] Fix MKL stream leak in irk_set_state_mkl

---
 mkl_random/src/randomkit.cpp | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/mkl_random/src/randomkit.cpp b/mkl_random/src/randomkit.cpp
index 612f8acc..7fcf2fb2 100644
--- a/mkl_random/src/randomkit.cpp
+++ b/mkl_random/src/randomkit.cpp
@@ -289,9 +289,20 @@ void irk_get_state_mkl(irk_state *state, char *buf)
 
 int irk_set_state_mkl(irk_state *state, char *buf)
 {
-    int err = vslLoadStreamM(&(state->stream), buf);
+    // vslLoadStreamM allocates a new stream
+    // free the old one to avoid a leak
+    VSLStreamStatePtr stream_loc = NULL;
+    int err = vslLoadStreamM(&stream_loc, buf);
 
-    return (err == VSL_STATUS_OK) ? 0 : 1;
+    if (err != VSL_STATUS_OK) {
+        return 1;
+    }
+    if (state->stream) {
+        vslDeleteStream(&(state->stream));
+    }
+    state->stream = stream_loc;
+
+    return 0;
 }
 
 int irk_leapfrog_stream_mkl(irk_state *state,

From 5d111e52f78e5ddb42879f75ed64445d1255674f Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 9 Sep 2026 04:27:09 -0700
Subject: [PATCH 15/30] Make GIL-import test robust to PYTHON_GIL

---
 mkl_random/tests/test_freethreading.py | 26 ++++++++++++++++++--------
 1 file changed, 18 insertions(+), 8 deletions(-)

diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index ed86c394..63f6f92d 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -30,14 +30,12 @@
 import threading
 from collections import Counter
 
-import numpy as np
-import pytest
-
-# Oversubscription: MKL spawns its own thread pool per calling thread, so
-# generating from many Python threads concurrently can spawn far more OS
-# threads than cores. Cap it before mkl_random/MKL initialize.
+# Cap MKL threads before numpy (may init MKL).
 os.environ.setdefault("MKL_NUM_THREADS", "1")
 
+import numpy as np  # noqa: E402
+import pytest  # noqa: E402
+
 import mkl_random  # noqa: E402
 
 FREE_THREADED = bool(sysconfig.get_config_var("Py_GIL_DISABLED"))
@@ -216,8 +214,20 @@ def run():
     assert done.wait(timeout=30), "shuffle deadlocked on re-entrant callback"
 
 
+_GIL_CHECK = "import sys, mkl_random; assert not sys._is_gil_enabled()"
+
+
 @pytest.mark.skipif(
     not FREE_THREADED, reason="requires a free-threaded CPython build"
 )
-def test_gil_not_reenabled_on_import():
-    assert not sys._is_gil_enabled()  # pylint: disable=no-member
+def test_import_does_not_reenable_gil():
+    # Import in a clean subprocess (no forced PYTHON_GIL); GIL must stay off.
+    env = {k: v for k, v in os.environ.items() if k != "PYTHON_GIL"}
+    proc = subprocess.run(
+        [sys.executable, "-c", _GIL_CHECK],
+        env=env,
+        timeout=60,
+        capture_output=True,
+        text=True,
+    )
+    assert proc.returncode == 0, proc.stderr[-2000:]

From 6e13efd80283d8d9b12d515275bb88049d8d8808 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 9 Sep 2026 04:30:52 -0700
Subject: [PATCH 16/30] Remove dead dummy_threading import fallback

---
 mkl_random/mklrand.pyx | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index 548268a1..4d1a6215 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -546,13 +546,10 @@ if (r < 0):
 import operator
 import warnings
 from collections.abc import Sequence
+from threading import Lock
 
 import numpy as np
 
-try:
-    from threading import Lock
-except ImportError:
-    from dummy_threading import Lock
 
 cdef object vec_cont0_array(
     irk_state *state, irk_cont0_vec func, object size, object lock

From 5d6d26deb819f971de3c39d2d17edcdfe22d3c1b Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 9 Sep 2026 04:36:45 -0700
Subject: [PATCH 17/30] Add TODO for uninitialized vslSaveStreamM save bytes

---
 mkl_random/src/randomkit.cpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/mkl_random/src/randomkit.cpp b/mkl_random/src/randomkit.cpp
index 7fcf2fb2..ff7aec33 100644
--- a/mkl_random/src/randomkit.cpp
+++ b/mkl_random/src/randomkit.cpp
@@ -278,6 +278,8 @@ int irk_get_stream_size(irk_state *state)
 
 void irk_get_state_mkl(irk_state *state, char *buf)
 {
+    // TODO: vslSaveStreamM leaves a few bytes uninitialized (e.g. MT19937
+    // offsets 6,7,14,15); check if oneMKL reserves them, else zero buf.
     int err = vslSaveStreamM(state->stream, buf);
 
     if (err != VSL_STATUS_OK) {

From 548d6255762fe5cc49453473ae78e3fec722de19 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 9 Sep 2026 04:54:08 -0700
Subject: [PATCH 18/30] Fix memory leak in irk_logseries_vec

---
 mkl_random/src/mkl_distributions.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/mkl_random/src/mkl_distributions.cpp b/mkl_random/src/mkl_distributions.cpp
index 01cd2564..5eb6071e 100644
--- a/mkl_random/src/mkl_distributions.cpp
+++ b/mkl_random/src/mkl_distributions.cpp
@@ -1590,6 +1590,7 @@ void irk_logseries_vec(irk_state *state,
     }
 
     mkl_free(Vvec);
+    mkl_free(Uvec);
 }
 
 /* samples discrete uniforms from [low, high) */

From c14aa6c5eb133ce52acb496a725aec5fe9da4fe6 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Wed, 9 Sep 2026 04:54:25 -0700
Subject: [PATCH 19/30] Emit patch print/warn outside the instance lock

---
 mkl_random/_patch_numpy.py             | 60 +++++++++++++++-----------
 mkl_random/tests/test_freethreading.py | 16 +++++++
 2 files changed, 50 insertions(+), 26 deletions(-)

diff --git a/mkl_random/_patch_numpy.py b/mkl_random/_patch_numpy.py
index 9879a5b7..38ecaa24 100644
--- a/mkl_random/_patch_numpy.py
+++ b/mkl_random/_patch_numpy.py
@@ -67,40 +67,48 @@ def _restore_func(self, name, verbose=False):
     def do_patch(self, verbose=False):
         with self._lock:
             local_count = getattr(self._tls, "local_count", 0)
-            if self._patch_count == 0:
-                if verbose:
-                    print(
-                        "Now patching NumPy random submodule with mkl_random "
-                        "NumPy interface."
-                    )
-                    print(
-                        "Please direct bug reports to "
-                        "https://github.com/IntelPython/mkl_random"
-                    )
+            first = self._patch_count == 0
+            if first:
                 for f in self._patched_functions:
                     self._register_func(f, getattr(_nrand, f))
             self._patch_count += 1
             self._tls.local_count = local_count + 1
+        # print outside the lock: it can run arbitrary Python.
+        if verbose and first:
+            print(
+                "Now patching NumPy random submodule with mkl_random "
+                "NumPy interface."
+            )
+            print(
+                "Please direct bug reports to "
+                "https://github.com/IntelPython/mkl_random"
+            )
 
     def do_restore(self, verbose=False):
+        restored = None
         with self._lock:
             local_count = getattr(self._tls, "local_count", 0)
-            if local_count <= 0:
-                warnings.warn(
-                    "restore_numpy_random called more times than "
-                    "patch_numpy_random in this thread.",
-                    RuntimeWarning,
-                    stacklevel=2,
-                )
-                return
-            self._tls.local_count -= 1
-            self._patch_count -= 1
-            if self._patch_count == 0:
-                if verbose:
-                    print("Now restoring original NumPy random submodule.")
-                for name in tuple(self._restore_dict):
-                    self._restore_func(name, verbose=verbose)
-                self._restore_dict.clear()
+            imbalanced = local_count <= 0
+            if not imbalanced:
+                self._tls.local_count -= 1
+                self._patch_count -= 1
+                if self._patch_count == 0:
+                    restored = tuple(self._restore_dict)
+                    for name in restored:
+                        self._restore_func(name)
+                    self._restore_dict.clear()
+        # warn/print outside the lock: they can run arbitrary Python.
+        if imbalanced:
+            warnings.warn(
+                "restore_numpy_random called more times than "
+                "patch_numpy_random in this thread.",
+                RuntimeWarning,
+                stacklevel=2,
+            )
+        elif verbose and restored is not None:
+            print("Now restoring original NumPy random submodule.")
+            for name in restored:
+                print(f"found and restoring {name}...")
 
     def is_patched(self):
         with self._lock:
diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index 63f6f92d..549a411c 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -28,6 +28,7 @@
 import sys
 import sysconfig
 import threading
+import warnings
 from collections import Counter
 
 # Cap MKL threads before numpy (may init MKL).
@@ -214,6 +215,21 @@ def run():
     assert done.wait(timeout=30), "shuffle deadlocked on re-entrant callback"
 
 
+def test_patch_restore_reentrancy():
+    # do_restore must warn outside the lock
+    done = threading.Event()
+
+    def run():
+        with warnings.catch_warnings():
+            warnings.simplefilter("always")
+            warnings.showwarning = lambda *a, **k: mkl_random.is_patched()
+            mkl_random.restore_numpy_random()  # imbalanced -> warns
+        done.set()
+
+    threading.Thread(target=run, daemon=True).start()
+    assert done.wait(timeout=30), "patch restore deadlocked in warn callback"
+
+
 _GIL_CHECK = "import sys, mkl_random; assert not sys._is_gil_enabled()"
 
 

From c322bbfae698f5f7d8565fd5bce2b9e79ffebe6f Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 03:28:34 -0700
Subject: [PATCH 20/30] Fix _seed_impl re-entrancy by coercing the seed outside
 the lock

---
 mkl_random/mklrand.pyx                 | 44 +++++++++++++++-----------
 mkl_random/tests/test_freethreading.py | 19 +++++++++++
 2 files changed, 44 insertions(+), 19 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index 4d1a6215..2ab82c32 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -1638,32 +1638,22 @@ cdef class _MKLRandomState:
         cdef irk_brng_t brng_token = MT19937
         cdef unsigned int stream_id
         cdef cnp.ndarray obj "arrayObject_obj"
+        cdef bint use_array = False
 
         if (brng):
             # Parse before the lock to avoid warn
             brng_token, stream_id = _parse_brng_argument(brng)
 
-        with self.lock:
-            if not brng:
-                # Reads state->stream, which a concurrent seed can free.
-                brng_token = <irk_brng_t> irk_get_brng_and_stream_mkl(
-                    self.internal_state, &stream_id
-                )
+        # Coerce the seed before the lock: operator.index/np.asarray/astype
+        # can run user code that re-enters the generator.
+        idx = 0
+        if seed is not None:
             try:
-                if seed is None:
-                    _errcode = irk_randomseed_mkl(
-                        self.internal_state, brng_token, stream_id
-                    )
-                else:
-                    idx = operator.index(seed)
-                    if idx > int(2**32 - 1) or idx < 0:
-                        raise ValueError(
-                            "Seed must be between 0 and 4294967295"
-                        )
-                    irk_seed_mkl(
-                        self.internal_state, idx, brng_token, stream_id
-                    )
+                idx = operator.index(seed)
+                if idx > int(2**32 - 1) or idx < 0:
+                    raise ValueError("Seed must be between 0 and 4294967295")
             except TypeError:
+                use_array = True
                 obj = np.asarray(seed)
                 if obj.size == 0:
                     raise ValueError("Seed must be non-empty")
@@ -1674,6 +1664,18 @@ cdef class _MKLRandomState:
                 if ((obj > int(2**32 - 1)) | (obj < 0)).any():
                     raise ValueError("Seed must be between 0 and 4294967295")
                 obj = obj.astype("uint32", casting="unsafe", order="C")
+
+        with self.lock:
+            if not brng:
+                # Reads state->stream, which a concurrent seed can free.
+                brng_token = <irk_brng_t> irk_get_brng_and_stream_mkl(
+                    self.internal_state, &stream_id
+                )
+            if seed is None:
+                _errcode = irk_randomseed_mkl(
+                    self.internal_state, brng_token, stream_id
+                )
+            elif use_array:
                 irk_seed_mkl_array(
                     self.internal_state,
                     <unsigned int *>cnp.PyArray_DATA(obj),
@@ -1681,6 +1683,10 @@ cdef class _MKLRandomState:
                     brng_token,
                     stream_id
                 )
+            else:
+                irk_seed_mkl(
+                    self.internal_state, idx, brng_token, stream_id
+                )
 
     def seed(self, seed=None, brng=None):
         """
diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index 549a411c..c74f4264 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -230,6 +230,25 @@ def run():
     assert done.wait(timeout=30), "patch restore deadlocked in warn callback"
 
 
+def test_seed_reentrancy():
+    # A re-entrant __index__ on the seed must not deadlock
+    rs = mkl_random.MKLRandomState(1)
+
+    class ReSeed:
+        def __index__(self):
+            rs.uniform(size=1)
+            return 42
+
+    done = threading.Event()
+
+    def run():
+        rs.seed(ReSeed())
+        done.set()
+
+    threading.Thread(target=run, daemon=True).start()
+    assert done.wait(timeout=30), "seed deadlocked on re-entrant __index__"
+
+
 _GIL_CHECK = "import sys, mkl_random; assert not sys._is_gil_enabled()"
 
 

From 5fba7304de3794367af1d612e324eba664420683 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 03:33:51 -0700
Subject: [PATCH 21/30] Strengthen free-threading tests per review

---
 mkl_random/tests/test_freethreading.py | 32 ++++++++------------------
 1 file changed, 10 insertions(+), 22 deletions(-)

diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index c74f4264..8ca1c59e 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -63,6 +63,12 @@ def wrapped(i):
     assert not errors
 
 
+def _sample_key(x):
+    # Exact key: repr truncates float arrays, so use full-precision bytes.
+    a = np.asarray(x)
+    return (a.dtype.str, a.shape, a.tobytes())
+
+
 def _draw_concurrently(rs, call, k):
     # k threads each draw once, released together by a barrier.
     out = [None] * k
@@ -70,7 +76,7 @@ def _draw_concurrently(rs, call, k):
 
     def body(i):
         barrier.wait()
-        out[i] = repr(call(rs))
+        out[i] = _sample_key(call(rs))
 
     threads = [threading.Thread(target=body, args=(i,)) for i in range(k)]
     for t in threads:
@@ -81,9 +87,8 @@ def body(i):
 
 
 def test_concurrent_sampling_per_instance():
-    # Each thread owns a private MKLRandomState seeded identically, so the
-    # per-instance lock + `nogil` sampling must reproduce the single-threaded
-    # result exactly regardless of concurrency.
+    # Independent instances: the same seed in each thread must reproduce the
+    # single-threaded result regardless of concurrency.
     n_threads = 4
     size = 10**5 + 1  # large enough that per-thread nogil sampling overlaps
     seed = 1234
@@ -99,23 +104,6 @@ def worker(i):
         np.testing.assert_array_equal(r, expected)
 
 
-def test_concurrent_shared_singleton():
-    # Module-level functions share a single lock-guarded RandomState. Hammering
-    # it from many threads must not corrupt state, crash, or return garbage.
-    n_threads = 8
-    size = 10**5 + 1
-    results = [None] * n_threads
-
-    def worker(i):
-        results[i] = mkl_random.uniform(size=size)
-
-    _run_on_threads(worker, n_threads)
-
-    for r in results:
-        assert r.shape == (size,)
-        assert np.all(np.isfinite(r))
-
-
 def test_concurrent_patch_restore():
     n_threads = 8
     n_iters = 20
@@ -151,7 +139,7 @@ def test_shared_stream_multiset_invariant(call):
     k, rounds, seed = 32, 20, 777
     rs = mkl_random.MKLRandomState(seed)
     rs.seed(seed)
-    ref = Counter(repr(call(rs)) for _ in range(k))
+    ref = Counter(_sample_key(call(rs)) for _ in range(k))
     for _ in range(rounds):
         rs.seed(seed)
         assert _draw_concurrently(rs, call, k) == ref

From dd3b15d2623759b9a234efcecad69182c9a65c1c Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 03:59:58 -0700
Subject: [PATCH 22/30] Hash the brng name outside the lock in set_state

---
 mkl_random/mklrand.pyx                 |  4 +++-
 mkl_random/tests/test_freethreading.py | 20 ++++++++++++++++++++
 2 files changed, 23 insertions(+), 1 deletion(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index 2ab82c32..17de80ba 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -1876,6 +1876,8 @@ cdef class _MKLRandomState:
                 "basic number generator algorithm must be one of ['"
                 + "', '".join(_brng_dict.keys()) + "']"
             )
+        # Hash the name outside the lock (re-entrant __hash__)
+        expected_brng = _brng_dict[algorithm_name]
 
         stream_buf = state[1]
         if not is_bytes_object(stream_buf):
@@ -1888,7 +1890,7 @@ cdef class _MKLRandomState:
             if(err):
                 raise ValueError("The stream state buffer is corrupted")
             brng_id = irk_get_brng_mkl(self.internal_state)
-            if (_brng_dict[algorithm_name] != brng_id):
+            if (expected_brng != brng_id):
                 raise ValueError(
                     "The algorithm name does not match content of the buffer"
                 )
diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index 8ca1c59e..7cdd8d8c 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -237,6 +237,26 @@ def run():
     assert done.wait(timeout=30), "seed deadlocked on re-entrant __index__"
 
 
+def test_set_state_reentrancy():
+    # A re-entrant __hash__ on the brng name must not deadlock.
+    rs = mkl_random.MKLRandomState(1)
+    st = rs.get_state()
+
+    class ReStr(str):
+        def __hash__(self):
+            rs.uniform(size=1)
+            return str.__hash__(self)
+
+    done = threading.Event()
+
+    def run():
+        rs.set_state((ReStr(st[0]), st[1]))
+        done.set()
+
+    threading.Thread(target=run, daemon=True).start()
+    assert done.wait(timeout=30), "set_state deadlocked on re-entrant __hash__"
+
+
 _GIL_CHECK = "import sys, mkl_random; assert not sys._is_gil_enabled()"
 
 

From 2c8907d4eddb44d7872d50ac12fef92897ed960c Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 05:31:19 -0700
Subject: [PATCH 23/30] Grow the get_state buffer until it fits

---
 mkl_random/mklrand.pyx | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index 17de80ba..1b77635a 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -1760,21 +1760,24 @@ cdef class _MKLRandomState:
         MKL Documentation: https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html  # no-cython-lint
 
         """
-        cdef int state_buffer_size, cur_size
+        cdef int state_buffer_size = 0
+        cdef int cur_size
         cdef int brng_id
         cdef void *bytesPtr
 
-        # The BRNG size can change between sizing and saving, so retry.
+        bytestring = None
+        # Reseed can change the size; grow until the buffer fits, then save
         while True:
-            with self.lock:
-                state_buffer_size = irk_get_stream_size(self.internal_state)
-            bytestring = empty_py_bytes(state_buffer_size, &bytesPtr)
             with self.lock:
                 cur_size = irk_get_stream_size(self.internal_state)
-                if cur_size == state_buffer_size:
+                if cur_size <= state_buffer_size:
                     brng_id = irk_get_brng_mkl(self.internal_state)
                     irk_get_state_mkl(self.internal_state, <char *>bytesPtr)
                     break
+            state_buffer_size = cur_size
+            bytestring = empty_py_bytes(state_buffer_size, &bytesPtr)
+        if cur_size != state_buffer_size:
+            bytestring = bytestring[:cur_size]
 
         brng_name = _brng_id_to_name(brng_id)
         if legacy:

From 3fd62ab0b48ba7f470aae14e4d6da211f2770952 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 05:44:20 -0700
Subject: [PATCH 24/30] Extend test to cover all size=None distribution

---
 mkl_random/tests/test_freethreading.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index 7cdd8d8c..f9eae8b8 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -125,6 +125,14 @@ def worker(_i):
     "_rand_int32": lambda rs: rs._rand_int32(0, 2**30, None),
     "multinomial": lambda rs: rs.multinomial(8, [0.25] * 4),
     "mvn_cholesky": lambda rs: rs.multinormal_cholesky(np.zeros(3), np.eye(3)),
+    "random_sample": lambda rs: rs.random_sample(),
+    "exponential": lambda rs: rs.exponential(),
+    "triangular": lambda rs: rs.triangular(0.0, 1.0, 2.0),
+    "tomaxint": lambda rs: rs.tomaxint(),
+    "binomial": lambda rs: rs.binomial(10, 0.5),
+    "negative_binomial": lambda rs: rs.negative_binomial(5, 0.5),
+    "hypergeometric": lambda rs: rs.hypergeometric(10, 10, 5),
+    "zipf": lambda rs: rs.zipf(2.0),
 }
 
 

From fd1b15682d1ca1fdd8f3d9b9a184aa884b6e5841 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 05:55:07 -0700
Subject: [PATCH 25/30] Verify Py_GIL_DISABLED matches the target in build_pip

---
 .github/workflows/build_pip.yml | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/.github/workflows/build_pip.yml b/.github/workflows/build_pip.yml
index 0ed45f2a..4bddec40 100644
--- a/.github/workflows/build_pip.yml
+++ b/.github/workflows/build_pip.yml
@@ -11,6 +11,7 @@ permissions:
 
 jobs:
   build:
+    name: build (${{ matrix.python_tag || matrix.python }}${{ matrix.use_pre != '' && ', pre' || '' }})
     runs-on: ubuntu-latest
     defaults:
       run:
@@ -22,16 +23,22 @@ jobs:
         use_pre: ["", "--pre"]
         include:
           - python: "3.14"
+            python_tag: "3.14"
             python_spec: "3.14.* *_cp314"
             use_pre: ""
           - python: "3.14"
+            python_tag: "3.14"
             python_spec: "3.14.* *_cp314"
             use_pre: "--pre"
           - python: "3.14"
+            python_tag: "3.14t"
             python_spec: "3.14.* *_cp314t"
+            free_threaded: "true"
             use_pre: ""
           - python: "3.14"
+            python_tag: "3.14t"
             python_spec: "3.14.* *_cp314t"
+            free_threaded: "true"
             use_pre: "--pre"
 
     steps:
@@ -51,6 +58,23 @@ jobs:
         run: |
           conda install mkl-devel "python=${{ matrix.python_spec || matrix.python }}"
 
+      - name: Verify GIL vs free-threaded ABI
+        env:
+          EXPECT_FREE_THREADED: ${{ matrix.free_threaded || '' }}
+        run: |
+          python - <<'PY'
+          import os, sys, sysconfig
+          gil_disabled = bool(sysconfig.get_config_var("Py_GIL_DISABLED"))
+          expected = bool(os.environ.get("EXPECT_FREE_THREADED"))
+          assert gil_disabled == expected, (
+              f"ABI mismatch: Py_GIL_DISABLED={gil_disabled}, expected "
+              f"free-threaded={expected}; the conda ABI swap did not take effect"
+          )
+          if expected:
+              assert not sys._is_gil_enabled(), "GIL re-enabled on free-threaded build"
+          print(f"OK: Py_GIL_DISABLED={gil_disabled}")
+          PY
+
       - name: Build with pip
         run: |
           pip install --no-cache-dir meson-python ninja cmake "cython>=3.1.0"

From 90e9cb2aec67738fac1913c6035a0073fd1d5a33 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 06:15:55 -0700
Subject: [PATCH 26/30] Add a thread-safety section to the how-to guide

---
 CHANGELOG.md           |  3 +++
 docs/source/how_to.rst | 34 ++++++++++++++++++++++++++++++++++
 2 files changed, 37 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 356619b3..29083d6f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 * Added support for `array_like` (broadcastable) `low`/`high` bounds in `randint` [gh-168](https://github.com/IntelPython/mkl_random/pull/168)
 * Added support for free-threaded (GIL-disabled) CPython builds: the Cython extension is compiled with `freethreading_compatible=True`, so importing `mkl_random` no longer re-enables the GIL [gh-159](https://github.com/IntelPython/mkl_random/pull/159)
+* Added a thread-safety section to the how-to guide for free-threaded Python [gh-159](https://github.com/IntelPython/mkl_random/pull/159)
 
 ### Changed
 * Pinned Cython in the Coverity Scan workflow so generated code stays stable between scans, and added `coverity/README.md` documenting the known Cython-boilerplate false positives and the scan review checklist [gh-164](https://github.com/IntelPython/mkl_random/pull/164)
@@ -16,6 +17,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 * Fixed `uniform` to return a Python `float` for scalar bounds with `size=None` instead of a 0-d array [gh-167](https://github.com/IntelPython/mkl_random/pull/167)
+* Fixed a memory leak in `set_state`, which leaked the previous MKL stream on every call [gh-159](https://github.com/IntelPython/mkl_random/pull/159)
+* Fixed a memory leak in `logseries`, which leaked a temporary buffer on every call [gh-159](https://github.com/IntelPython/mkl_random/pull/159)
 
 ### Removed
 * Removed the `python-gil` constraint from the conda recipes, which pinned `mkl_random` to GIL-enabled Python 3.14 builds [gh-159](https://github.com/IntelPython/mkl_random/pull/159)
diff --git a/docs/source/how_to.rst b/docs/source/how_to.rst
index 985d5748..e7bafadd 100644
--- a/docs/source/how_to.rst
+++ b/docs/source/how_to.rst
@@ -85,3 +85,37 @@ indistinguishable from independent.
 randomness stasistically indistunguishable from independent. To use such families in parallel computation, assign
 difference family generators to different parallel workers and sample those assigned generators in each parallel worker.
 Please refer to "examples/" folder in the `GitHub repo <https://github.com/IntelPython/mkl_random>`_ for more details.
+
+
+Thread safety and free-threaded Python
+---------------------------------------
+
+:mod:`mkl_random` supports free-threaded (GIL-disabled) CPython.
+Every ``MKLRandomState`` instance owns a lock guarding its underlying MKL stream,
+so concurrent calls on a single instance cannot corrupt its state.
+A few properties are worth keeping in mind when sampling from several threads:
+
+* **Prefer one generator per thread.** Sharing one instance across threads is
+  safe -- draws are serialized by the instance lock -- but the order in which
+  concurrent calls interleave is not reproducible from run to run. For
+  reproducible parallel streams, give each thread its own generator and
+  partition the streams as described in the parallel section above
+  (:meth:`skipahead`, :meth:`leapfrog`, or the ``mt2203`` / ``wh`` families).
+
+* **The module-level functions share one global generator.** Calls such as
+  ``mkl_random.normal(...)`` all delegate to a single hidden ``MKLRandomState``.
+  They are thread-safe, but concurrent calls draw from the same stream with a
+  non-reproducible interleaving. Instantiate your own generators when you need
+  control over the streams.
+
+* **Patching NumPy is process-global.** ``patch_numpy_random`` replaces
+  functions on ``numpy.random`` for the whole process, so it affects every
+  thread, not just the calling one. Prefer the ``mkl_random`` context manager
+  for balanced patch / restore handling, but keep in mind that the patch itself
+  is still process-global.
+
+* **In-place shuffles on a shared array are a data race.** ``shuffle`` and other
+  in-place operations mutate the array passed to them. The instance lock
+  protects the random-number stream, not the array itself, so shuffling one
+  array from several threads is a user-level data race. Give each thread its
+  own array.

From 590e460e85e461fcbfe1f4840a88368b7c9bc1d3 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 06:23:06 -0700
Subject: [PATCH 27/30] Add generated mklrand.cpp to .gitignore

---
 .gitignore | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/.gitignore b/.gitignore
index d16f217a..d7e11b01 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,5 +3,8 @@ build/
 .mesonpy-*
 mkl_random.egg-info/
 
+# Cython-generated sources
+mkl_random/mklrand.cpp
+
 # Byte-compiled / optimized / DLL files
 __pycache__/

From a1007259a533d15878a2f1a28b21943ed858f27a Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 12:03:02 -0700
Subject: [PATCH 28/30] Re-lock shuffles typed fast path

---
 mkl_random/mklrand.pyx | 27 ++++++++++++++-------------
 1 file changed, 14 insertions(+), 13 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index 1b77635a..17123dff 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -6808,8 +6808,7 @@ cdef class _MKLRandomState:
 
         u = <cnp.ndarray>self.random_sample(n - 1)
         u_data = <double*>cnp.PyArray_DATA(u)
-        # Indices are already drawn under the lock; the swaps touch no stream
-        # state and run unlocked (locking across the callback would deadlock).
+        # Object/untyped swaps run unlocked: __setitem__ can re-enter the lock.
 
         if type(x) is np.ndarray and x.ndim == 1 and x.size:
             # Fast, statically typed path: shuffle the underlying buffer.
@@ -6824,17 +6823,19 @@ cdef class _MKLRandomState:
             # when the function exits.
             buf = np.empty(itemsize, dtype=np.int8)  # GC'd at function exit
             buf_ptr = cnp.PyArray_BYTES(buf)
-            # We trick gcc into providing a specialized implementation for
-            # the most common case, yielding a ~33% performance improvement.
-            # Note that apparently, only one branch can ever be specialized.
-            if itemsize == sizeof(cnp.npy_intp):
-                self._shuffle_raw(
-                    n, sizeof(cnp.npy_intp), stride, x_ptr, buf_ptr, u_data
-                )
-            else:
-                self._shuffle_raw(
-                    n, itemsize, stride, x_ptr, buf_ptr, u_data
-                )
+            # Pure-C swaps, no callback: safe to lock.
+            with self.lock:
+                # We trick gcc into providing a specialized implementation for
+                # the most common case, yielding a ~33% performance improvement.
+                # Note that apparently, only one branch can ever be specialized.
+                if itemsize == sizeof(cnp.npy_intp):
+                    self._shuffle_raw(
+                        n, sizeof(cnp.npy_intp), stride, x_ptr, buf_ptr, u_data
+                    )
+                else:
+                    self._shuffle_raw(
+                        n, itemsize, stride, x_ptr, buf_ptr, u_data
+                    )
         elif isinstance(x, np.ndarray):
             if x.size == 0:
                 # shuffling is a no-op

From 552d01ba93bc335d0539d4a1ffd0c2ceca34cfe2 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 12:14:51 -0700
Subject: [PATCH 29/30] Make set_state atomic on a BRNG-name mismatch

---
 mkl_random/mklrand.pyx          | 21 +++++++++++----------
 mkl_random/src/randomkit.cpp    | 11 ++++++++++-
 mkl_random/src/randomkit.h      |  3 ++-
 mkl_random/tests/test_random.py | 11 +++++++++++
 4 files changed, 34 insertions(+), 12 deletions(-)

diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx
index 17123dff..af96d88b 100644
--- a/mkl_random/mklrand.pyx
+++ b/mkl_random/mklrand.pyx
@@ -106,7 +106,7 @@ cdef extern from "randomkit.h":
     )
     int irk_get_stream_size(irk_state * state) noexcept nogil
     void irk_get_state_mkl(irk_state * state, char * buf)
-    int irk_set_state_mkl(irk_state * state, char * buf)
+    int irk_set_state_mkl(irk_state * state, char * buf, int expected_brng)
     int irk_get_brng_mkl(irk_state *state) noexcept nogil
     int irk_get_brng_and_stream_mkl(
         irk_state *state, unsigned int * stream_id
@@ -1835,7 +1835,7 @@ cdef class _MKLRandomState:
 
         """
         cdef char *bytes_ptr
-        cdef int brng_id
+        cdef int err
         cdef cnp.ndarray obj "arrayObject_obj"
 
         if isinstance(state, dict):
@@ -1889,14 +1889,15 @@ cdef class _MKLRandomState:
         bytes_ptr = py_bytes_DataPtr(stream_buf)
 
         with self.lock:
-            err = irk_set_state_mkl(self.internal_state, bytes_ptr)
-            if(err):
-                raise ValueError("The stream state buffer is corrupted")
-            brng_id = irk_get_brng_mkl(self.internal_state)
-            if (expected_brng != brng_id):
-                raise ValueError(
-                    "The algorithm name does not match content of the buffer"
-                )
+            err = irk_set_state_mkl(
+                self.internal_state, bytes_ptr, expected_brng
+            )
+        if err == 1:
+            raise ValueError("The stream state buffer is corrupted")
+        if err == 2:
+            raise ValueError(
+                "The algorithm name does not match content of the buffer"
+            )
 
     # Pickling support:
     def __getstate__(self):
diff --git a/mkl_random/src/randomkit.cpp b/mkl_random/src/randomkit.cpp
index ff7aec33..18dfb482 100644
--- a/mkl_random/src/randomkit.cpp
+++ b/mkl_random/src/randomkit.cpp
@@ -289,16 +289,25 @@ void irk_get_state_mkl(irk_state *state, char *buf)
     }
 }
 
-int irk_set_state_mkl(irk_state *state, char *buf)
+int irk_set_state_mkl(irk_state *state, char *buf, int expected_brng)
 {
     // vslLoadStreamM allocates a new stream
     // free the old one to avoid a leak
+    irk_state probe;
     VSLStreamStatePtr stream_loc = NULL;
     int err = vslLoadStreamM(&stream_loc, buf);
 
     if (err != VSL_STATUS_OK) {
         return 1;
     }
+
+    // check the BRNG before publishing; a mismatch leaves state unchanged
+    probe.stream = stream_loc;
+    if (irk_get_brng_mkl(&probe) != expected_brng) {
+        vslDeleteStream(&stream_loc);
+        return 2;
+    }
+
     if (state->stream) {
         vslDeleteStream(&(state->stream));
     }
diff --git a/mkl_random/src/randomkit.h b/mkl_random/src/randomkit.h
index 329be08e..965e00c9 100644
--- a/mkl_random/src/randomkit.h
+++ b/mkl_random/src/randomkit.h
@@ -107,7 +107,8 @@ extern "C"
                                         const unsigned int stream_id);
     extern int irk_get_stream_size(irk_state *state);
     extern void irk_get_state_mkl(irk_state *state, char *buf);
-    extern int irk_set_state_mkl(irk_state *state, char *buf);
+    extern int
+        irk_set_state_mkl(irk_state *state, char *buf, int expected_brng);
     extern int irk_get_brng_mkl(irk_state *state);
     extern int irk_get_brng_and_stream_mkl(irk_state *state,
                                            unsigned int *stream_id);
diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py
index b14fdb30..5fc47cd4 100644
--- a/mkl_random/tests/test_random.py
+++ b/mkl_random/tests/test_random.py
@@ -214,6 +214,17 @@ def test_set_state_negative_binomial(rng_state):
     assert isinstance(v, int)
 
 
+def test_set_state_mismatched_name_is_atomic(rng_state):
+    # A wrong-name buffer must raise and leave the generator unchanged.
+    prng = rng_state.prng
+    ref = prng.tomaxint(8)
+    prng.set_state(rng_state.state)  # rewind
+    other_buf = rnd.MKLRandomState(1, brng="SFMT19937").get_state()[1]
+    with assert_raises(ValueError):
+        prng.set_state(("MT19937", other_buf))
+    assert_equal(prng.tomaxint(8), ref)
+
+
 class RandIntData(NamedTuple):
     rfunc: object
     itype: list

From 65a4a6960061cf614e1351c10c3a2deb1feb9303 Mon Sep 17 00:00:00 2001
From: Vladislav Perevezentsev <vladislav.perevezentsev@intel.com>
Date: Thu, 10 Sep 2026 12:28:48 -0700
Subject: [PATCH 30/30] Make the size=None race test reliably catch a dropped
 lock

---
 mkl_random/tests/test_freethreading.py | 20 ++++++++++++--------
 1 file changed, 12 insertions(+), 8 deletions(-)

diff --git a/mkl_random/tests/test_freethreading.py b/mkl_random/tests/test_freethreading.py
index f9eae8b8..9da72139 100644
--- a/mkl_random/tests/test_freethreading.py
+++ b/mkl_random/tests/test_freethreading.py
@@ -69,21 +69,25 @@ def _sample_key(x):
     return (a.dtype.str, a.shape, a.tobytes())
 
 
-def _draw_concurrently(rs, call, k):
-    # k threads each draw once, released together by a barrier.
-    out = [None] * k
+def _draw_concurrently(rs, call, k, draws):
+    # k threads each draw `draws` samples in a tight loop, released together
+    # by a barrier, so the threads overlap and expose an unlocked stream.
+    results = [None] * k
     barrier = threading.Barrier(k)
 
     def body(i):
         barrier.wait()
-        out[i] = _sample_key(call(rs))
+        results[i] = [_sample_key(call(rs)) for _ in range(draws)]
 
     threads = [threading.Thread(target=body, args=(i,)) for i in range(k)]
     for t in threads:
         t.start()
     for t in threads:
         t.join()
-    return Counter(out)
+    counts = Counter()
+    for local in results:
+        counts.update(local)
+    return counts
 
 
 def test_concurrent_sampling_per_instance():
@@ -144,13 +148,13 @@ def worker(_i):
 )
 def test_shared_stream_multiset_invariant(call):
     # Concurrent draws must match the serial multiset; a mismatch = race.
-    k, rounds, seed = 32, 20, 777
+    k, draws, rounds, seed = 32, 32, 3, 777
     rs = mkl_random.MKLRandomState(seed)
     rs.seed(seed)
-    ref = Counter(_sample_key(call(rs)) for _ in range(k))
+    ref = Counter(_sample_key(call(rs)) for _ in range(k * draws))
     for _ in range(rounds):
         rs.seed(seed)
-        assert _draw_concurrently(rs, call, k) == ref
+        assert _draw_concurrently(rs, call, k, draws) == ref
 
 
 _GET_STATE_RACE = """