diff --git a/.github/workflows/build-with-clang.yml b/.github/workflows/build-with-clang.yml
index 15d15578..87478d91 100644
--- a/.github/workflows/build-with-clang.yml
+++ b/.github/workflows/build-with-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:
       ONEAPI_ROOT: /opt/intel/oneapi
@@ -55,7 +55,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 2a903b50..c2eed53e 100644
--- a/.github/workflows/build-with-standard-clang.yml
+++ b/.github/workflows/build-with-standard-clang.yml
@@ -14,7 +14,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:
@@ -49,7 +49,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 0a145d11..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:
@@ -18,8 +19,27 @@ 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_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:
       - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1
@@ -36,11 +56,28 @@ jobs:
 
       - name: Install MKL
         run: |
-          conda install mkl-devel
+          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
+          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 b76793d6..2dea67fb 100644
--- a/.github/workflows/conda-package-cf.yml
+++ b/.github/workflows/conda-package-cf.yml
@@ -31,6 +31,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
@@ -68,7 +74,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 \
@@ -80,7 +86,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:
@@ -101,6 +107,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\
@@ -130,9 +142,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
@@ -146,12 +158,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:
@@ -171,6 +183,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
@@ -179,7 +197,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"
@@ -202,7 +220,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
@@ -216,9 +234,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
@@ -226,7 +244,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
 
@@ -255,6 +273,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
@@ -263,14 +287,13 @@ 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:
           miniforge-version: latest
           channels: conda-forge
-          activate-environment: ${{ env.TEST_ENV_NAME }}
-          python-version: ${{ matrix.python }}
+          activate-environment: base
 
       - name: Install conda-index
         run: |
@@ -295,7 +318,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 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
@@ -304,9 +327,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
@@ -318,7 +341,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
 
@@ -342,6 +365,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
@@ -368,9 +397,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
@@ -388,7 +417,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 \
@@ -400,7 +429,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:
@@ -419,12 +448,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:
@@ -452,7 +487,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
@@ -462,9 +497,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
@@ -472,7 +507,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 c45bf389..fea31707 100644
--- a/.github/workflows/conda-package.yml
+++ b/.github/workflows/conda-package.yml
@@ -33,6 +33,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
@@ -56,9 +62,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
@@ -87,7 +93,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 \
@@ -99,7 +105,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:
@@ -120,6 +126,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
@@ -148,9 +160,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
@@ -167,12 +179,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:
@@ -181,9 +193,21 @@ 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
+    continue-on-error: ${{ matrix.experimental }}
 
     env:
       CHANNELS: -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels
@@ -192,7 +216,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"
@@ -228,7 +252,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
@@ -242,9 +266,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
@@ -252,7 +276,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
@@ -270,9 +294,21 @@ 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
+    continue-on-error: ${{ matrix.experimental }}
 
     env:
       CHANNELS: -c https://software.repos.intel.com/python/conda -c conda-forge --override-channels
