Add: Stateful Python random generator API

include/numkong/numkong.h

@@ -29,6 +29,7 @@
 #include "numkong/spatials.h"     // Batched Angular & Euclidean distances, like `nk_angulars_packed_f32`
 #include "numkong/maxsim.h"       // MaxSim: Multi-Vector Maximum Similarity, like `nk_maxsim_packed_f32`
 #include "numkong/trigonometry.h" // Sin, Cos, Atan, like `nk_each_sin_f64`
+#include "numkong/random.h"       // Stateful pseudo-random number generation
 
 #if defined(__cplusplus)
 extern "C" {

include/numkong/random.h

@@ -36,13 +36,65 @@
 #ifndef NK_RANDOM_H
 #define NK_RANDOM_H
 
+#include <math.h>
+#include <stdint.h>
+
 #include "numkong/types.h"
 #include "numkong/cast.h"
 
 #if defined(__cplusplus)
 extern "C" {
 #endif // defined(__cplusplus)
 
+/** @brief Stateful generator used by the language bindings and native callers. */
+typedef struct nk_random_generator_t {
+    uint64_t state;
+} nk_random_generator_t;
+
+/** @brief Initialize a generator from an explicit seed. */
+NK_PUBLIC void nk_random_seed(nk_random_generator_t *generator, uint64_t seed) { generator->state = seed; }
+
+/** @brief Advance the generator and return a deterministic 64-bit value. */
+NK_PUBLIC uint64_t nk_random_next_u64(nk_random_generator_t *generator) {
+    uint64_t value = (generator->state += UINT64_C(0x9E3779B97F4A7C15));
+    value = (value ^ (value >> 30)) * UINT64_C(0xBF58476D1CE4E5B9);
+    value = (value ^ (value >> 27)) * UINT64_C(0x94D049BB133111EB);
+    return value ^ (value >> 31);
+}
+
+/** @brief Return a uniformly distributed double in [0, 1). */
+NK_PUBLIC double nk_random_uniform_f64(nk_random_generator_t *generator) {
+    return (double)(nk_random_next_u64(generator) >> 11) * 0x1.0p-53;
+}
+
+/** @brief Return a uniformly distributed float in [0, 1). */
+NK_PUBLIC float nk_random_uniform_f32(nk_random_generator_t *generator) {
+    return (float)((double)(nk_random_next_u64(generator) >> 40) * 0x1.0p-24);
+}
+
+/** @brief Return a uniformly distributed integer in [0, bound). */
+NK_PUBLIC uint64_t nk_random_bounded_u64(nk_random_generator_t *generator, uint64_t bound) {
+    if (!bound) return nk_random_next_u64(generator);
+    uint64_t const threshold = (uint64_t)(0 - bound) % bound;
+    uint64_t value;
+    do { value = nk_random_next_u64(generator); } while (value < threshold);
+    return value % bound;
+}
+
+/** @brief Return a standard normal sample using the Box-Muller transform. */
+NK_PUBLIC double nk_random_standard_normal_f64(nk_random_generator_t *generator) {
+    double const first = ((double)(nk_random_next_u64(generator) >> 11) + 1.0) * 0x1.0p-53;
+    double const second = nk_random_uniform_f64(generator);
+    double const radius = sqrt(-2.0 * log(first));
+    double const angle = 2.0 * 3.14159265358979323846264338327950288 * second;
+    return radius * cos(angle);
+}
+
+/** @brief Return a normal sample with the requested location and scale. */
+NK_PUBLIC double nk_random_normal_f64(nk_random_generator_t *generator, double location, double scale) {
+    return location + scale * nk_random_standard_normal_f64(generator);
+}
+
 #if defined(__cplusplus)
 } // extern "C"
 #endif // defined(__cplusplus)

python/README.md

@@ -35,6 +35,32 @@ dot = nk.dot(a, b)  # widened accumulation, not same-dtype
 print(dot)
 ```
 
+## Stateful Random Generation
+
+`nk.Generator` provides an independent, reproducible source of random values for pipelines that must not share
+process-global random state. It supports `random`, `uniform`, `normal`, `standard_normal`, and
+`integers`. Every method requires `size` as an integer or tuple and accepts optional `dtype=` and writable
+C-contiguous `out=` buffers.
+
+```python
+import numpy as np
+import numkong as nk
+
+generator = nk.Generator(seed=137)
+images = generator.uniform(0.0, 1.0, size=(8, 224, 224, 3), dtype="float32")
+labels = generator.integers(0, 10, size=32, dtype="int32")
+
+same_images = nk.Generator(seed=137).uniform(0.0, 1.0, size=(8, 224, 224, 3), dtype="float32")
+np.testing.assert_array_equal(np.asarray(images), np.asarray(same_images))
+```
+
+Each generator serializes fills and state access across threads. It releases the GIL while filling native output
+and returns a NumKong `Tensor` unless `out=` is supplied. With `out=`, the method writes in place and returns `None`.
+NumPy arrays, `memoryview` objects, and NumKong tensors are accepted when their shape, dtype, and C-contiguous layout
+match the request.
+
+The implementation uses deterministic scalar kernels. SIMD kernels and additional distributions are not included.
+
 ## Installation
 
 From PyPI:

python/numkong.c

@@ -98,6 +98,7 @@
 #include "each.h"
 #include "mesh.h"
 #include "maxsim.h"
+#include "random.h"
 #include "numpy_interop.h"
 #include "dlpack_interop.h"
 
@@ -1212,6 +1213,7 @@ PyMODINIT_FUNC PyInit__numkong(void) {
     if (PyType_Ready(&PackedMatrixType) < 0) return NULL;
     if (PyType_Ready(&MaxSimPackedMatrixType) < 0) return NULL;
     if (PyType_Ready(&MeshAlignmentResultType) < 0) return NULL;
+    if (PyType_Ready(&RandomGeneratorType) < 0) return NULL;
 
     m = PyModule_Create(&nk_module);
     if (m == NULL) return NULL;
@@ -1264,6 +1266,13 @@ PyMODINIT_FUNC PyInit__numkong(void) {
         return NULL;
     }
 
+    Py_INCREF(&RandomGeneratorType);
+    if (PyModule_AddObject(m, "Generator", (PyObject *)&RandomGeneratorType) < 0) {
+        Py_XDECREF(&RandomGeneratorType);
+        Py_XDECREF(m);
+        return NULL;
+    }
+
     static_capabilities = nk_capabilities();
     nk_configure_thread(static_capabilities);
 

python/numkong/__init__.pyi

@@ -1,4 +1,4 @@
-from typing import Any, Literal, TypeAlias
+from typing import Any, Literal, TypeAlias, overload
 
 # Many annotation features depend on the Python version:
 # - `typing.TypeAlias` type aliases are supported in Python 3.10-3.11.
@@ -369,6 +369,118 @@ class Tensor(memoryview):
         """Return the index of the maximum element, or None if all elements are NaN."""
         ...
 
+class Generator:
+    """Independent stateful pseudo-random number generator."""
+
+    def __init__(self, seed: int | None = None) -> None:
+        """Create a generator with an explicit integer seed."""
+        ...
+
+    @property
+    def state(self) -> int:
+        """Read or replace the generator's reproducible state."""
+        ...
+
+    @state.setter
+    def state(self, value: int) -> None: ...
+    def random(
+        self,
+        size: int | tuple[int, ...],
+        *,
+        dtype: Literal["f32", "float32", "f64", "float64"] = "float32",
+        out: _BufferType | None = None,
+    ) -> Tensor | None:
+        """Generate uniform samples in [0, 1)."""
+        ...
+
+    @overload
+    def uniform(
+        self,
+        low: float,
+        high: float,
+        size: int | tuple[int, ...],
+        *,
+        dtype: Literal["f32", "float32", "f64", "float64"] = "float32",
+        out: _BufferType | None = None,
+    ) -> Tensor | None:
+        """Generate uniform samples in [low, high)."""
+        ...
+
+    @overload
+    def uniform(
+        self,
+        low: float = 0.0,
+        high: float = 1.0,
+        *,
+        size: int | tuple[int, ...],
+        dtype: Literal["f32", "float32", "f64", "float64"] = "float32",
+        out: _BufferType | None = None,
+    ) -> Tensor | None:
+        """Generate uniform samples in [low, high)."""
+        ...
+
+    @overload
+    def normal(
+        self,
+        loc: float,
+        scale: float,
+        size: int | tuple[int, ...],
+        *,
+        dtype: Literal["f32", "float32", "f64", "float64"] = "float32",
+        out: _BufferType | None = None,
+    ) -> Tensor | None:
+        """Generate normal samples with the requested location and scale."""
+        ...
+
+    @overload
+    def normal(
+        self,
+        loc: float = 0.0,
+        scale: float = 1.0,
+        *,
+        size: int | tuple[int, ...],
+        dtype: Literal["f32", "float32", "f64", "float64"] = "float32",
+        out: _BufferType | None = None,
+    ) -> Tensor | None:
+        """Generate normal samples with the requested location and scale."""
+        ...
+
+    def standard_normal(
+        self,
+        size: int | tuple[int, ...],
+        *,
+        dtype: Literal["f32", "float32", "f64", "float64"] = "float32",
+        out: _BufferType | None = None,
+    ) -> Tensor | None:
+        """Generate standard normal samples."""
+        ...
+
+    @overload
+    def integers(
+        self,
+        low: int,
+        high: int | None,
+        size: int | tuple[int, ...],
+        *,
+        dtype: Literal["int8", "int16", "int32", "int64", "uint8", "uint16", "uint32", "uint64"] = "int64",
+        out: _BufferType | None = None,
+    ) -> Tensor | None:
+        """Generate uniformly distributed integers in [low, high)."""
+        ...
+
+    @overload
+    def integers(
+        self,
+        low: int,
+        high: int | None = None,
+        *,
+        size: int | tuple[int, ...],
+        dtype: Literal["int8", "int16", "int32", "int64", "uint8", "uint16", "uint32", "uint64"] = "int64",
+        out: _BufferType | None = None,
+    ) -> Tensor | None:
+        """Generate uniformly distributed integers in [low, high)."""
+        ...
+
 class PackedMatrix:
     """Opaque pre-packed matrix for repeated cross operations.