@@ -281,15 +317,13 @@ 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:
-          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: |
@@ -333,7 +367,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 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
@@ -342,9 +376,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
@@ -359,7 +393,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
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__/
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 61e54d3d..29083d6f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,12 +8,20 @@ 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)
+* 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)
 
 ### 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)
 
 ## [1.5.0] (08/12/2026)
 
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:
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
 
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.
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/mklrand.pyx b/mkl_random/mklrand.pyx
index a3259170..af96d88b 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)
@@ -105,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
@@ -545,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
@@ -562,7 +560,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)
@@ -582,7 +581,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)
@@ -678,7 +678,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)
@@ -777,7 +778,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)
@@ -884,7 +886,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)
@@ -909,7 +912,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)
@@ -1008,7 +1012,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)
@@ -1108,7 +1113,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)
@@ -1219,7 +1225,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)
@@ -1243,7 +1250,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"))
@@ -1630,29 +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)
-        else:
-            brng_token = <irk_brng_t> irk_get_brng_and_stream_mkl(
-                self.internal_state, &stream_id
-            )
-        with self.lock:
+
+        # 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")
@@ -1663,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),
@@ -1670,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):
         """
@@ -1743,16 +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
+        cdef int state_buffer_size = 0
+        cdef int 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)
+        bytestring = None
+        # Reseed can change the size; grow until the buffer fits, then save
+        while True:
+            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
+            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:
@@ -1810,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):
@@ -1854,6 +1879,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):
@@ -1862,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 (_brng_dict[algorithm_name] != 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):
@@ -1972,13 +2000,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
 
@@ -1995,13 +2024,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
 
@@ -2018,13 +2048,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
 
@@ -2062,13 +2093,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
 
@@ -2085,13 +2117,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
 
@@ -2108,13 +2141,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
 
@@ -2131,13 +2165,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
 
@@ -2154,13 +2189,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
 
@@ -2177,13 +2213,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
 
@@ -2194,7 +2231,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
             )
@@ -2205,7 +2242,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
             )
@@ -2216,7 +2253,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
             )
@@ -2227,7 +2264,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
             )
@@ -2238,7 +2275,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
             )
@@ -2249,7 +2286,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
             )
@@ -2260,7 +2297,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
             )
@@ -2271,7 +2308,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
             )
@@ -2282,7 +2319,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
             )
@@ -2326,8 +2363,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
 
@@ -2446,8 +2482,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)
@@ -6598,7 +6633,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
 
@@ -6773,6 +6809,7 @@ cdef class _MKLRandomState:
 
         u = <cnp.ndarray>self.random_sample(n - 1)
         u_data = <double*>cnp.PyArray_DATA(u)
+        # 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.
@@ -6787,6 +6824,7 @@ 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)
+            # 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.
@@ -6812,13 +6850,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):
@@ -6829,10 +6866,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,
@@ -6975,7 +7011,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")
@@ -6996,7 +7033,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")
@@ -7079,9 +7117,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)
@@ -7098,9 +7137,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)
@@ -7283,35 +7323,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
 
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) */
diff --git a/mkl_random/src/randomkit.cpp b/mkl_random/src/randomkit.cpp
index 612f8acc..18dfb482 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) {
@@ -287,11 +289,31 @@ 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 err = vslLoadStreamM(&(state->stream), buf);
+    // 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);
 
-    return (err == VSL_STATUS_OK) ? 0 : 1;
+    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));
+    }
+    state->stream = stream_loc;
+
+    return 0;
 }
 
 int irk_leapfrog_stream_mkl(irk_state *state,
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_freethreading.py b/mkl_random/tests/test_freethreading.py
new file mode 100644
index 00000000..9da72139
--- /dev/null
+++ b/mkl_random/tests/test_freethreading.py
@@ -0,0 +1,288 @@
+# 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 subprocess
+import sys
+import sysconfig
+import threading
+import warnings
+from collections import Counter
+
+# 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"))
+
+
+def _run_on_threads(worker, n_threads):
+    # Run worker(i) on n_threads and fail if any thread raised.
+    errors = []
+
+    def wrapped(i):
+        try:
+            worker(i)
+        except Exception as exc:  # pylint: disable=broad-except
+            errors.append(exc)
+
+    threads = [
+        threading.Thread(target=wrapped, args=(i,)) for i in range(n_threads)
+    ]
+    for t in threads:
+        t.start()
+    for t in threads:
+        t.join()
+
+    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, 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()
+        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()
+    counts = Counter()
+    for local in results:
+        counts.update(local)
+    return counts
+
+
+def test_concurrent_sampling_per_instance():
+    # 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
+    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_patch_restore():
+    n_threads = 8
+    n_iters = 20
+
+    def worker(_i):
+        for _ in range(n_iters):
+            mkl_random.patch_numpy_random()
+            mkl_random.restore_numpy_random()
+
+    _run_on_threads(worker, n_threads)
+
+    assert not mkl_random.is_patched()
+
+
+_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)),
+    "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),
+}
+
+
+@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, 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 * draws))
+    for _ in range(rounds):
+        rs.seed(seed)
+        assert _draw_concurrently(rs, call, k, draws) == 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)
+
+    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"
+
+
+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"
+
+
+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__"
+
+
+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()"
+
+
+@pytest.mark.skipif(
+    not FREE_THREADED, reason="requires a free-threaded CPython build"
+)
+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:]
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
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",