diff --git a/CHANGELOG.md b/CHANGELOG.md
index 596f483596a..9dee06ee9ca 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Added `dpnp.top_k` to return the `k` largest or smallest elements and their indices along an axis, with a `sorted` keyword to return them in sorted order [#3133](https://github.com/IntelPython/dpnp/pull/3133)
 * Added the `dpnp.random.Generator` API with `default_rng`, `BitGenerator`, `SeedSequence` and the `Philox4x3210` bit generator, sampling on the device with the oneMKL device RNG API [#3157](https://github.com/IntelPython/dpnp/pull/3157)
 * Added the continuous distributions of `numpy.random.Generator`, such as `gamma`, `beta` and `vonmises`, to `dpnp.random.Generator` [#3158](https://github.com/IntelPython/dpnp/pull/3158)
+* Added `integers` and the discrete distributions of `numpy.random.Generator`, such as `poisson`, `binomial` and `hypergeometric`, to `dpnp.random.Generator`, along with `noncentral_chisquare` and `noncentral_f` [#3159](https://github.com/IntelPython/dpnp/pull/3159)
 
 ### Changed
 
diff --git a/doc/reference/random.rst b/doc/reference/random.rst
index 69cf865d7fb..b6b989aa2e8 100644
--- a/doc/reference/random.rst
+++ b/doc/reference/random.rst
@@ -17,19 +17,21 @@ host-side oneMKL RNG.
 
     rng = dpnp.random.default_rng(seed=42)
     x = rng.standard_normal((3, 4))                # on the default device
-    y = rng.random(5, usm_type="shared")
+    y = rng.integers(0, 10, size=5, usm_type="shared")
 
     rng = dpnp.random.default_rng(42, device="cpu")
 
 A generator draws all its samples on one SYCL queue, chosen with the
 ``device`` or ``sycl_queue`` keyword when it is created, and the samples
-depend only on the seed, not on the device. All samples are
+depend only on the seed, not on the device, except that the discrete
+distributions also depend on its support of double precision. All samples are
 :class:`dpnp.ndarray`, zero-dimensional when ``size`` is ``None``.
 
 Differences from NumPy
 ----------------------
 
-* The same seed gives the same stream on every device.
+* The same seed gives the same stream on every device, but for the discrete
+  distributions noted below.
   The streams differ from NumPy's.
 * The floating-point distributions take a ``dtype`` keyword, as in CuPy,
   which defaults to the default floating-point type of the device. All
@@ -39,7 +41,20 @@ Differences from NumPy
   not, since that would wait for the device.
   Invalid ones give undefined values, such as NaN or the smallest ``int64``,
   but the kernels handle them safely: they never loop without end or access
-  memory out of bounds.
+  memory out of bounds. Floating-point arrays given for integer parameters
+  are converted as by :meth:`dpnp.ndarray.astype`, whose result for NaN and
+  out-of-range values is undefined.
+* The discrete distributions return ``int64`` values. They compute in double
+  precision on devices supporting it, and in single precision otherwise.
+  There, the streams differ from the ones of devices with double precision,
+  and some samples are only multiples of a power of two: the ones above
+  ``2**24`` of ``geometric``, ``logseries`` and ``zipf``, and the ones of
+  ``poisson`` and ``binomial`` with a standard deviation above about
+  ``2**22``.
+* ``zipf`` follows the Zipf distribution truncated to the largest ``int64``.
+  NumPy's rejection test rounds to accepting every sample above ``2**53``,
+  which makes its tail too heavy for ``a`` close to 1.
+* ``hypergeometric`` requires ``ngood`` and ``nbad`` below ``10**9``.
 
 Random Generator
 ----------------
@@ -58,6 +73,7 @@ Simple random data
    :toctree: generated/
    :nosignatures:
 
+   Generator.integers
    Generator.random
 
 Distributions
@@ -68,16 +84,24 @@ Distributions
    :nosignatures:
 
    Generator.beta
+   Generator.binomial
    Generator.chisquare
    Generator.exponential
    Generator.f
    Generator.gamma
+   Generator.geometric
    Generator.gumbel
+   Generator.hypergeometric
    Generator.laplace
    Generator.logistic
    Generator.lognormal
+   Generator.logseries
+   Generator.negative_binomial
+   Generator.noncentral_chisquare
+   Generator.noncentral_f
    Generator.normal
    Generator.pareto
+   Generator.poisson
    Generator.power
    Generator.rayleigh
    Generator.standard_cauchy
@@ -90,6 +114,7 @@ Distributions
    Generator.vonmises
    Generator.wald
    Generator.weibull
+   Generator.zipf
 
 Bit generators
 --------------
diff --git a/dpnp/AGENTS.md b/dpnp/AGENTS.md
index 7efb8ae7b9b..f87cde80760 100644
--- a/dpnp/AGENTS.md
+++ b/dpnp/AGENTS.md
@@ -21,7 +21,7 @@
 - Extension calls are ordered with `_manager = dpctl.utils.SequentialOrderManager[exec_q]`: pass `depends=_manager.submitted_events`, then `_manager.add_event_pair(ht_ev, ev)`. Don't add blocking waits.
 - Use `dpnp.default_float_type(...)` or `map_dtype_to_device(...)` for the default float type.
 - Deprecated ufuncs use `DPNPDeprecatedUnaryFunc(..., deprecated_msg=...)`.
-- A `dpnp.random` distribution is a `Generator` method: it calls `self._check_queues(...)` on its array parameters (raises `ExecutionPlacementError`), checks host parameters as NumPy does, and draws through `_cont` or `_float_sample` into `BitGenerator._fill`. Parameters on the device are not checked, as the check waits for the device. Temporaries take the output's `usm_type`.
+- A `dpnp.random` distribution is a `Generator` method: it calls `self._check_queues(...)` on its array parameters (raises `ExecutionPlacementError`), checks host parameters as NumPy does, and draws through `_cont`, `_float_sample` or `_discrete_sample` into `BitGenerator._fill`. Parameters on the device are not checked, as the check waits for the device. Temporaries take the output's `usm_type`.
 - Unsupported NumPy keywords raise `NotImplementedError` (`dpnp.check_limitations`) and are listed under `Limitations` in the docstring.
 - The Array API spec version is set in three places that change together: `__array_api_version__` in `dpnp/tensor/_array_api.py`, `ARRAY_API_TESTS_VERSION` in `.github/workflows/conda-package.yml` and `doc/reference/array_api.rst` (plus the version test in `dpnp/tests/test_ndarray.py`). Inspection utilities come from `__array_namespace_info__` in `dpnp_array_api_info.py`.
 - `dpnp.tensor` replaces `dpctl.tensor`, which no longer exists here; `SyclQueue`, `SyclDevice`, `dpctl.memory` and `SequentialOrderManager` still come from `dpctl`.
diff --git a/dpnp/libtensor/AGENTS.md b/dpnp/libtensor/AGENTS.md
index 0bb184f73f9..a0e98b23381 100644
--- a/dpnp/libtensor/AGENTS.md
+++ b/dpnp/libtensor/AGENTS.md
@@ -23,7 +23,7 @@ A new extension's `CMakeLists.txt` lists its sources and calls `build_dpnp_pybin
 - Guard double, complex<double> and half with `dpnp::tensor::type_utils::validate_type_for_device<T>(q)`; don't use double literals or casts in kernels that must run without fp64.
 - `blas`, `fft`, `lapack` and `sparse` build against either oneMKL or oneMath (`USE_ONEMATH`); `vm` has kernels only with oneMKL: with oneMath it builds empty and `_is_available()` returns `False`; `rng` uses the header-only oneMKL device RNG API and links no oneMKL library.
 - Release the GIL around blocking oneMKL calls (`py::gil_scoped_release`).
-- `rng` values depend only on the engine state and the number of elements, never on the device or the layout chosen for it (`extensions/rng/kernels/sample.hpp`); each engine's `items_per_work_item` (`kernels/engines.hpp`) is part of its stream definition. Parameters on the device are unchecked, so a rejection sampler returns NaN for the inputs that would make it loop forever, including values that are valid but round badly in fp32.
+- `rng` values depend only on the engine state and the number of elements, never on the device or the layout chosen for it (`extensions/rng/kernels/sample.hpp`); each engine's `items_per_work_item` (`kernels/engines.hpp`) is part of its stream definition. Parameters on the device are unchecked, so a rejection sampler returns NaN or the smallest `int64` for the inputs that would make it loop forever, including values that are valid but round badly in fp32. Use `umul_hi` from `kernels/umul_hi.hpp` instead of `sycl::mul_hi`, which DPC++ 2026.2 miscompiles on Intel GPUs.
 - An extension whose kernels may be launched over more than `INT_MAX` work-items passes `LARGE_RANGES` to `build_dpnp_pybind11_ext`, which adds `-fno-sycl-id-queries-fit-in-int` at compile and link time (DPC++ 2026.2 and later need it at link time). `source/CMakeLists.txt` and `dpnp_backend_c` pass it to both `target_compile_options` and `target_link_options`.
 - clang-format 23 formats C++ (`.clang-format`: LLVM base, 4-space indent, 80 columns); gersemi formats CMake at 88 columns, except `dpnp/backend/cmake/Modules/`.
 
diff --git a/dpnp/libtensor/extensions/rng/CMakeLists.txt b/dpnp/libtensor/extensions/rng/CMakeLists.txt
index 7e41305df9d..e2289b19fd9 100644
--- a/dpnp/libtensor/extensions/rng/CMakeLists.txt
+++ b/dpnp/libtensor/extensions/rng/CMakeLists.txt
@@ -29,9 +29,14 @@
 
 set(python_module_name _rng_impl)
 set(_module_src
+    ${CMAKE_CURRENT_SOURCE_DIR}/binomial.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/discrete.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/exponential.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/gamma.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/integers.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/noncentral.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/normal.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/poisson.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/rng_py.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/sample.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/uniform.cpp
diff --git a/dpnp/libtensor/extensions/rng/binomial.cpp b/dpnp/libtensor/extensions/rng/binomial.cpp
new file mode 100644
index 00000000000..58c8bce687e
--- /dev/null
+++ b/dpnp/libtensor/extensions/rng/binomial.cpp
@@ -0,0 +1,44 @@
+//*****************************************************************************
+// Copyright (c) 2026, Intel Corporation
+// All rights reserved.
+//
+// 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 the copyright holder 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 HOLDER 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.
+//*****************************************************************************
+#include <pybind11/pybind11.h>
+
+#include "common.hpp"
+#include "distributions.hpp"
+
+namespace dpnp::extensions::rng
+{
+static discrete_dispatch_table_t binomial_dispatch_table;
+
+void init_binomial(py::module_ m)
+{
+    def_discrete_sampler<rng_ns::BinomialSampler>(
+        m, "_binomial",
+        "Fill the result with binomial samples of n trials of probability p",
+        binomial_dispatch_table);
+}
+} // namespace dpnp::extensions::rng
diff --git a/dpnp/libtensor/extensions/rng/common.hpp b/dpnp/libtensor/extensions/rng/common.hpp
index 1e72cd95017..abdd80042eb 100644
--- a/dpnp/libtensor/extensions/rng/common.hpp
+++ b/dpnp/libtensor/extensions/rng/common.hpp
@@ -81,6 +81,10 @@ typedef sycl::event (*sample_fn_ptr_t)(sycl::queue &,
 using sample_dispatch_table_t =
     sample_fn_ptr_t[rng_ns::num_engines][td_ns::num_types];
 
+// Discrete samplers, whose result is always int64, indexed by the engine and
+// by whether they compute in double (1) or in float (0)
+using discrete_dispatch_table_t = sample_fn_ptr_t[rng_ns::num_engines][2];
+
 namespace detail
 {
 template <typename T>
@@ -138,8 +142,11 @@ sycl::event submit_sample(sycl::queue &exec_q,
 } // namespace detail
 
 // Checks the parameters and the result against SamplerT, then submits the
-// kernel, unless there is nothing to sample.
-template <typename EngineTraitsT, typename SamplerT>
+// kernel, unless there is nothing to sample. ComputeT is the type the
+// sampler computes in, which the device must support.
+template <typename EngineTraitsT,
+          typename SamplerT,
+          typename ComputeT = typename SamplerT::result_type>
 sycl::event sample_impl(sycl::queue &exec_q,
                         const rng_ns::EngineState &state,
                         const std::size_t nelems,
@@ -162,6 +169,7 @@ sycl::event sample_impl(sycl::queue &exec_q,
     }
 
     detail::validate_type_for_device<T>(exec_q);
+    detail::validate_type_for_device<ComputeT>(exec_q);
     detail::validate_types_for_device<param_types>(
         exec_q, std::make_index_sequence<n_params>{});
 
@@ -195,6 +203,12 @@ struct RealTypes
 {
 };
 
+// Types the integer samplers support
+template <typename T>
+struct IntegerTypes : std::bool_constant<std::is_integral_v<T>>
+{
+};
+
 // The raw output of a bit generator, whole or as 32-bit halves
 template <typename T>
 struct RawTypes : std::bool_constant<std::is_same_v<T, std::uint64_t> ||
@@ -238,11 +252,27 @@ void init_sample_dispatch_table(sample_dispatch_table_t &table)
         table[philox]);
 }
 
+template <template <typename> class SamplerT>
+void init_discrete_dispatch_table(discrete_dispatch_table_t &table)
+{
+    static_assert(rng_ns::num_engines == 1,
+                  "Every engine must have a row in the dispatch table");
+
+    using Philox = rng_ns::Philox4x32x10;
+    constexpr int philox = static_cast<int>(rng_ns::engine_id::philox4x32x10);
+    table[philox][0] = sample_impl<Philox, SamplerT<float>, float>;
+    table[philox][1] = sample_impl<Philox, SamplerT<double>, double>;
+}
+
 // Looks up the sampler of `result`'s type, and throws if there is none
 extern sample_fn_ptr_t lookup_sampler(const sample_dispatch_table_t &table,
                                       const int engine,
                                       const dpnp::tensor::usm_ndarray &result);
 
+extern sample_fn_ptr_t lookup_sampler(const discrete_dispatch_table_t &table,
+                                      const int engine,
+                                      const bool double_precision);
+
 extern std::pair<sycl::event, sycl::event>
     py_sample(sycl::queue &exec_q,
               const std::vector<std::uint64_t> &seed,
@@ -281,4 +311,34 @@ void def_sampler(py::module_ &m,
           py::arg("params"), py::arg("result"),
           py::arg("depends") = py::list());
 }
+
+// As def_sampler, for a discrete sampler, whose Python function also takes
+// whether to compute in double precision
+template <template <typename> class SamplerT>
+void def_discrete_sampler(py::module_ &m,
+                          const char *name,
+                          const char *doc,
+                          discrete_dispatch_table_t &table)
+{
+    init_discrete_dispatch_table<SamplerT>(table);
+
+    auto sample_pyapi = [&table](sycl::queue &exec_q, const int engine,
+                                 const std::vector<std::uint64_t> &seed,
+                                 const std::uint64_t counter_lo,
+                                 const std::uint64_t counter_hi,
+                                 const py::sequence &params,
+                                 const dpnp::tensor::usm_ndarray &result,
+                                 const bool double_precision,
+                                 const std::vector<sycl::event> &depends) {
+        const sample_fn_ptr_t fn =
+            lookup_sampler(table, engine, double_precision);
+        return py_sample(exec_q, seed, counter_lo, counter_hi, params, result,
+                         depends, fn);
+    };
+
+    m.def(name, sample_pyapi, doc, py::arg("sycl_queue"), py::arg("engine"),
+          py::arg("seed"), py::arg("counter_lo"), py::arg("counter_hi"),
+          py::arg("params"), py::arg("result"), py::arg("double_precision"),
+          py::arg("depends") = py::list());
+}
 } // namespace dpnp::extensions::rng
diff --git a/dpnp/libtensor/extensions/rng/discrete.cpp b/dpnp/libtensor/extensions/rng/discrete.cpp
new file mode 100644
index 00000000000..a3c05a207ac
--- /dev/null
+++ b/dpnp/libtensor/extensions/rng/discrete.cpp
@@ -0,0 +1,62 @@
+//*****************************************************************************
+// Copyright (c) 2026, Intel Corporation
+// All rights reserved.
+//
+// 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 the copyright holder 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 HOLDER 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.
+//*****************************************************************************
+#include <pybind11/pybind11.h>
+
+#include "common.hpp"
+#include "distributions.hpp"
+
+namespace dpnp::extensions::rng
+{
+static discrete_dispatch_table_t geometric_dispatch_table;
+static discrete_dispatch_table_t logseries_dispatch_table;
+static discrete_dispatch_table_t zipf_dispatch_table;
+static discrete_dispatch_table_t hypergeometric_dispatch_table;
+
+void init_discrete(py::module_ m)
+{
+    def_discrete_sampler<rng_ns::GeometricSampler>(
+        m, "_geometric",
+        "Fill the result with geometric samples of probability p",
+        geometric_dispatch_table);
+
+    def_discrete_sampler<rng_ns::LogseriesSampler>(
+        m, "_logseries",
+        "Fill the result with logarithmic series samples of shape p",
+        logseries_dispatch_table);
+
+    def_discrete_sampler<rng_ns::ZipfSampler>(
+        m, "_zipf", "Fill the result with Zipf samples of exponent a",
+        zipf_dispatch_table);
+
+    def_discrete_sampler<rng_ns::HypergeometricSampler>(
+        m, "_hypergeometric",
+        "Fill the result with hypergeometric samples of nsample draws from "
+        "ngood good and nbad bad items",
+        hypergeometric_dispatch_table);
+}
+} // namespace dpnp::extensions::rng
diff --git a/dpnp/libtensor/extensions/rng/distributions.hpp b/dpnp/libtensor/extensions/rng/distributions.hpp
index 6d41e71b39d..c78e93bc100 100644
--- a/dpnp/libtensor/extensions/rng/distributions.hpp
+++ b/dpnp/libtensor/extensions/rng/distributions.hpp
@@ -38,4 +38,9 @@ void init_uniform(py::module_ m);
 void init_normal(py::module_ m);
 void init_exponential(py::module_ m);
 void init_gamma(py::module_ m);
+void init_noncentral(py::module_ m);
+void init_integers(py::module_ m);
+void init_poisson(py::module_ m);
+void init_binomial(py::module_ m);
+void init_discrete(py::module_ m);
 } // namespace dpnp::extensions::rng
diff --git a/dpnp/libtensor/extensions/rng/integers.cpp b/dpnp/libtensor/extensions/rng/integers.cpp
new file mode 100644
index 00000000000..5174e2f68ee
--- /dev/null
+++ b/dpnp/libtensor/extensions/rng/integers.cpp
@@ -0,0 +1,45 @@
+//*****************************************************************************
+// Copyright (c) 2026, Intel Corporation
+// All rights reserved.
+//
+// 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 the copyright holder 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 HOLDER 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.
+//*****************************************************************************
+#include <pybind11/pybind11.h>
+
+#include "common.hpp"
+#include "distributions.hpp"
+
+namespace dpnp::extensions::rng
+{
+static sample_dispatch_table_t integers_dispatch_table;
+
+void init_integers(py::module_ m)
+{
+    def_sampler<rng_ns::IntegersSampler, IntegerTypes>(
+        m, "_integers",
+        "Fill the result with uniform integers on the closed interval [low, "
+        "high]",
+        integers_dispatch_table);
+}
+} // namespace dpnp::extensions::rng
diff --git a/dpnp/libtensor/extensions/rng/kernels/distributions.hpp b/dpnp/libtensor/extensions/rng/kernels/distributions.hpp
index a39bf73e25e..37ab100746c 100644
--- a/dpnp/libtensor/extensions/rng/kernels/distributions.hpp
+++ b/dpnp/libtensor/extensions/rng/kernels/distributions.hpp
@@ -64,10 +64,14 @@
 #pragma once
 
 #include <cstddef>
+#include <cstdint>
 #include <limits>
+#include <type_traits>
 
 #include <sycl/sycl.hpp>
 
+#include "umul_hi.hpp"
+
 namespace dpnp::kernels::rng
 {
 // A distribution parameter: either one value for every sample, or a
@@ -241,6 +245,67 @@ T beta(BitGenT &bg, StandardNormal<T> &normal, const T a, const T b)
     return (res > T(1)) ? T(1) : res;
 }
 
+// Lemire's method for a uniform integer on [0, rng], rng < 2**32 - 1
+template <typename BitGenT>
+std::uint32_t bounded_lemire_uint32(BitGenT &bg, const std::uint32_t rng)
+{
+    const std::uint32_t rng_excl = rng + 1;
+
+    std::uint64_t m = static_cast<std::uint64_t>(bg.next_uint32()) * rng_excl;
+    std::uint32_t leftover = static_cast<std::uint32_t>(m);
+
+    if (leftover < rng_excl) {
+        const std::uint32_t threshold =
+            (std::numeric_limits<std::uint32_t>::max() - rng) % rng_excl;
+
+        while (leftover < threshold) {
+            m = static_cast<std::uint64_t>(bg.next_uint32()) * rng_excl;
+            leftover = static_cast<std::uint32_t>(m);
+        }
+    }
+    return static_cast<std::uint32_t>(m >> 32);
+}
+
+// Lemire's method for a uniform integer on [0, rng], rng < 2**64 - 1
+template <typename BitGenT>
+std::uint64_t bounded_lemire_uint64(BitGenT &bg, const std::uint64_t rng)
+{
+    const std::uint64_t rng_excl = rng + 1;
+
+    std::uint64_t x = bg.next_uint64();
+    std::uint64_t leftover = x * rng_excl;
+
+    if (leftover < rng_excl) {
+        const std::uint64_t threshold =
+            (std::numeric_limits<std::uint64_t>::max() - rng) % rng_excl;
+
+        while (leftover < threshold) {
+            x = bg.next_uint64();
+            leftover = x * rng_excl;
+        }
+    }
+    return umul_hi(x, rng_excl);
+}
+
+// A uniform integer on [0, rng], as in NumPy's random_bounded_uint64()
+template <typename BitGenT>
+std::uint64_t bounded_uint64(BitGenT &bg, const std::uint64_t rng)
+{
+    if (rng == 0) {
+        return 0;
+    }
+    else if (rng <= std::numeric_limits<std::uint32_t>::max()) {
+        if (rng == std::numeric_limits<std::uint32_t>::max()) {
+            return bg.next_uint32();
+        }
+        return bounded_lemire_uint32(bg, static_cast<std::uint32_t>(rng));
+    }
+    else if (rng == std::numeric_limits<std::uint64_t>::max()) {
+        return bg.next_uint64();
+    }
+    return bounded_lemire_uint64(bg, rng);
+}
+
 template <typename T, typename BitGenT>
 T chisquare(BitGenT &bg, StandardNormal<T> &normal, const T df)
 {
@@ -269,6 +334,215 @@ T f(BitGenT &bg, StandardNormal<T> &normal, const T dfnum, const T dfden)
     return num / (chisquare<T>(bg, normal, dfden) * dfnum);
 }
 
+// The error of Stirling's approximation, log(n!) - log(sqrt(2*pi*n)*(n/e)**n)
+template <typename T>
+T stirlerr(const T n)
+{
+    constexpr T S0 = T(1) / T(12);
+    constexpr T S1 = T(1) / T(360);
+    constexpr T S2 = T(1) / T(1260);
+    constexpr T S3 = T(1) / T(1680);
+    constexpr T S4 = T(1) / T(1188);
+    // log(sqrt(2*pi))
+    constexpr T ln_sqrt_2pi = T(0.918938533204672741780329736406);
+
+    if (n <= T(15)) {
+        return sycl::lgamma(n + T(1)) - (n + T(0.5)) * sycl::log(n) + n -
+               ln_sqrt_2pi;
+    }
+    const T nn = n * n;
+    return (S0 - (S1 - (S2 - (S3 - S4 / nn) / nn) / nn) / nn) / n;
+}
+
+// Loader's deviance term x*log(x/np) + np - x, given d = x - np, which does
+// not cancel for x close to np
+template <typename T>
+T bd0(const T x, const T np, const T d)
+{
+    if (sycl::fabs(d) < T(0.1) * (x + np)) {
+        T v = d / (x + np);
+        T s = d * v;
+        T ej = T(2) * x * v;
+        v = v * v;
+        for (int j = 1; j < 1000; ++j) {
+            ej *= v;
+            const T s1 = s + ej / static_cast<T>(2 * j + 1);
+            if (s1 == s) {
+                return s1;
+            }
+            s = s1;
+        }
+        return s;
+    }
+    return x * sycl::log(x / np) - d;
+}
+
+// log(a!) - log(b!) - (a - b)*log(a), for a and b of at least 16, with
+// log(n!) = stirlerr(n) + (n + 1/2)*log(n) - n + log(sqrt(2*pi)). Its terms
+// are small for b close to a, so a sum of such differences can gather the
+// large (a - b)*log(a) terms into one which does not cancel.
+template <typename T>
+T logfactorial_diff_rest(const std::int64_t a, const std::int64_t b)
+{
+    const T af = static_cast<T>(a);
+    const T bf = static_cast<T>(b);
+    const T d = static_cast<T>(a - b);
+    return stirlerr(af) - stirlerr(bf) - bd0(bf, af, -d) +
+           T(0.5) * sycl::log1p(d / bf);
+}
+
+// log(a!) - log(b!). The terms of the difference are of the order of
+// (a - b) * log(a), not of a * log(a), so the result keeps its precision
+// for large a and b, where lgamma(a + 1) - lgamma(b + 1) loses all of it.
+template <typename T>
+T logfactorial_diff(const std::int64_t a, const std::int64_t b)
+{
+    if (a == b) {
+        return T(0);
+    }
+    if (sycl::min(a, b) < 16) {
+        return sycl::lgamma(static_cast<T>(a) + T(1)) -
+               sycl::lgamma(static_cast<T>(b) + T(1));
+    }
+    return logfactorial_diff_rest<T>(a, b) +
+           static_cast<T>(a - b) * sycl::log(static_cast<T>(a));
+}
+
+// Poisson variates by multiplication of uniforms, for lam < 10
+template <typename T, typename BitGenT>
+std::int64_t poisson_mult(BitGenT &bg, const T lam)
+{
+    const T enlam = sycl::exp(-lam);
+    std::int64_t X = 0;
+    T prod = T(1);
+    for (;;) {
+        prod *= bg.template next_uniform<T>();
+        if (prod > enlam) {
+            X += 1;
+        }
+        else {
+            return X;
+        }
+    }
+}
+
+// Hormann's transformed rejection with squeeze (PTRS), for lam >= 10.
+//
+// Unlike NumPy's random_poisson_ptrs(), this keeps k as the integer part of
+// lam plus an offset, and evaluates the log-pmf in Loader's form, as with
+// floats both lam + offset and -lam + k*log(lam) - log(k!) lose too much.
+template <typename T, typename BitGenT>
+std::int64_t poisson_ptrs(BitGenT &bg, const T lam)
+{
+    constexpr T ln_2pi = T(1.83787706640934548356065947281);
+    // the largest T below 2**63, which lam may round up to in float
+    constexpr T lam_cap =
+        T(0x1p63) - T(0x1p63) * std::numeric_limits<T>::epsilon() / T(2);
+
+    const T slam = sycl::sqrt(lam);
+    const T b = T(0.931) + T(2.53) * slam;
+    const T a = T(-0.059) + T(0.02483) * b;
+    const T invalpha = T(1.1239) + T(1.1328) / (b - T(3.4));
+    const T vr = T(0.9277) - T(3.6224) / (b - T(2));
+
+    const T lam_fl = sycl::fmin(sycl::floor(lam), lam_cap);
+    const std::int64_t lam_i = static_cast<std::int64_t>(lam_fl);
+    const T lam_f = lam - lam_fl;
+
+    for (;;) {
+        const T U = bg.template next_uniform<T>() - T(0.5);
+        const T V = bg.template next_uniform<T>();
+        const T us = T(0.5) - sycl::fabs(U);
+        const T dk = sycl::floor((T(2) * a / us + b) * U + lam_f + T(0.43));
+
+        // k = lam_i + dk must be a non-negative int64; NumPy casts first,
+        // which is undefined out of range. This also rejects us == 0.
+        if (!(sycl::fabs(dk) < T(0x1p62))) {
+            continue;
+        }
+        const std::int64_t dk_i = static_cast<std::int64_t>(dk);
+        if ((dk_i < -lam_i) ||
+            (dk_i > std::numeric_limits<std::int64_t>::max() - lam_i)) {
+            continue;
+        }
+        const std::int64_t k = lam_i + dk_i;
+
+        if ((us >= T(0.07)) && (V <= vr)) {
+            return k;
+        }
+        if ((us < T(0.013)) && (V > us)) {
+            continue;
+        }
+
+        // log(lam**k * exp(-lam) / k!)
+        T log_pmf;
+        if (k == 0) {
+            log_pmf = -lam;
+        }
+        else {
+            const T kf = static_cast<T>(k);
+            log_pmf = -stirlerr(kf) - bd0(kf, lam, static_cast<T>(dk) - lam_f) -
+                      T(0.5) * (ln_2pi + sycl::log(kf));
+        }
+
+        // log(V) == log(0.0) is fine here
+        if ((sycl::log(V) + sycl::log(invalpha) -
+             sycl::log(a / (us * us) + b)) <= log_pmf) {
+            return k;
+        }
+    }
+}
+
+template <typename T, typename BitGenT>
+std::int64_t poisson(BitGenT &bg, const T lam)
+{
+    if (lam >= T(0x1p63)) {
+        // invalid input, which PTRS would reject forever; this is the value
+        // of the out of range cast in NumPy's legacy code. NaN gives 0 below,
+        // as there.
+        return std::numeric_limits<std::int64_t>::min();
+    }
+    if (lam >= T(10)) {
+        return poisson_ptrs<T>(bg, lam);
+    }
+    else if (lam == T(0)) {
+        return 0;
+    }
+    return poisson_mult<T>(bg, lam);
+}
+
+template <typename T, typename BitGenT>
+T noncentral_chisquare(BitGenT &bg,
+                       StandardNormal<T> &normal,
+                       const T df,
+                       const T nonc)
+{
+    if (sycl::isnan(nonc)) {
+        return std::numeric_limits<T>::quiet_NaN();
+    }
+    if (nonc == T(0)) {
+        return chisquare<T>(bg, normal, df);
+    }
+    if (T(1) < df) {
+        const T Chi2 = chisquare<T>(bg, normal, df - T(1));
+        const T n = normal(bg) + sycl::sqrt(nonc);
+        return Chi2 + n * n;
+    }
+    const std::int64_t i = poisson<T>(bg, nonc / T(2));
+    return chisquare<T>(bg, normal, df + T(2) * static_cast<T>(i));
+}
+
+template <typename T, typename BitGenT>
+T noncentral_f(BitGenT &bg,
+               StandardNormal<T> &normal,
+               const T dfnum,
+               const T dfden,
+               const T nonc)
+{
+    const T t = noncentral_chisquare<T>(bg, normal, dfnum, nonc) * dfden;
+    return t / (chisquare<T>(bg, normal, dfden) * dfnum);
+}
+
 template <typename T, typename BitGenT>
 T laplace(BitGenT &bg, const T loc, const T scale)
 {
@@ -408,4 +682,526 @@ T vonmises(BitGenT &bg, StandardNormal<T> &normal, const T mu, const T kappa)
     mod = sycl::fmod(mod + pi, T(2) * pi) - pi;
     return neg ? -mod : mod;
 }
+
+template <typename T, typename BitGenT>
+std::int64_t binomial_inversion(BitGenT &bg, const std::int64_t n, const T p)
+{
+    const T q = T(1) - p;
+    const T qn = sycl::exp(static_cast<T>(n) * sycl::log(q));
+    const T np = static_cast<T>(n) * p;
+    const T bound_f =
+        sycl::fmin(static_cast<T>(n), np + T(10) * sycl::sqrt(np * q + T(1)));
+    const std::int64_t bound = static_cast<std::int64_t>(bound_f);
+
+    std::int64_t X = 0;
+    T px = qn;
+    T U = bg.template next_uniform<T>();
+    while (U > px) {
+        X++;
+        if (X > bound) {
+            X = 0;
+            px = qn;
+            U = bg.template next_uniform<T>();
+        }
+        else {
+            U -= px;
+            px = (static_cast<T>(n - X + 1) * p * px) / (static_cast<T>(X) * q);
+        }
+    }
+    return X;
+}
+
+// floor((n + 1) * r) for 0 < r < 1, with its fractional part in `frac`,
+// exactly: r = M * 2**-k, and the 128-bit product (n + 1) * M is shifted
+template <typename T>
+std::int64_t binomial_mode(const std::int64_t n, const T r, T &frac)
+{
+    using UIntT = std::conditional_t<std::is_same_v<T, double>, std::uint64_t,
+                                     std::uint32_t>;
+    constexpr int mant_bits = std::numeric_limits<T>::digits - 1;
+    constexpr int bias = std::numeric_limits<T>::max_exponent - 1;
+
+    const UIntT bits = sycl::bit_cast<UIntT>(r);
+    const int e = static_cast<int>(bits >> mant_bits);
+    std::uint64_t M =
+        static_cast<std::uint64_t>(bits & ((UIntT(1) << mant_bits) - 1));
+    if (e != 0) {
+        M |= std::uint64_t(1) << mant_bits;
+    }
+    // k >= mant_bits + 1, as r < 1
+    const int k = bias + mant_bits - ((e != 0) ? e : 1);
+
+    const std::uint64_t n1 = static_cast<std::uint64_t>(n) + 1;
+    const std::uint64_t lo = n1 * M;
+    const std::uint64_t hi = umul_hi(n1, M);
+
+    std::uint64_t m;
+    std::uint64_t rem_hi = 0;
+    std::uint64_t rem_lo = lo;
+    if (k < 64) {
+        m = (hi << (64 - k)) | (lo >> k);
+        rem_lo = lo & ((std::uint64_t(1) << k) - 1);
+    }
+    else if (k < 128) {
+        m = hi >> (k - 64);
+        rem_hi = hi & ((std::uint64_t(1) << (k - 64)) - 1);
+    }
+    else {
+        m = 0;
+        rem_hi = hi;
+    }
+    frac = sycl::ldexp(static_cast<T>(rem_hi), 64 - k) +
+           sycl::ldexp(static_cast<T>(rem_lo), -k);
+    return static_cast<std::int64_t>(m);
+}
+
+// Kachitvichyanukul and Schmeiser's BTPE algorithm, for p <= 0.5. The steps
+// are numbered as in NumPy's random_binomial_btpe().
+//
+// Unlike NumPy's, the mode m is computed exactly in integers, every draw as
+// m plus an offset computed in floats, and the final acceptance test from
+// log-factorials whose large terms are gathered into one. Otherwise, above
+// 2**24 in single precision, or 2**53 in double, the draws would only take
+// the values of floats, and the test would cancel to noise.
+template <typename T, typename BitGenT>
+std::int64_t binomial_btpe(BitGenT &bg, const std::int64_t n, const T p)
+{
+    const T nf = static_cast<T>(n);
+    const T r = sycl::fmin(p, T(1) - p);
+    const T q = T(1) - r;
+    // fm = (n + 1) * r = m + g
+    T g;
+    const std::int64_t m = binomial_mode<T>(n, r, g);
+    const T mf = static_cast<T>(m);
+    const T nmf = static_cast<T>(n - m);
+    const T p1 =
+        sycl::floor(T(2.195) * sycl::sqrt(nf * r * q) - T(4.6) * q) + T(0.5);
+    // xm = m + 1/2, xl = xm - p1 and xr = xm + p1, all relative to m
+    const T xm = T(0.5);
+    const T xl = xm - p1;
+    const T xr = xm + p1;
+    const T c = T(0.134) + T(20.5) / (T(15.3) + mf);
+    // fm - xl = p1 + dm, xr - fm = p1 - dm and fm - xl * r = (n + 1 - xl) * r
+    // in absolute terms
+    const T dm = g - xm;
+    T a = (p1 + dm) / ((nmf + T(1) - xl) * r);
+    const T laml = a * (T(1) + a / T(2));
+    a = (p1 - dm) / ((mf + xr) * q);
+    const T lamr = a * (T(1) + a / T(2));
+    const T p2 = p1 * (T(1) + T(2) * c);
+    const T p3 = p2 + c / laml;
+    const T p4 = p3 + c / lamr;
+    const T nrq = nf * r * q;
+
+    for (;;) {
+        // Step 10
+        const T u = bg.template next_uniform<T>() * p4;
+        T v = bg.template next_uniform<T>();
+        std::int64_t y;
+        if (u <= p1) {
+            return m + static_cast<std::int64_t>(sycl::floor(xm - p1 * v + u));
+        }
+        else if (u <= p2) {
+            // Step 20
+            const T x = xl + (u - p1) / c;
+            v = v * c + T(1) - sycl::fabs(xm - x) / p1;
+            if (v > T(1)) {
+                continue;
+            }
+            y = m + static_cast<std::int64_t>(sycl::floor(x));
+        }
+        else if (u <= p3) {
+            // Step 30; reject v == 0 before the undefined cast
+            if (v == T(0)) {
+                continue;
+            }
+            const T dy = sycl::floor(xl + sycl::log(v) / laml);
+            if (!(dy >= -mf - T(1))) {
+                continue;
+            }
+            y = m + static_cast<std::int64_t>(dy);
+            if (y < 0) {
+                continue;
+            }
+            v = v * (u - p2) * laml;
+        }
+        else {
+            // Step 40; as above
+            if (v == T(0)) {
+                continue;
+            }
+            const T dy = sycl::floor(xr - sycl::log(v) / lamr);
+            if (!(dy <= nmf + T(1))) {
+                continue;
+            }
+            y = m + static_cast<std::int64_t>(dy);
+            if (y > n) {
+                continue;
+            }
+            v = v * (u - p3) * lamr;
+        }
+
+        // Step 50
+        const std::int64_t k = (y > m) ? y - m : m - y;
+        const T kf = static_cast<T>(k);
+        if (!((k > 20) && (kf < nrq / T(2) - T(1)))) {
+            const T s = r / q;
+            const T aa = s * (nf + T(1));
+            T F = T(1);
+            if (m < y) {
+                for (std::int64_t i = m + 1; i <= y; ++i) {
+                    F *= (aa / static_cast<T>(i) - s);
+                }
+            }
+            else if (m > y) {
+                for (std::int64_t i = y + 1; i <= m; ++i) {
+                    F /= (aa / static_cast<T>(i) - s);
+                }
+            }
+            if (v > F) {
+                continue;
+            }
+            return y;
+        }
+
+        // Step 52
+        const T rho =
+            (kf / nrq) *
+            ((kf * (kf / T(3) + T(0.625)) + T(0.16666666666666666)) / nrq +
+             T(0.5));
+        const T t = -kf * kf / (T(2) * nrq);
+        // log(0.0) is fine here
+        const T A = sycl::log(v);
+        if (A < (t - rho)) {
+            return y;
+        }
+        if (A > (t + rho)) {
+            continue;
+        }
+
+        // log(f(y) / f(m)) = log(m! (n - m)! / (y! (n - y)!)) + d*log(r / q)
+        // with d = y - m, whose large terms make d*log((n - m)*r / (m*q)),
+        // where (n - m)*r - m*q = g - r
+        const T d = static_cast<T>(y - m);
+        T log_ratio;
+        if (sycl::min(sycl::min(m, y), sycl::min(n - m, n - y)) >= 16) {
+            log_ratio = logfactorial_diff_rest<T>(m, y) +
+                        logfactorial_diff_rest<T>(n - m, n - y) +
+                        d * sycl::log1p((g - r) / (mf * q));
+        }
+        else {
+            log_ratio = logfactorial_diff<T>(m, y) +
+                        logfactorial_diff<T>(n - m, n - y) +
+                        d * sycl::log(r / q);
+        }
+        if (A > log_ratio) {
+            continue;
+        }
+        return y;
+    }
+}
+
+template <typename T, typename BitGenT>
+std::int64_t binomial(BitGenT &bg, const std::int64_t n, const T p)
+{
+    // Invalid params are only possible on the device, where they are not
+    // checked. BTPE rejects forever for a NaN p.
+    if (!(n >= 0) || !(p >= T(0) && p <= T(1))) {
+        return std::numeric_limits<std::int64_t>::min();
+    }
+    if ((n == 0) || (p == T(0))) {
+        return 0;
+    }
+
+    const T nf = static_cast<T>(n);
+    if (p <= T(0.5)) {
+        if (p * nf <= T(30)) {
+            return binomial_inversion<T>(bg, n, p);
+        }
+        return binomial_btpe<T>(bg, n, p);
+    }
+
+    const T q = T(1) - p;
+    if (q * nf <= T(30)) {
+        return n - binomial_inversion<T>(bg, n, q);
+    }
+    return n - binomial_btpe<T>(bg, n, q);
+}
+
+template <typename T, typename BitGenT>
+std::int64_t negative_binomial(BitGenT &bg,
+                               StandardNormal<T> &normal,
+                               const T n,
+                               const T p)
+{
+    const T Y = standard_gamma<T>(bg, normal, n) * ((T(1) - p) / p);
+    return poisson<T>(bg, Y);
+}
+
+template <typename T, typename BitGenT>
+std::int64_t geometric(BitGenT &bg, const T p)
+{
+    if (!(p > T(0))) {
+        // invalid, and the inversion would convert NaN or -inf to int64
+        return std::numeric_limits<std::int64_t>::min();
+    }
+    if (p >= T(0.333333333333333333333333)) {
+        // search
+        std::int64_t X = 1;
+        T sum = p;
+        T prod = p;
+        const T q = T(1) - p;
+        const T U = bg.template next_uniform<T>();
+        // rounding can stall the sum below U, until prod underflows
+        while ((U > sum) && (prod > T(0))) {
+            prod *= q;
+            sum += prod;
+            X++;
+        }
+        return X;
+    }
+
+    // inversion
+    const T z = sycl::ceil(-standard_exponential<T>(bg) / sycl::log1p(-p));
+    // the smallest value above INT64_MAX in both float and double
+    if (z >= T(9.223372036854775808e18)) {
+        return std::numeric_limits<std::int64_t>::max();
+    }
+    return static_cast<std::int64_t>(z);
+}
+
+template <typename T, typename BitGenT>
+std::int64_t logseries(BitGenT &bg, const T p)
+{
+    if (!(p < T(1))) {
+        // invalid input; for p == 1 every draw would be rejected
+        return std::numeric_limits<std::int64_t>::min();
+    }
+    const T r = sycl::log1p(-p);
+    for (;;) {
+        const T V = bg.template next_uniform<T>();
+        if (V >= p) {
+            return 1;
+        }
+        const T U = bg.template next_uniform<T>();
+        const T q = -sycl::expm1(r * U);
+        if (V <= (q * q)) {
+            if (V == T(0)) {
+                continue;
+            }
+            const T result = sycl::floor(T(1) + sycl::log(V) / sycl::log(q));
+            if (!(result >= T(1) && result < T(9.223372036854775808e18))) {
+                continue;
+            }
+            return static_cast<std::int64_t>(result);
+        }
+        if (V >= q) {
+            return 1;
+        }
+        return 2;
+    }
+}
+
+template <typename T, typename BitGenT>
+std::int64_t zipf(BitGenT &bg, const T a)
+{
+    if (!(a > T(1))) {
+        // invalid, and every draw would be rejected
+        return 0;
+    }
+    if (a >= T(std::numeric_limits<T>::max_exponent + 1)) {
+        // every draw is 1, and 2**(a - 1) would overflow, so that every
+        // draw would be rejected
+        return 1;
+    }
+
+    // the smallest value above INT64_MAX in both float and double
+    constexpr T int64_bound = T(9.223372036854775808e18);
+
+    const T am1 = a - T(1);
+    const T b = sycl::pow(T(2), am1);
+    // reject values that do not fit an int64, so this models a Zipf
+    // distribution truncated to INT64_MAX
+    const T Umin = sycl::pow(int64_bound, -am1);
+    for (;;) {
+        // U is on (Umin, 1], and Umin may be zero
+        const T U01 = bg.template next_uniform<T>();
+        const T U = U01 * Umin + (T(1) - U01);
+        const T V = bg.template next_uniform<T>();
+        const T X = sycl::floor(sycl::pow(U, T(-1) / am1));
+        if (!(X >= T(1) && X < int64_bound)) {
+            continue;
+        }
+        // t = (1 + 1/X)**(a - 1), whose t - 1 would round to 0 for large X
+        const T tm1 = sycl::expm1(am1 * sycl::log1p(T(1) / X));
+        const T t = T(1) + tm1;
+        if (V * X * tm1 / (b - T(1)) <= t / b) {
+            return static_cast<std::int64_t>(X);
+        }
+    }
+}
+
+// Draw `sample` items without replacement one by one, for small samples
+template <typename T, typename BitGenT>
+std::int64_t hypergeometric_sample(BitGenT &bg,
+                                   const std::int64_t good,
+                                   const std::int64_t bad,
+                                   const std::int64_t sample)
+{
+    const std::int64_t total = good + bad;
+    std::int64_t computed_sample =
+        (sample > total / 2) ? total - sample : sample;
+
+    std::int64_t remaining_total = total;
+    std::int64_t remaining_good = good;
+    while ((computed_sample > 0) && (remaining_good > 0) &&
+           (remaining_total > remaining_good)) {
+        --remaining_total;
+        if (static_cast<std::int64_t>(bounded_uint64(
+                bg, static_cast<std::uint64_t>(remaining_total))) <
+            remaining_good) {
+            --remaining_good;
+        }
+        --computed_sample;
+    }
+
+    if (remaining_total == remaining_good) {
+        // only good items are left
+        remaining_good -= computed_sample;
+    }
+
+    return (sample > total / 2) ? remaining_good : good - remaining_good;
+}
+
+// Stadlober's ratio-of-uniforms method (HRUA).
+//
+// Unlike NumPy's hypergeometric_hrua(), the mode m and the candidate K are
+// computed in integers, as the integer part of the mean plus an offset, and
+// g - gp as a sum of stable differences of log-factorials, so that neither
+// loses its precision with floats and large populations. NumPy limits good
+// and bad to below 10**9, so the integer products cannot overflow.
+template <typename T, typename BitGenT>
+std::int64_t hypergeometric_hrua(BitGenT &bg,
+                                 const std::int64_t good,
+                                 const std::int64_t bad,
+                                 const std::int64_t sample)
+{
+    // D1 = 2*sqrt(2/e), D2 = 3 - 2*sqrt(3/e)
+    constexpr T D1 = T(1.7155277699214135);
+    constexpr T D2 = T(0.8989161620588988);
+
+    const std::int64_t popsize = good + bad;
+    const std::int64_t computed_sample = sycl::min(sample, popsize - sample);
+    const std::int64_t mingoodbad = sycl::min(good, bad);
+    const std::int64_t maxgoodbad = sycl::max(good, bad);
+
+    const T p = static_cast<T>(mingoodbad) / static_cast<T>(popsize);
+    const T q = static_cast<T>(maxgoodbad) / static_cast<T>(popsize);
+
+    // the mean mu = computed_sample * p, and a = mu + 1/2 = a_i + a_f
+    const std::int64_t mu_num = computed_sample * mingoodbad;
+    const std::int64_t a_i = mu_num / popsize;
+    const T a_f =
+        static_cast<T>(mu_num % popsize) / static_cast<T>(popsize) + T(0.5);
+
+    const T var = static_cast<T>(popsize - computed_sample) *
+                  static_cast<T>(computed_sample) * p * q /
+                  static_cast<T>(popsize - 1);
+    const T c = sycl::sqrt(var + T(0.5));
+    const T h = D1 * c + D2;
+
+    const std::int64_t m =
+        (computed_sample + 1) * (mingoodbad + 1) / (popsize + 2);
+    const std::int64_t b = sycl::min(
+        sycl::min(computed_sample, mingoodbad) + 1,
+        a_i + static_cast<std::int64_t>(sycl::floor(a_f + T(16) * c)));
+
+    std::int64_t K;
+    for (;;) {
+        const T U = bg.template next_uniform<T>();
+        const T V = bg.template next_uniform<T>();
+        // X = a + h*(V - 1/2)/U = a_i + dX
+        const T dX = a_f + h * (V - T(0.5)) / U;
+
+        // fast rejection of X < 0 or X >= b; this also rejects U == 0
+        if (!(dX >= -static_cast<T>(a_i) - T(1) &&
+              dX < static_cast<T>(b - a_i) + T(1))) {
+            continue;
+        }
+        K = a_i + static_cast<std::int64_t>(sycl::floor(dX));
+        if ((K < 0) || (K >= b)) {
+            continue;
+        }
+
+        // g - gp, where g and gp are sums of four log-factorials at m and K
+        const std::int64_t a1 = m;
+        const std::int64_t a2 = mingoodbad - m;
+        const std::int64_t a3 = computed_sample - m;
+        const std::int64_t a4 = maxgoodbad - computed_sample + m;
+        const std::int64_t b1 = K;
+        const std::int64_t b2 = mingoodbad - K;
+        const std::int64_t b3 = computed_sample - K;
+        const std::int64_t b4 = maxgoodbad - computed_sample + K;
+        T t;
+        if (sycl::min(sycl::min(sycl::min(a1, a2), sycl::min(a3, a4)),
+                      sycl::min(sycl::min(b1, b2), sycl::min(b3, b4))) >= 16) {
+            // The terms (m - K)*log(a1 * a4 / (a2 * a3)) of the differences
+            // are each about |m - K|*log(popsize), while their sum is small,
+            // as m is the mode; gather them with the exact difference of the
+            // integer products, below 2 * 10**18 for good and bad below 10**9
+            const std::int64_t num = a1 * a4;
+            const std::int64_t den = a2 * a3;
+            t = logfactorial_diff_rest<T>(a1, b1) +
+                logfactorial_diff_rest<T>(a2, b2) +
+                logfactorial_diff_rest<T>(a3, b3) +
+                logfactorial_diff_rest<T>(a4, b4) +
+                static_cast<T>(m - K) * sycl::log1p(static_cast<T>(num - den) /
+                                                    static_cast<T>(den));
+        }
+        else {
+            t = logfactorial_diff<T>(a1, b1) + logfactorial_diff<T>(a2, b2) +
+                logfactorial_diff<T>(a3, b3) + logfactorial_diff<T>(a4, b4);
+        }
+
+        // fast acceptance
+        if ((U * (T(4) - U) - T(3)) <= t) {
+            break;
+        }
+        // fast rejection
+        if (U * (U - t) >= T(1)) {
+            continue;
+        }
+        if (T(2) * sycl::log(U) <= t) {
+            break;
+        }
+    }
+
+    if (good > bad) {
+        K = computed_sample - K;
+    }
+    if (computed_sample < sample) {
+        K = good - K;
+    }
+    return K;
+}
+
+template <typename T, typename BitGenT>
+std::int64_t hypergeometric(BitGenT &bg,
+                            const std::int64_t good,
+                            const std::int64_t bad,
+                            const std::int64_t sample)
+{
+    // Invalid params are only possible on the device, where they are not
+    // checked. HRUA's integer products overflow for good or bad of 10**9 or
+    // more, after which it may reject forever.
+    constexpr std::int64_t max_good_bad = 1000000000;
+    if (!((good >= 0) && (good < max_good_bad) && (bad >= 0) &&
+          (bad < max_good_bad) && (sample >= 0) && (sample <= good + bad))) {
+        return std::numeric_limits<std::int64_t>::min();
+    }
+    if ((sample >= 10) && (sample <= good + bad - 10)) {
+        return hypergeometric_hrua<T>(bg, good, bad, sample);
+    }
+    return hypergeometric_sample<T>(bg, good, bad, sample);
+}
 } // namespace dpnp::kernels::rng
diff --git a/dpnp/libtensor/extensions/rng/kernels/samplers.hpp b/dpnp/libtensor/extensions/rng/kernels/samplers.hpp
index c63cf4f2b8d..7e772e5e1dc 100644
--- a/dpnp/libtensor/extensions/rng/kernels/samplers.hpp
+++ b/dpnp/libtensor/extensions/rng/kernels/samplers.hpp
@@ -29,7 +29,9 @@
 #pragma once
 
 #include <cstddef>
+#include <cstdint>
 #include <tuple>
+#include <type_traits>
 
 #include <sycl/sycl.hpp>
 
@@ -107,6 +109,37 @@ class UniformSampler : public SamplerBase<T, T, T>
     }
 };
 
+// Uniform integers on the closed interval [low, high]
+template <typename T>
+class IntegersSampler : public SamplerBase<T, T, T>
+{
+private:
+    Param<T> low_;
+    Param<T> high_;
+
+public:
+    IntegersSampler(const Param<T> &low, const Param<T> &high)
+        : low_(low), high_(high)
+    {
+    }
+
+    template <typename BitGenT>
+    T operator()(BitGenT &bg, const std::size_t i)
+    {
+        // Widen to 64 bits with the sign of T, then do the arithmetic
+        // modulo 2**64, which is exact for any low <= high.
+        using wideT = std::conditional_t<std::is_signed_v<T>, std::int64_t,
+                                         std::uint64_t>;
+
+        const std::uint64_t off =
+            static_cast<std::uint64_t>(static_cast<wideT>(low_[i]));
+        const std::uint64_t rng =
+            static_cast<std::uint64_t>(static_cast<wideT>(high_[i])) - off;
+
+        return static_cast<T>(off + bounded_uint64(bg, rng));
+    }
+};
+
 // Normal with mean `loc` and standard deviation `scale`
 template <typename T>
 class NormalSampler : public SamplerBase<T, T, T>
@@ -404,6 +437,54 @@ class FSampler : public SamplerBase<T, T, T>
     }
 };
 
+// Noncentral chi-square with `df` degrees of freedom and noncentrality `nonc`
+template <typename T>
+class NoncentralChisquareSampler : public SamplerBase<T, T, T>
+{
+private:
+    Param<T> df_;
+    Param<T> nonc_;
+    StandardNormal<T> normal_{};
+
+public:
+    NoncentralChisquareSampler(const Param<T> &df, const Param<T> &nonc)
+        : df_(df), nonc_(nonc)
+    {
+    }
+
+    template <typename BitGenT>
+    T operator()(BitGenT &bg, const std::size_t i)
+    {
+        return noncentral_chisquare<T>(bg, normal_, df_[i], nonc_[i]);
+    }
+};
+
+// Noncentral F with `dfnum` and `dfden` degrees of freedom and noncentrality
+// `nonc`
+template <typename T>
+class NoncentralFSampler : public SamplerBase<T, T, T, T>
+{
+private:
+    Param<T> dfnum_;
+    Param<T> dfden_;
+    Param<T> nonc_;
+    StandardNormal<T> normal_{};
+
+public:
+    NoncentralFSampler(const Param<T> &dfnum,
+                       const Param<T> &dfden,
+                       const Param<T> &nonc)
+        : dfnum_(dfnum), dfden_(dfden), nonc_(nonc)
+    {
+    }
+
+    template <typename BitGenT>
+    T operator()(BitGenT &bg, const std::size_t i)
+    {
+        return noncentral_f<T>(bg, normal_, dfnum_[i], dfden_[i], nonc_[i]);
+    }
+};
+
 // Laplace with location `loc` and scale `scale`
 template <typename T>
 class LaplaceSampler : public SamplerBase<T, T, T>
@@ -515,4 +596,143 @@ class VonmisesSampler : public SamplerBase<T, T, T>
         return vonmises<T>(bg, normal_, mu_[i], kappa_[i]);
     }
 };
+
+// The discrete samplers below give int64 samples, and compute in the
+// floating-point type T.
+
+// Poisson with rate `lam`
+template <typename T>
+class PoissonSampler : public SamplerBase<std::int64_t, T>
+{
+private:
+    Param<T> lam_;
+
+public:
+    explicit PoissonSampler(const Param<T> &lam) : lam_(lam) {}
+
+    template <typename BitGenT>
+    std::int64_t operator()(BitGenT &bg, const std::size_t i)
+    {
+        return poisson<T>(bg, lam_[i]);
+    }
+};
+
+// The number of failures before `n` successes of probability `p`
+template <typename T>
+class NegativeBinomialSampler : public SamplerBase<std::int64_t, T, T>
+{
+private:
+    Param<T> n_;
+    Param<T> p_;
+    StandardNormal<T> normal_{};
+
+public:
+    NegativeBinomialSampler(const Param<T> &n, const Param<T> &p) : n_(n), p_(p)
+    {
+    }
+
+    template <typename BitGenT>
+    std::int64_t operator()(BitGenT &bg, const std::size_t i)
+    {
+        return negative_binomial<T>(bg, normal_, n_[i], p_[i]);
+    }
+};
+
+// The number of successes of `n` trials of probability `p`
+template <typename T>
+class BinomialSampler : public SamplerBase<std::int64_t, std::int64_t, T>
+{
+private:
+    Param<std::int64_t> n_;
+    Param<T> p_;
+
+public:
+    BinomialSampler(const Param<std::int64_t> &n, const Param<T> &p)
+        : n_(n), p_(p)
+    {
+    }
+
+    template <typename BitGenT>
+    std::int64_t operator()(BitGenT &bg, const std::size_t i)
+    {
+        return binomial<T>(bg, n_[i], p_[i]);
+    }
+};
+
+// The number of trials of probability `p` up to the first success
+template <typename T>
+class GeometricSampler : public SamplerBase<std::int64_t, T>
+{
+private:
+    Param<T> p_;
+
+public:
+    explicit GeometricSampler(const Param<T> &p) : p_(p) {}
+
+    template <typename BitGenT>
+    std::int64_t operator()(BitGenT &bg, const std::size_t i)
+    {
+        return geometric<T>(bg, p_[i]);
+    }
+};
+
+// Logarithmic series with shape `p`
+template <typename T>
+class LogseriesSampler : public SamplerBase<std::int64_t, T>
+{
+private:
+    Param<T> p_;
+
+public:
+    explicit LogseriesSampler(const Param<T> &p) : p_(p) {}
+
+    template <typename BitGenT>
+    std::int64_t operator()(BitGenT &bg, const std::size_t i)
+    {
+        return logseries<T>(bg, p_[i]);
+    }
+};
+
+// Zipf with exponent `a`
+template <typename T>
+class ZipfSampler : public SamplerBase<std::int64_t, T>
+{
+private:
+    Param<T> a_;
+
+public:
+    explicit ZipfSampler(const Param<T> &a) : a_(a) {}
+
+    template <typename BitGenT>
+    std::int64_t operator()(BitGenT &bg, const std::size_t i)
+    {
+        return zipf<T>(bg, a_[i]);
+    }
+};
+
+// The number of good items of `sample` drawn without replacement from `good`
+// good and `bad` bad items
+template <typename T>
+class HypergeometricSampler
+    : public SamplerBase<std::int64_t, std::int64_t, std::int64_t, std::int64_t>
+{
+private:
+    Param<std::int64_t> good_;
+    Param<std::int64_t> bad_;
+    Param<std::int64_t> sample_;
+
+public:
+    HypergeometricSampler(const Param<std::int64_t> &good,
+                          const Param<std::int64_t> &bad,
+                          const Param<std::int64_t> &sample)
+        : good_(good), bad_(bad), sample_(sample)
+    {
+    }
+
+    template <typename BitGenT>
+    std::int64_t operator()(BitGenT &bg, const std::size_t i)
+    {
+        return hypergeometric<T>(bg, good_[i], bad_[i], sample_[i]);
+    }
+};
 } // namespace dpnp::kernels::rng
diff --git a/dpnp/libtensor/extensions/rng/kernels/umul_hi.hpp b/dpnp/libtensor/extensions/rng/kernels/umul_hi.hpp
new file mode 100644
index 00000000000..5399077ec82
--- /dev/null
+++ b/dpnp/libtensor/extensions/rng/kernels/umul_hi.hpp
@@ -0,0 +1,55 @@
+//*****************************************************************************
+// Copyright (c) 2026, Intel Corporation
+// All rights reserved.
+//
+// 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 the copyright holder 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 HOLDER 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.
+//*****************************************************************************
+
+#pragma once
+
+#include <cstdint>
+
+namespace dpnp::kernels::rng
+{
+// The high 64 bits of the 128-bit product a * b, from 32-bit halves.
+// sycl::mul_hi is not used: optimized DPC++ 2026.2 builds get it wrong on
+// Intel GPUs when one 64-bit operand is a compile-time constant (all results
+// wrong at -O2, right at -O0). oneMKL's device RNG uses 32-bit halves too.
+inline std::uint64_t umul_hi(const std::uint64_t a, const std::uint64_t b)
+{
+    const std::uint64_t a_lo = a & 0xffffffffULL;
+    const std::uint64_t a_hi = a >> 32;
+    const std::uint64_t b_lo = b & 0xffffffffULL;
+    const std::uint64_t b_hi = b >> 32;
+
+    const std::uint64_t lo_lo = a_lo * b_lo;
+    const std::uint64_t hi_lo = a_hi * b_lo;
+    const std::uint64_t lo_hi = a_lo * b_hi;
+
+    // the carry out of the low 64 bits
+    const std::uint64_t mid =
+        (lo_lo >> 32) + (hi_lo & 0xffffffffULL) + (lo_hi & 0xffffffffULL);
+    return a_hi * b_hi + (hi_lo >> 32) + (lo_hi >> 32) + (mid >> 32);
+}
+} // namespace dpnp::kernels::rng
diff --git a/dpnp/libtensor/extensions/rng/noncentral.cpp b/dpnp/libtensor/extensions/rng/noncentral.cpp
new file mode 100644
index 00000000000..991f21efc73
--- /dev/null
+++ b/dpnp/libtensor/extensions/rng/noncentral.cpp
@@ -0,0 +1,52 @@
+//*****************************************************************************
+// Copyright (c) 2026, Intel Corporation
+// All rights reserved.
+//
+// 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 the copyright holder 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 HOLDER 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.
+//*****************************************************************************
+#include <pybind11/pybind11.h>
+
+#include "common.hpp"
+#include "distributions.hpp"
+
+namespace dpnp::extensions::rng
+{
+static sample_dispatch_table_t noncentral_chisquare_dispatch_table;
+static sample_dispatch_table_t noncentral_f_dispatch_table;
+
+void init_noncentral(py::module_ m)
+{
+    def_sampler<rng_ns::NoncentralChisquareSampler, RealTypes>(
+        m, "_noncentral_chisquare",
+        "Fill the result with noncentral chi-square samples of df degrees "
+        "of freedom and noncentrality nonc",
+        noncentral_chisquare_dispatch_table);
+
+    def_sampler<rng_ns::NoncentralFSampler, RealTypes>(
+        m, "_noncentral_f",
+        "Fill the result with noncentral F samples of dfnum and dfden "
+        "degrees of freedom and noncentrality nonc",
+        noncentral_f_dispatch_table);
+}
+} // namespace dpnp::extensions::rng
diff --git a/dpnp/libtensor/extensions/rng/poisson.cpp b/dpnp/libtensor/extensions/rng/poisson.cpp
new file mode 100644
index 00000000000..5d6ba347069
--- /dev/null
+++ b/dpnp/libtensor/extensions/rng/poisson.cpp
@@ -0,0 +1,50 @@
+//*****************************************************************************
+// Copyright (c) 2026, Intel Corporation
+// All rights reserved.
+//
+// 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 the copyright holder 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 HOLDER 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.
+//*****************************************************************************
+#include <pybind11/pybind11.h>
+
+#include "common.hpp"
+#include "distributions.hpp"
+
+namespace dpnp::extensions::rng
+{
+static discrete_dispatch_table_t poisson_dispatch_table;
+static discrete_dispatch_table_t negative_binomial_dispatch_table;
+
+void init_poisson(py::module_ m)
+{
+    def_discrete_sampler<rng_ns::PoissonSampler>(
+        m, "_poisson", "Fill the result with Poisson samples of rate lam",
+        poisson_dispatch_table);
+
+    def_discrete_sampler<rng_ns::NegativeBinomialSampler>(
+        m, "_negative_binomial",
+        "Fill the result with negative binomial samples of n successes of "
+        "probability p",
+        negative_binomial_dispatch_table);
+}
+} // namespace dpnp::extensions::rng
diff --git a/dpnp/libtensor/extensions/rng/rng_py.cpp b/dpnp/libtensor/extensions/rng/rng_py.cpp
index e7f058d9770..512243045ec 100644
--- a/dpnp/libtensor/extensions/rng/rng_py.cpp
+++ b/dpnp/libtensor/extensions/rng/rng_py.cpp
@@ -58,4 +58,9 @@ PYBIND11_MODULE(_rng_impl, m, py::mod_gil_not_used())
     rng_ns::init_normal(m);
     rng_ns::init_exponential(m);
     rng_ns::init_gamma(m);
+    rng_ns::init_noncentral(m);
+    rng_ns::init_integers(m);
+    rng_ns::init_poisson(m);
+    rng_ns::init_binomial(m);
+    rng_ns::init_discrete(m);
 }
diff --git a/dpnp/libtensor/extensions/rng/sample.cpp b/dpnp/libtensor/extensions/rng/sample.cpp
index 57cbbabcae8..2dc5524785d 100644
--- a/dpnp/libtensor/extensions/rng/sample.cpp
+++ b/dpnp/libtensor/extensions/rng/sample.cpp
@@ -105,6 +105,14 @@ sample_fn_ptr_t lookup_sampler(const sample_dispatch_table_t &table,
     return fn;
 }
 
+sample_fn_ptr_t lookup_sampler(const discrete_dispatch_table_t &table,
+                               const int engine,
+                               const bool double_precision)
+{
+    check_engine(engine);
+    return table[engine][double_precision ? 1 : 0];
+}
+
 std::pair<sycl::event, sycl::event>
     py_sample(sycl::queue &exec_q,
               const std::vector<std::uint64_t> &seed,
diff --git a/dpnp/random/_bit_generator.py b/dpnp/random/_bit_generator.py
index 9c38a48110c..ee90b2d6fb8 100644
--- a/dpnp/random/_bit_generator.py
+++ b/dpnp/random/_bit_generator.py
@@ -168,13 +168,14 @@ def _launch(self, kernel, dtype, size):
         """
         raise NotImplementedError
 
-    def _fill(self, kernel, params, result):
+    def _fill(self, kernel, params, result, **kwargs):
         """
         Fill the C-contiguous array `result` by `kernel` of the extension,
         with the distribution parameters `params`.
 
         Each parameter is either a Python scalar or a C-contiguous array of
         the size of `result` and of the data type the kernel expects for it.
+        Keyword arguments are passed to `kernel` as is.
 
         """
 
@@ -199,6 +200,7 @@ def _fill(self, kernel, params, result):
             params,
             result_usm,
             depends=_manager.submitted_events,
+            **kwargs,
         )
         _manager.add_event_pair(ht_ev, ev)
         return result
diff --git a/dpnp/random/_generator.py b/dpnp/random/_generator.py
index 0367d239ace..4274b1788d1 100644
--- a/dpnp/random/_generator.py
+++ b/dpnp/random/_generator.py
@@ -54,11 +54,44 @@
 __all__ = ["Generator", "default_rng"]
 
 _float_types = (numpy.dtype(numpy.float32), numpy.dtype(numpy.float64))
+_integer_types = tuple(
+    numpy.dtype(t)
+    for t in (
+        numpy.bool_,
+        numpy.int8,
+        numpy.uint8,
+        numpy.int16,
+        numpy.uint16,
+        numpy.int32,
+        numpy.uint32,
+        numpy.int64,
+        numpy.uint64,
+    )
+)
+_int64 = numpy.dtype(numpy.int64)
+_INT64_MAX = int(numpy.iinfo(numpy.int64).max)
+
 # The parameter constraints of NumPy's distributions, see
 # numpy/random/_common.pyx
 _CONS_NONE = 0
 _CONS_NON_NEGATIVE = 1
 _CONS_POSITIVE = 2
+_CONS_POSITIVE_NOT_NAN = 3
+_CONS_BOUNDED_0_1 = 4
+_CONS_BOUNDED_GT_0_1 = 5
+_CONS_BOUNDED_LT_0_1 = 6
+_CONS_GT_1 = 7
+_CONS_POISSON = 8
+
+_POISSON_LAM_MAX = _INT64_MAX - numpy.sqrt(_INT64_MAX) * 10
+_HYPERGEOM_MAX = 10**9
+
+# bounds keeping valid host parameters valid once rounded to float32; the
+# smallest normal also survives a device flushing subnormals to zero
+_F32_GT_0 = (numpy.finfo(numpy.float32).tiny, None)
+_F32_GT_1 = (numpy.nextafter(numpy.float32(1), numpy.float32(2)), None)
+_F32_LT_1 = (None, numpy.nextafter(numpy.float32(1), numpy.float32(0)))
+_F32_LT_2_63 = (None, numpy.nextafter(numpy.float32(2**63), numpy.float32(0)))
 
 
 def _is_array(x):
@@ -154,6 +187,7 @@ def _check_arrays(params, dtype=None):
     return numpy, [numpy.asarray(p, dtype) for p in params]
 
 
+# pylint: disable=too-many-branches
 def _check_constraint(value, name, cons, scalar):
     """
     Raise ValueError, as NumPy does, when the parameter `value` violates the
@@ -172,12 +206,64 @@ def _check_constraint(value, name, cons, scalar):
     if cons == _CONS_NONE:
         return
 
+    nan_msg = "is NaN" if scalar else "contains NaNs"
     if cons == _CONS_NON_NEGATIVE:
         if numpy.any(~numpy.isnan(val) & numpy.signbit(val)):
             raise ValueError(f"{name} < 0")
-    elif cons == _CONS_POSITIVE:
+    elif cons in (_CONS_POSITIVE, _CONS_POSITIVE_NOT_NAN):
+        if cons == _CONS_POSITIVE_NOT_NAN and numpy.any(numpy.isnan(val)):
+            raise ValueError(f"{name} must not be NaN")
         if numpy.any(val <= 0):
             raise ValueError(f"{name} <= 0")
+    elif cons == _CONS_BOUNDED_0_1:
+        if not numpy.all((val >= 0) & (val <= 1)):
+            raise ValueError(f"{name} < 0, {name} > 1 or {name} {nan_msg}")
+    elif cons == _CONS_BOUNDED_GT_0_1:
+        if not numpy.all((val > 0) & (val <= 1)):
+            raise ValueError(f"{name} <= 0, {name} > 1 or {name} contains NaNs")
+    elif cons == _CONS_BOUNDED_LT_0_1:
+        if not numpy.all((val >= 0) & (val < 1)):
+            raise ValueError(f"{name} < 0, {name} >= 1 or {name} {nan_msg}")
+    elif cons == _CONS_GT_1:
+        if not numpy.all(val > 1):
+            raise ValueError(f"{name} <= 1 or {name} {nan_msg}")
+    elif cons == _CONS_POISSON:
+        # the array path checks the upper bound first, so NaN fails it
+        if scalar:
+            if not val >= 0:
+                raise ValueError(f"{name} < 0 or {name} is NaN")
+            if val > _POISSON_LAM_MAX:
+                raise ValueError(f"{name} value too large")
+        else:
+            if not numpy.all(val <= _POISSON_LAM_MAX):
+                raise ValueError(f"{name} value too large")
+            if not numpy.all(val >= 0):
+                raise ValueError(f"{name} < 0 or {name} contains NaNs")
+
+
+def _host_float32(value, bounds):
+    """
+    Round a host parameter to ``float32`` and clip it to `bounds`, which
+    rounding a valid parameter to the nearest ``float32`` may cross.
+
+    """
+
+    if dpnp.is_supported_array_type(value):
+        return value
+    with numpy.errstate(over="ignore"):
+        value = numpy.asarray(value, dtype=numpy.float64).astype(numpy.float32)
+    value = numpy.clip(value, *bounds)
+    return value if value.ndim else float(value)
+
+
+def _as_int64(x):
+    """Truncate an integer parameter to int64, as NumPy does."""
+
+    if dpnp.is_supported_array_type(x):
+        return dpnp.asarray(x).astype(_int64)
+    if numpy.ndim(x) > 0:
+        return numpy.asarray(x).astype(_int64)
+    return int(x)
 
 
 # pylint: disable=too-many-public-methods
@@ -206,8 +292,14 @@ class Generator:
 
     Parameters on the host (scalars, lists and NumPy arrays) are checked as
     in NumPy. Parameters on the device are not, since that would wait for
-    the device. Invalid ones give undefined values, such as NaN, but never
-    make the kernels loop without end or access memory out of bounds.
+    the device. Invalid ones give undefined values, such as NaN or the
+    smallest ``int64``, but never make the kernels loop without end or
+    access memory out of bounds.
+
+    The discrete distributions return ``int64`` values, and compute in
+    double precision on devices supporting it, and in single precision
+    otherwise. On the latter, the samples for very large parameters, such
+    as a Poisson `lam` beyond about ``1e7``, are less accurate than NumPy's.
 
     """
 
@@ -291,6 +383,13 @@ def _float_dtype(self, dtype, name):
             raise TypeError(f"Unsupported dtype {dtype!r} for {name}")
         return dtype
 
+    def _discrete_float_type(self):
+        """The floating-point type the discrete distributions compute in."""
+
+        if self.sycl_device.has_aspect_fp64:
+            return numpy.dtype(numpy.float64)
+        return numpy.dtype(numpy.float32)
+
     def _prepare_out(self, out, size, dtype, name):
         """Resolve `size` and `dtype` against a user-supplied `out`."""
 
@@ -323,6 +422,8 @@ def _sample(
         *,
         out=None,
         kernel_params=None,
+        param_types=None,
+        **kwargs,
     ):
         """
         Draw samples by `kernel` of the extension into a new array, or into
@@ -330,8 +431,9 @@ def _sample(
 
         `params` are the user's distribution parameters, which set the output
         shape along with `size`. `kernel_params` are the parameters passed to
-        the kernel, `params` by default, of data type `dtype`. Every kernel
-        parameter that is an array is broadcast to the output shape.
+        the kernel, `params` by default, and `param_types` their data types,
+        `dtype` by default. Every kernel parameter that is an array is
+        broadcast to the output shape. `kwargs` are passed to the kernel.
 
         """
 
@@ -357,24 +459,32 @@ def _sample(
 
         if kernel_params is None:
             kernel_params = params
+        if param_types is None:
+            param_types = (dtype,) * len(kernel_params)
 
-        def _to_kernel_param(p):
+        def _to_kernel_param(p, dt):
             if _is_array(p):
                 if not dpnp.is_supported_array_type(p):
                     # NumPy converts None to NaN, dpnp.asarray does not
-                    p = numpy.asarray(p, dtype=dtype)
-                p = dpnp.asarray(p, dtype=dtype, sycl_queue=queue)
+                    p = numpy.asarray(p, dtype=dt)
+                p = dpnp.asarray(p, dtype=dt, sycl_queue=queue)
                 return dpnp.ascontiguousarray(dpnp.broadcast_to(p, shape))
-            # as in NumPy, None is NaN along with array parameters
-            return float(numpy.asarray(p, dtype=numpy.float64))
-
-        kernel_params = tuple(_to_kernel_param(p) for p in kernel_params)
+            if dt.kind == "f":
+                # as in NumPy, None is NaN along with array parameters
+                return float(numpy.asarray(p, dtype=numpy.float64))
+            if dt.kind == "b":
+                return bool(p)
+            return int(p)
+
+        kernel_params = tuple(
+            _to_kernel_param(p, dt) for p, dt in zip(kernel_params, param_types)
+        )
 
         if out is None:
             out = dpnp.empty(
                 shape, dtype=dtype, usm_type=usm_type, sycl_queue=queue
             )
-        return self._bit_generator._fill(kernel, kernel_params, out)
+        return self._bit_generator._fill(kernel, kernel_params, out, **kwargs)
 
     def _float_sample(
         self,
@@ -419,6 +529,50 @@ def _cont(self, kernel, name, size, dtype, usm_type, *params, out=None):
             kernel, values, size, dtype, usm_type, name, out=out
         )
 
+    def _discrete_sample(
+        self,
+        kernel,
+        params,
+        size,
+        usm_type,
+        *,
+        kernel_params=None,
+        param_types=None,
+        f32_bounds=None,
+    ):
+        """
+        Draw ``int64`` samples of a discrete distribution.
+
+        A ``None`` in `param_types` stands for the floating-point type the
+        distribution computes in, which is also the default for every
+        parameter. When that type is ``float32``, the host kernel parameters
+        are clipped to the matching `f32_bounds`, unless they are ``None``.
+
+        """
+
+        validate_usm_type(usm_type, allow_none=False)
+        ft = self._discrete_float_type()
+        if kernel_params is None:
+            kernel_params = params
+        if f32_bounds is not None and ft == numpy.float32:
+            kernel_params = tuple(
+                p if b is None else _host_float32(p, b)
+                for p, b in zip(kernel_params, f32_bounds)
+            )
+        if param_types is None:
+            param_types = (None,) * len(kernel_params)
+        param_types = tuple(ft if t is None else t for t in param_types)
+        return self._sample(
+            kernel,
+            params,
+            size,
+            _int64,
+            usm_type,
+            kernel_params=kernel_params,
+            param_types=param_types,
+            double_precision=ft == numpy.float64,
+        )
+
     def random(self, size=None, dtype=None, out=None, *, usm_type="device"):
         """
         Return random floats in the half-open interval [0.0, 1.0).
@@ -925,6 +1079,30 @@ def chisquare(self, df, size=None, dtype=None, *, usm_type="device"):
             (df, "df", _CONS_POSITIVE),
         )
 
+    def noncentral_chisquare(
+        self, df, nonc, size=None, dtype=None, *, usm_type="device"
+    ):
+        """
+        Draw samples from a noncentral chi-square distribution.
+
+        For full documentation refer to
+        :obj:`numpy.random.Generator.noncentral_chisquare`.
+
+        Limitations
+        -----------
+        `dtype` and `usm_type` are not in NumPy's signature.
+
+        """
+        return self._cont(
+            ri._noncentral_chisquare,
+            "noncentral_chisquare",
+            size,
+            dtype,
+            usm_type,
+            (df, "df", _CONS_POSITIVE),
+            (nonc, "nonc", _CONS_NON_NEGATIVE),
+        )
+
     def standard_t(self, df, size=None, dtype=None, *, usm_type="device"):
         """
         Draw samples from a standard Student's t distribution with `df`
@@ -968,6 +1146,31 @@ def f(self, dfnum, dfden, size=None, dtype=None, *, usm_type="device"):
             (dfden, "dfden", _CONS_POSITIVE),
         )
 
+    def noncentral_f(
+        self, dfnum, dfden, nonc, size=None, dtype=None, *, usm_type="device"
+    ):
+        """
+        Draw samples from the noncentral F distribution.
+
+        For full documentation refer to
+        :obj:`numpy.random.Generator.noncentral_f`.
+
+        Limitations
+        -----------
+        `dtype` and `usm_type` are not in NumPy's signature.
+
+        """
+        return self._cont(
+            ri._noncentral_f,
+            "noncentral_f",
+            size,
+            dtype,
+            usm_type,
+            (dfnum, "dfnum", _CONS_POSITIVE),
+            (dfden, "dfden", _CONS_POSITIVE),
+            (nonc, "nonc", _CONS_NON_NEGATIVE),
+        )
+
     def beta(self, a, b, size=None, dtype=None, *, usm_type="device"):
         """
         Draw samples from a beta distribution.
@@ -989,6 +1192,319 @@ def beta(self, a, b, size=None, dtype=None, *, usm_type="device"):
             (b, "b", _CONS_POSITIVE),
         )
 
+    def poisson(self, lam=1.0, size=None, *, usm_type="device"):
+        """
+        Draw samples from a Poisson distribution.
+
+        For full documentation refer to :obj:`numpy.random.Generator.poisson`.
+
+        Limitations
+        -----------
+        `usm_type` is not in NumPy's signature.
+
+        """
+
+        self._check_queues(lam)
+        _check_constraint(
+            lam,
+            "lam",
+            _CONS_POISSON,
+            _is_scalar_path(lam),
+        )
+        return self._discrete_sample(
+            ri._poisson, (lam,), size, usm_type, f32_bounds=(_F32_LT_2_63,)
+        )
+
+    def binomial(self, n, p, size=None, *, usm_type="device"):
+        """
+        Draw samples from a binomial distribution.
+
+        For full documentation refer to
+        :obj:`numpy.random.Generator.binomial`.
+
+        Limitations
+        -----------
+        `usm_type` is not in NumPy's signature.
+
+        """
+
+        self._check_queues(n, p)
+        scalar = _is_scalar_path(n, p)
+        # NumPy converts both parameters before checking either
+        n = _as_int64(n)
+        _check_constraint(p, "p", _CONS_BOUNDED_0_1, scalar)
+        _check_constraint(n, "n", _CONS_NON_NEGATIVE, scalar)
+        return self._discrete_sample(
+            ri._binomial, (n, p), size, usm_type, param_types=(_int64, None)
+        )
+
+    def negative_binomial(self, n, p, size=None, *, usm_type="device"):
+        """
+        Draw samples from a negative binomial distribution.
+
+        For full documentation refer to
+        :obj:`numpy.random.Generator.negative_binomial`.
+
+        Limitations
+        -----------
+        `usm_type` is not in NumPy's signature.
+
+        """
+
+        self._check_queues(n, p)
+        # NumPy takes scalar parameters as arrays here, so None gives NaN
+        _check_constraint(n, "n", _CONS_POSITIVE_NOT_NAN, False)
+        _check_constraint(p, "p", _CONS_BOUNDED_GT_0_1, False)
+        if not (
+            dpnp.is_supported_array_type(n) or dpnp.is_supported_array_type(p)
+        ):
+            dn = numpy.asarray(n, dtype=numpy.float64)
+            dp = numpy.asarray(p, dtype=numpy.float64)
+            with numpy.errstate(all="ignore"):
+                max_lam = (1 - dp) / dp * (dn + 10 * numpy.sqrt(dn))
+            if numpy.any(max_lam > _POISSON_LAM_MAX):
+                raise ValueError(
+                    "n too large or p too small, see "
+                    "Generator.negative_binomial Notes"
+                )
+        return self._discrete_sample(
+            ri._negative_binomial, (n, p), size, usm_type
+        )
+
+    def geometric(self, p, size=None, *, usm_type="device"):
+        """
+        Draw samples from the geometric distribution.
+
+        For full documentation refer to
+        :obj:`numpy.random.Generator.geometric`.
+
+        Limitations
+        -----------
+        `usm_type` is not in NumPy's signature.
+
+        """
+
+        self._check_queues(p)
+        _check_constraint(
+            p,
+            "p",
+            _CONS_BOUNDED_GT_0_1,
+            _is_scalar_path(p),
+        )
+        return self._discrete_sample(
+            ri._geometric, (p,), size, usm_type, f32_bounds=(_F32_GT_0,)
+        )
+
+    def logseries(self, p, size=None, *, usm_type="device"):
+        """
+        Draw samples from a logarithmic series distribution.
+
+        For full documentation refer to
+        :obj:`numpy.random.Generator.logseries`.
+
+        Limitations
+        -----------
+        `usm_type` is not in NumPy's signature.
+
+        """
+
+        self._check_queues(p)
+        _check_constraint(
+            p,
+            "p",
+            _CONS_BOUNDED_LT_0_1,
+            _is_scalar_path(p),
+        )
+        return self._discrete_sample(
+            ri._logseries, (p,), size, usm_type, f32_bounds=(_F32_LT_1,)
+        )
+
+    def zipf(self, a, size=None, *, usm_type="device"):
+        """
+        Draw samples from a Zipf distribution.
+
+        For full documentation refer to :obj:`numpy.random.Generator.zipf`.
+
+        Limitations
+        -----------
+        `usm_type` is not in NumPy's signature. The samples follow the Zipf
+        distribution truncated to the largest ``int64``, while NumPy's
+        rejection test rounds to accepting every sample above ``2**53``,
+        which makes its tail too heavy for `a` close to 1.
+
+        """
+
+        self._check_queues(a)
+        _check_constraint(a, "a", _CONS_GT_1, _is_scalar_path(a))
+        return self._discrete_sample(
+            ri._zipf, (a,), size, usm_type, f32_bounds=(_F32_GT_1,)
+        )
+
+    def hypergeometric(
+        self, ngood, nbad, nsample, size=None, *, usm_type="device"
+    ):
+        """
+        Draw samples from a hypergeometric distribution.
+
+        For full documentation refer to
+        :obj:`numpy.random.Generator.hypergeometric`.
+
+        Limitations
+        -----------
+        `usm_type` is not in NumPy's signature.
+
+        """
+
+        self._check_queues(ngood, nbad, nsample)
+        params = tuple(_as_int64(x) for x in (ngood, nbad, nsample))
+        names = ("ngood", "nbad", "nsample")
+        scalar = _is_scalar_path(*params)
+
+        for x in params[:2]:
+            xp, (xa,) = _check_arrays((x,))
+            if xp is not None and xp.any(xa >= _HYPERGEOM_MAX):
+                raise ValueError(
+                    "both ngood and nbad must be less than " f"{_HYPERGEOM_MAX}"
+                )
+        xp, (good, bad, sample) = _check_arrays(params)
+        if xp is not None and xp.any(good + bad < sample):
+            raise ValueError("ngood + nbad < nsample")
+        for x, name in zip(params, names):
+            _check_constraint(x, name, _CONS_NON_NEGATIVE, scalar)
+
+        return self._discrete_sample(
+            ri._hypergeometric,
+            params,
+            size,
+            usm_type,
+            param_types=(_int64,) * 3,
+        )
+
+    def integers(
+        self,
+        low,
+        high=None,
+        size=None,
+        dtype=numpy.int64,
+        endpoint=False,
+        *,
+        usm_type="device",
+    ):
+        """
+        Return random integers from `low` (inclusive) to `high` (exclusive),
+        or if ``endpoint=True``, `low` (inclusive) to `high` (inclusive).
+
+        For full documentation refer to
+        :obj:`numpy.random.Generator.integers`.
+
+        Limitations
+        -----------
+        `usm_type` is not in NumPy's signature. The values drawn differ from
+        NumPy's for the same range, as there is no buffering of 32-bit
+        draws for the narrow types.
+
+        """
+
+        return self._integers(
+            low, high, size, dtype, endpoint, usm_type, "integers"
+        )
+
+    def _integers(self, low, high, size, dtype, endpoint, usm_type, name):
+        """The implementation of :meth:`integers`, for the method `name`."""
+
+        self._check_queues(low, high)
+        validate_usm_type(usm_type, allow_none=False)
+
+        if high is None:
+            low, high = 0, low
+
+        dtype = numpy.dtype(dtype)
+        if dtype not in _integer_types:
+            raise TypeError(f"Unsupported dtype {dtype!r} for {name}")
+
+        # The bounds are passed to the kernel as a closed interval
+        if dtype == numpy.bool_:
+            lo_lim, hi_lim = 0, 1
+        else:
+            info = numpy.iinfo(dtype)
+            lo_lim, hi_lim = int(info.min), int(info.max)
+        dt_name = "bool" if dtype == numpy.bool_ else str(dtype)
+
+        def _check_bounds(low, high):
+            # the checks of NumPy, for host bounds only
+            if low is not None and numpy.any(numpy.less(low, lo_lim)):
+                raise ValueError(f"low is out of bounds for {dt_name}")
+            if high is not None and numpy.any(
+                numpy.greater(high, hi_lim if endpoint else hi_lim + 1)
+            ):
+                raise ValueError(f"high is out of bounds for {dt_name}")
+            if low is None or high is None:
+                return
+            if numpy.any(
+                numpy.greater(low, high)
+                if endpoint
+                else numpy.greater_equal(low, high)
+            ):
+                if not numpy.any(low):
+                    # low is 0 when only one bound is given
+                    raise ValueError("high < 0" if endpoint else "high <= 0")
+                raise ValueError("low > high" if endpoint else "low >= high")
+
+        if _is_array(low) or _is_array(high):
+            # Widen to 64 bits with the sign of the result type, so that an
+            # exclusive high of max + 1 survives the subtraction
+            wide = (
+                numpy.uint64
+                if dtype.kind == "u" or dtype == numpy.bool_
+                else numpy.int64
+            )
+            queue = self.sycl_queue
+            params = (low, high)
+
+            def _on_host(x):
+                # bounds on the device are not checked
+                return None if dpnp.is_supported_array_type(x) else x
+
+            _check_bounds(
+                *(
+                    None if x is None else numpy.asarray(x)
+                    for x in map(_on_host, params)
+                )
+            )
+
+            def _to_bound(x, shift):
+                if dpnp.is_supported_array_type(x):
+                    x = dpnp.asarray(x, dtype=wide, sycl_queue=queue)
+                    return (x - shift if shift else x).astype(dtype)
+                if not _is_array(x):
+                    return int(x) - shift
+                # exact integers, truncated as in NumPy, since max + 1 fits
+                # the wide type only after the shift
+                x = numpy.frompyfunc(int, 1, 1)(numpy.asarray(x)) - shift
+                x = dpnp.asarray(x.astype(wide), sycl_queue=queue)
+                return x.astype(dtype)
+
+            low = _to_bound(low, 0)
+            high = _to_bound(high, 0 if endpoint else 1)
+            return self._sample(
+                ri._integers,
+                params,
+                size,
+                dtype,
+                usm_type,
+                kernel_params=(low, high),
+            )
+
+        low = int(low)
+        high = int(high)
+        _check_bounds(low, high)
+        if not endpoint:
+            high -= 1
+
+        return self._sample(
+            ri._integers, (), size, dtype, usm_type, kernel_params=(low, high)
+        )
+
 
 def default_rng(seed=None, *, device=None, sycl_queue=None):
     """
diff --git a/dpnp/tests/AGENTS.md b/dpnp/tests/AGENTS.md
index 704d8373cd1..81fa736ac91 100644
--- a/dpnp/tests/AGENTS.md
+++ b/dpnp/tests/AGENTS.md
@@ -20,7 +20,7 @@
 - Error cases run on both libraries: parametrize over `xp` in `[dpnp, numpy]` so both must raise. Gate NumPy-version-dependent behaviour with `numpy_version()` (or `testing.with_requires` in CuPy tests).
 - Compute-follows-data error tests expect `(ValueError, ExecutionPlacementError)`: dpnp raises either, depending on the helper.
 - A new function is added to the matching parametrized lists in `test_sycl_queue.py` and `test_usm_type.py` (`test_1in_1out`, `test_2in_1out`, ...); elementwise ones also to `test_strides.py`.
-- `dpnp.random.Generator` values differ from NumPy's; test the samplers statistically, as the `test_moments` tests do. Restore any NumPy global random state a test changes, even when NumPy raises.
+- `dpnp.random.Generator` values differ from NumPy's; test the samplers statistically, as the `test_moments` tests do. Cover the fp32 path of the discrete samplers with the `discrete_precision` fixture of `test_random_generator.py`, and restore any NumPy global random state a test changes, even when NumPy raises.
 - A test that may exceed the device's maximum allocation size skips on `USMAllocationError` (from `dpnp.exceptions`), not `MemoryError`, which it does not derive from.
 - Markers are strict and only `slow` and `multi_gpu` are registered. Tests that need huge allocations are `@pytest.mark.slow`; they never run in CI by design.
 - Skip lists are for device- or backend-specific failures only; use the full parametrized node ID relative to `dpnp/` (`tests/test_x.py::TestY::test_z[param]`; `conftest.py` matches it exactly, so a `dpnp/tests/...` entry matches nothing), and update or drop the entries of a test you rename or remove.
diff --git a/dpnp/tests/test_random_generator.py b/dpnp/tests/test_random_generator.py
index fe68a855ec5..37c96704101 100644
--- a/dpnp/tests/test_random_generator.py
+++ b/dpnp/tests/test_random_generator.py
@@ -6,6 +6,7 @@
 from numpy.testing import (
     assert_allclose,
     assert_array_equal,
+    assert_equal,
     assert_raises,
     assert_raises_regex,
 )
@@ -13,9 +14,20 @@
 import dpnp
 from dpnp.random import Generator, Philox4x3210, default_rng
 
-from .helper import get_float_dtypes
+from .helper import get_float_dtypes, has_support_aspect64
 
 float_dtypes = get_float_dtypes(no_float16=True)
+integer_dtypes = [
+    numpy.bool_,
+    numpy.int8,
+    numpy.uint8,
+    numpy.int16,
+    numpy.uint16,
+    numpy.int32,
+    numpy.uint32,
+    numpy.int64,
+    numpy.uint64,
+]
 
 # the number of output elements each work-item of a kernel produces, with
 # Philox4x3210
@@ -203,8 +215,10 @@ class TestFloatMethods:
         ("logistic", (1.0, 2.0)),
         ("triangular", (0.0, 1.0, 3.0)),
         ("vonmises", (0.5, 2.0)),
+        ("noncentral_chisquare", (3.0, 2.0)),
         ("standard_t", (5.0,)),
         ("f", (5.0, 7.0)),
+        ("noncentral_f", (5.0, 7.0, 2.0)),
     ]
 
     @pytest.mark.parametrize("dtype", float_dtypes)
@@ -295,8 +309,17 @@ def test_random_range(self, dtype):
             ("triangular", (0.0, 0.0, 1.0), 1 / 3, 1 / 18),
             # uniform on the circle, whatever mu is, as in NumPy
             ("vonmises", (0.5, 0.0), 0.0, math.pi**2 / 3),
+            ("noncentral_chisquare", (3.0, 2.0), 5.0, 14.0),
+            ("noncentral_chisquare", (0.5, 1.0), 1.5, 5.0),
+            ("noncentral_chisquare", (3.0, 0.0), 3.0, 6.0),
             ("standard_t", (10.0,), 0.0, 1.25),
             ("f", (10.0, 30.0), 30 / 28, 2 * 900 * 38 / (10 * 28**2 * 26)),
+            (
+                "noncentral_f",
+                (10.0, 30.0, 2.0),
+                30 * 12 / (10 * 28),
+                2 * 9 * (144 + 14 * 28) / (28**2 * 26),
+            ),
         ],
     )
     def test_moments(self, name, params, mean, var, dtype):
@@ -344,6 +367,8 @@ def test_vonmises(self, kappa, dtype):
             ("wald", (1e-10, 1e10), 0, numpy.inf),
             ("triangular", (0.0, 1.0, 1.0), 0, 1),
             ("vonmises", (-3.0, 1e-3), -numpy.pi, numpy.pi),
+            ("noncentral_chisquare", (1e-3, 1e-3), 0, numpy.inf),
+            ("noncentral_f", (1e-3, 1e-3, 1e3), 0, numpy.inf),
         ],
     )
     def test_extreme_range(self, name, params, lo, hi, dtype):
@@ -364,8 +389,10 @@ def test_extreme_range(self, name, params, lo, hi, dtype):
             ("logistic", (0.0, 1.0)),
             ("vonmises", (0.0, 1.0)),
             ("chisquare", (1.0,)),
+            ("noncentral_chisquare", (1.0, 1.0)),
             ("standard_t", (1.0,)),
             ("f", (1.0, 1.0)),
+            ("noncentral_f", (1.0, 1.0, 1.0)),
             ("gamma", (1.0, 1.0)),
             ("beta", (1.0, 1.0)),
         ],
@@ -545,6 +572,425 @@ def test_host_params_terminate(self):
         assert (dpnp.abs(res) <= numpy.pi).all()
 
 
+class TestIntegers:
+    @pytest.mark.parametrize("dtype", integer_dtypes)
+    def test_full_range(self, dtype):
+        if dtype == numpy.bool_:
+            lo, hi = 0, 1
+        else:
+            info = numpy.iinfo(dtype)
+            lo, hi = int(info.min), int(info.max)
+
+        rng = default_rng(6)
+        res = rng.integers(lo, hi, size=10**5, dtype=dtype, endpoint=True)
+        assert res.dtype == dtype
+        res = dpnp.asnumpy(res)
+        assert res.min() >= lo and res.max() <= hi
+
+        if numpy.dtype(dtype).itemsize == 1:
+            # every value of a small type is drawn
+            assert len(numpy.unique(res)) == hi - lo + 1
+
+        # an exclusive high may be one past the maximum
+        res = rng.integers(lo, hi + 1, size=10, dtype=dtype)
+        assert res.dtype == dtype
+
+    @pytest.mark.parametrize("dtype", integer_dtypes[1:])
+    def test_bounds(self, dtype):
+        rng = default_rng(7)
+        res = dpnp.asnumpy(rng.integers(3, 10, size=10**5, dtype=dtype))
+        assert_equal(numpy.unique(res), numpy.arange(3, 10, dtype=dtype))
+
+        res = dpnp.asnumpy(
+            rng.integers(3, 10, size=10**5, dtype=dtype, endpoint=True)
+        )
+        assert_equal(numpy.unique(res), numpy.arange(3, 11, dtype=dtype))
+
+        res = dpnp.asnumpy(rng.integers(5, size=1000, dtype=dtype))
+        assert_equal(numpy.unique(res), numpy.arange(5, dtype=dtype))
+
+        # a single-value range
+        res = dpnp.asnumpy(rng.integers(4, 5, size=10, dtype=dtype))
+        assert (res == 4).all()
+
+    def test_uniformity(self):
+        # a range that is not a power of two is not biased
+        n, k = 10**6, 3
+        res = dpnp.asnumpy(default_rng(8).integers(0, k, size=n))
+        counts = numpy.bincount(res, minlength=k)
+        assert (abs(counts - n / k) < 5 * numpy.sqrt(n / k)).all()
+
+    def test_large_range(self):
+        rng = default_rng(9)
+        hi = (1 << 64) - 12345
+        res = dpnp.asnumpy(rng.integers(0, hi, size=10**4, dtype=numpy.uint64))
+        assert (res < hi).all()
+        # about half of the values are at least 2**63
+        assert 0.4 < (res >= (1 << 63)).mean() < 0.6
+
+    def test_default_dtype(self):
+        res = default_rng(0).integers(10, size=3)
+        assert res.dtype == dpnp.int64
+
+    def test_array_bounds(self):
+        rng = default_rng(10)
+        high = dpnp.asarray([1, 10, 100])
+        res = rng.integers(0, high, size=(1000, 3))
+        assert res.shape == (1000, 3)
+        res = dpnp.asnumpy(res)
+        assert (res < dpnp.asnumpy(high)).all()
+        assert res[:, 0].max() == 0
+        assert res[:, 2].max() > 50
+
+        # an exclusive high of max + 1 does not overflow
+        res = rng.integers(numpy.array([0, 250]), 256, dtype=numpy.uint8)
+        assert res.dtype == numpy.uint8
+        assert dpnp.asnumpy(res)[1] >= 250
+
+        res = rng.integers(-1, [3, 2**63], size=(1000, 2))
+        res = dpnp.asnumpy(res)
+        assert res[:, 0].min() == -1 and res[:, 0].max() == 2
+        assert (res[:, 1] >= -1).all() and (res[:, 1] > 2**62).any()
+
+    def test_float_array_bounds(self):
+        # truncated to integers, as in NumPy
+        res = default_rng(0).integers([0.5, -0.5], 2.5, size=(1000, 2))
+        assert_equal(numpy.unique(dpnp.asnumpy(res)), [0, 1])
+
+    @pytest.mark.parametrize(
+        "args, kwargs, msg",
+        [
+            ((5, 5), {}, "low >= high"),
+            ((5, 4), {"endpoint": True}, "low > high"),
+            ((-1, 5), {"dtype": numpy.uint8}, "low is out of bounds for uint8"),
+            (
+                (0, 257),
+                {"dtype": numpy.uint8},
+                "high is out of bounds for uint8",
+            ),
+            ((0, 256), {"dtype": numpy.uint8, "endpoint": True}, "high is out"),
+            ((0, 3), {"dtype": numpy.bool_}, "high is out of bounds for bool"),
+        ],
+    )
+    def test_errors(self, args, kwargs, msg):
+        rng = default_rng(0)
+        assert_raises_regex(ValueError, msg, rng.integers, *args, **kwargs)
+
+    def test_bad_dtype(self):
+        rng = default_rng(0)
+        assert_raises(TypeError, rng.integers, 0, 5, dtype=numpy.float32)
+
+
+def _zeta(a, n=10**5):
+    """The Riemann zeta function, by a partial sum and Euler-Maclaurin."""
+    k = numpy.arange(1, n, dtype=numpy.float64)
+    return (k**-a).sum() + n ** (1 - a) / (a - 1) + 0.5 * n**-a
+
+
+def _lchoose(n, k):
+    return math.lgamma(n + 1) - math.lgamma(k + 1) - math.lgamma(n - k + 1)
+
+
+def _xlog1py(x, y):
+    """x * log1p(y), taking 0 * log(0) to be 0."""
+    if x == 0:
+        return 0.0
+    return -math.inf if y == -1 else x * math.log1p(y)
+
+
+# log-pmfs of the discrete distributions, for the support values k
+_log_pmfs = {
+    "poisson": lambda k, lam: k * math.log(lam) - lam - math.lgamma(k + 1),
+    "binomial": lambda k, n, p: (
+        _lchoose(n, k) + k * math.log(p) + _xlog1py(n - k, -p)
+        if k <= n
+        else -math.inf
+    ),
+    "negative_binomial": lambda k, n, p: (
+        math.lgamma(k + n)
+        - math.lgamma(k + 1)
+        - math.lgamma(n)
+        + n * math.log(p)
+        + k * math.log1p(-p)
+    ),
+    "geometric": lambda k, p: (
+        _xlog1py(k - 1, -p) + math.log(p) if k >= 1 else -math.inf
+    ),
+    "logseries": lambda k, p: (
+        k * math.log(p) - math.log(k) - math.log(-math.log1p(-p))
+        if k >= 1
+        else -math.inf
+    ),
+    "zipf": lambda k, a: (
+        -a * math.log(k) - math.log(_zeta(a)) if k >= 1 else -math.inf
+    ),
+    "hypergeometric": lambda k, g, b, s: (
+        _lchoose(g, k) + _lchoose(b, s - k) - _lchoose(g + b, s)
+        if max(0, s - b) <= k <= min(g, s)
+        else -math.inf
+    ),
+}
+
+
+@pytest.fixture(params=["double", "single"])
+def discrete_precision(request, monkeypatch):
+    """Run a test with both floating-point types of discrete sampling."""
+
+    if request.param == "double" and not has_support_aspect64():
+        pytest.skip("no fp64 support on the device")
+    if request.param == "single":
+        monkeypatch.setattr(
+            Generator,
+            "_discrete_float_type",
+            lambda self: numpy.dtype(numpy.float32),
+        )
+    return request.param
+
+
+class TestDiscreteMethods:
+    cases = [
+        ("poisson", (3.0,)),
+        ("binomial", (10, 0.3)),
+        ("negative_binomial", (5.0, 0.3)),
+        ("geometric", (0.5,)),
+        ("logseries", (0.6,)),
+        ("zipf", (2.5,)),
+        ("hypergeometric", (5, 7, 3)),
+    ]
+
+    @pytest.mark.parametrize("name, params", cases)
+    def test_shape_dtype(self, name, params):
+        rng = default_rng(0)
+        meth = getattr(rng, name)
+
+        res = meth(*params, size=(3, 4))
+        assert isinstance(res, dpnp.ndarray)
+        assert res.shape == (3, 4)
+        assert res.dtype == dpnp.int64
+        assert meth(*params).shape == ()
+        assert meth(*params, size=0).shape == (0,)
+
+        res = meth(*params, size=2, usm_type="host")
+        assert res.usm_type == "host"
+
+    @pytest.mark.parametrize(
+        "name, params",
+        [
+            ("poisson", (0.5,)),
+            ("poisson", (5.0,)),
+            ("poisson", (10.0,)),
+            ("poisson", (15.5,)),
+            ("poisson", (100.0,)),
+            ("binomial", (10, 0.3)),
+            ("binomial", (20, 0.9)),
+            ("binomial", (100, 0.4)),
+            ("binomial", (1000, 0.5)),
+            ("binomial", (5, 1.0)),
+            ("negative_binomial", (5.0, 0.3)),
+            ("negative_binomial", (0.5, 0.8)),
+            ("negative_binomial", (100.0, 0.5)),
+            ("geometric", (0.5,)),
+            ("geometric", (0.9,)),
+            ("geometric", (0.1,)),
+            ("geometric", (1.0,)),
+            ("logseries", (0.3,)),
+            ("logseries", (0.9,)),
+            ("zipf", (1.5,)),
+            ("zipf", (2.5,)),
+            ("zipf", (6.0,)),
+            ("hypergeometric", (5, 7, 3)),
+            ("hypergeometric", (7, 5, 10)),
+            ("hypergeometric", (100, 200, 50)),
+            ("hypergeometric", (1000, 5, 100)),
+            ("hypergeometric", (300, 300, 590)),
+        ],
+    )
+    def test_pmf(self, name, params, discrete_precision):
+        n = 10**6
+        res = dpnp.asnumpy(getattr(default_rng(3), name)(*params, size=n))
+        counts = numpy.bincount(res[res < 10**4], minlength=1)
+
+        for k in range(len(counts) + 1):
+            prob = math.exp(_log_pmfs[name](k, *params))
+            freq = counts[k] / n if k < len(counts) else 0.0
+            # five standard errors of the frequency
+            assert (
+                abs(freq - prob) <= 5 * math.sqrt(prob * (1 - prob) / n) + 5e-6
+            )
+
+    @pytest.mark.parametrize(
+        "name, params, mean, var",
+        [
+            ("poisson", (1e7,), 1e7, 1e7),
+            ("poisson", (1e15,), 1e15, 1e15),
+            ("binomial", (10**9, 0.3), 3e8, 2.1e8),
+            ("binomial", (2**40, 0.5), 2.0**39, 2.0**38),
+            ("binomial", (10**7, 1e-6), 10.0, 10.0),
+            ("negative_binomial", (1e6, 0.5), 1e6, 2e6),
+            ("geometric", (1e-6,), 1e6, (1 - 1e-6) / 1e-12),
+            ("logseries", (0.999,), None, None),
+            (
+                "hypergeometric",
+                (10**8, 3 * 10**8, 10**6),
+                2.5e5,
+                10**6 * 0.25 * 0.75 * (4e8 - 1e6) / (4e8 - 1),
+            ),
+            (
+                "hypergeometric",
+                (5 * 10**8, 4 * 10**8, 3 * 10**8),
+                3e8 * 5 / 9,
+                3e8 * (5 / 9) * (4 / 9) * 6e8 / (9e8 - 1),
+            ),
+        ],
+    )
+    def test_large_params(self, name, params, mean, var, discrete_precision):
+        n = 10**6
+        res = dpnp.asnumpy(getattr(default_rng(4), name)(*params, size=n))
+        assert (res >= 0).all()
+        # the draws are not limited to the integers of the compute type,
+        # but for a deviation from the mean drawn by a float32 sampler
+        # above about 2**22
+        large = res[res >= 2**25]
+        if large.size > 1000 and (
+            discrete_precision == "double" or var < 2**44
+        ):
+            odd = (large % 2).mean()
+            assert abs(odd - 0.5) <= 5 * math.sqrt(0.25 / large.size)
+        if name == "logseries":
+            p = params[0]
+            mean = -p / ((1 - p) * math.log1p(-p))
+            var = -p * (p + math.log1p(-p)) / ((1 - p) * math.log1p(-p)) ** 2
+            # a heavy tail
+            var_tol = 0.2
+        else:
+            var_tol = 5 * math.sqrt(10 / n)
+
+        res = res.astype(numpy.float64)
+        # the rounding of the parameter to the compute type
+        eps = 2**-23 if discrete_precision == "single" else 2**-52
+        assert abs(res.mean() - mean) <= 5 * math.sqrt(var / n) + eps * mean
+        assert abs(res.var() - var) <= var_tol * var
+
+    @pytest.mark.parametrize("a", [1.1, 1.5])
+    def test_zipf_tail(self, a, discrete_precision):
+        # NumPy's (1 + 1/X)**(a - 1) - 1 rounds to 0 above 2**53, which makes
+        # its tail too heavy for a close to 1; this is the Zipf distribution
+        # truncated to INT64_MAX
+        n, x = 2**20, 2**25
+        # the sum of k**-a over k >= x, by Euler-Maclaurin
+        tail = lambda x: x ** (1 - a) / (a - 1) + 0.5 * x**-a
+        prob = (tail(x) - tail(2.0**63)) / (_zeta(a) - tail(2.0**63))
+        res = dpnp.asnumpy(default_rng(7).zipf(a, size=n))
+        count = (res >= x).sum()
+        assert abs(count - n * prob) <= 5 * math.sqrt(n * prob * (1 - prob))
+
+    def test_edge_values(self, discrete_precision):
+        rng = default_rng(5)
+        assert (dpnp.asnumpy(rng.poisson(0.0, size=100)) == 0).all()
+        assert (dpnp.asnumpy(rng.binomial(0, 0.5, size=100)) == 0).all()
+        assert (dpnp.asnumpy(rng.binomial(7, 0.0, size=100)) == 0).all()
+        assert (dpnp.asnumpy(rng.binomial(7, 1.0, size=100)) == 7).all()
+        assert (dpnp.asnumpy(rng.geometric(1.0, size=100)) == 1).all()
+        assert (dpnp.asnumpy(rng.logseries(0.0, size=100)) == 1).all()
+        assert (dpnp.asnumpy(rng.zipf(1025.0, size=100)) == 1).all()
+        assert (dpnp.asnumpy(rng.hypergeometric(5, 0, 3, size=100)) == 3).all()
+        assert (dpnp.asnumpy(rng.hypergeometric(0, 5, 3, size=100)) == 0).all()
+        assert (dpnp.asnumpy(rng.hypergeometric(5, 5, 10, size=100)) == 5).all()
+
+        # values that do not fit an int64 are clamped or rejected
+        res = dpnp.asnumpy(rng.geometric(1e-300, size=1000))
+        assert (res > 0).all()
+        res = dpnp.asnumpy(rng.zipf(1.0001, size=1000))
+        assert (res > 0).all()
+
+    @pytest.mark.parametrize(
+        "name, param",
+        [
+            ("geometric", [1e-300, 0.5]),
+            ("logseries", 1 - 1e-9),
+            ("zipf", 1 + 1e-9),
+            # NumPy's largest valid lam
+            ("poisson", 2**63 - 1 - 10 * math.sqrt(2**63 - 1)),
+        ],
+    )
+    def test_param_near_bound(self, name, param, discrete_precision):
+        # rounding to float32 must not take a valid parameter out of bounds
+        res = dpnp.asnumpy(getattr(default_rng(5), name)(param, size=(100, 2)))
+        assert (res > 0).all()
+
+    def test_array_params(self, discrete_precision):
+        rng = default_rng(6)
+        lam = dpnp.asarray([1.0, 100.0, 1e6])
+        res = rng.poisson(lam, size=(10**5, 3))
+        assert_allclose(
+            dpnp.asnumpy(res.mean(axis=0)), [1.0, 100.0, 1e6], rtol=0.02
+        )
+
+        res = rng.binomial([10, 1000], dpnp.asarray([[0.5], [0.1]]))
+        assert res.shape == (2, 2)
+
+        res = rng.hypergeometric(numpy.array([5, 500]), 5, [3, 7], size=(4, 2))
+        assert res.shape == (4, 2)
+        res = dpnp.asnumpy(res)
+        assert (res[:, 0] <= 3).all() and (res[:, 1] <= 7).all()
+
+        # integer parameters are truncated, as in NumPy
+        res = dpnp.asnumpy(rng.binomial([2.9], 1.0, size=3))
+        assert (res == 2).all()
+        res = dpnp.asnumpy(rng.binomial(dpnp.asarray([2.9]), 1.0, size=3))
+        assert (res == 2).all()
+
+    def test_device_params_not_checked(self):
+        # invalid parameters on the device give undefined values
+        rng = default_rng(0)
+        rng.poisson(dpnp.asarray([-1.0, numpy.nan]))
+        rng.binomial(dpnp.asarray([-1]), 0.5)
+        rng.zipf(dpnp.asarray([0.5, numpy.nan]))
+        rng.hypergeometric(dpnp.asarray([10**10]), 1, 5)
+        rng.normal(0, dpnp.asarray([-1.0]))
+
+    def test_device_params_terminate(self):
+        # the parameters on which the rejection samplers could loop forever
+        # give INT64_MIN
+        rng = default_rng(0)
+        nan, inf = numpy.nan, numpy.inf
+        int64_min = numpy.iinfo(numpy.int64).min
+        finfo = numpy.finfo(dpnp.default_float_type(sycl_queue=rng.sycl_queue))
+        fmax, tiny = float(finfo.max), float(finfo.tiny)
+
+        res = rng.poisson(dpnp.asarray([inf, fmax, 2.0**63]), size=(8, 3))
+        assert_array_equal(res, int64_min, strict=False)
+        res = rng.binomial(5, dpnp.asarray([nan, -1.0, 2.0]), size=(8, 3))
+        assert_array_equal(res, int64_min, strict=False)
+        n = dpnp.asarray([2**40, -1], dtype=dpnp.int64)
+        res = rng.hypergeometric(n, n.max(), 10**9, size=(8, 2))
+        assert_array_equal(res, int64_min, strict=False)
+        res = rng.logseries(dpnp.asarray([1.0, 2.0, inf, nan]), size=(8, 4))
+        assert_array_equal(res, int64_min, strict=False)
+        rng.negative_binomial(dpnp.asarray([inf, nan]), 0.5)
+        rng.noncentral_chisquare(1.0, dpnp.asarray([inf, nan]))
+        rng.noncentral_f(1.0, 1.0, dpnp.asarray([inf, nan]))
+
+        # a shape just below zero would reject nearly every draw
+        shape = dpnp.asarray([-1e-8, -1e-12, -1e-30, -1.0])
+        rng.negative_binomial(shape, 0.5)
+        rng.noncentral_chisquare(shape, dpnp.asarray([1e-30] * 4))
+
+        p = dpnp.asarray([nan, -0.0, -tiny, 0.0, -inf])
+        res = rng.geometric(p, size=(8, 5))
+        assert_array_equal(res, int64_min, strict=False)
+
+    def test_host_params_terminate(self, discrete_precision):
+        # valid parameters on which rounding in single precision made every
+        # draw of a rejection sampler fail
+        rng = default_rng(0)
+        res = rng.zipf([129.0, 200.0, 1024.0], size=(8, 3))
+        assert_array_equal(res, 1, strict=False)
+        res = rng.binomial([25 * 10**14, 10**17], [0.5, 0.05], size=(8, 2))
+        mean = numpy.array([12.5e14, 5e15])
+        assert_allclose(dpnp.asnumpy(res).mean(axis=0), mean, rtol=1e-3)
+
+
 def _assert_same_error(np_call, dp_call):
     """Assert `dp_call` raises the same exception as `np_call`."""
 
@@ -578,9 +1024,14 @@ class TestChecks:
             ("gamma", ([-1.0], 1.0)),
             ("chisquare", (0.0,)),
             ("chisquare", ([1.0, 0.0],)),
+            ("noncentral_chisquare", (0.0, 1.0)),
+            ("noncentral_chisquare", (1.0, -1.0)),
             ("standard_t", (0.0,)),
             ("f", (0.0, 1.0)),
             ("f", (1.0, [0.0])),
+            ("noncentral_f", (0.0, 1.0, 1.0)),
+            ("noncentral_f", (1.0, 0.0, 1.0)),
+            ("noncentral_f", (1.0, 1.0, -1.0)),
             ("beta", (0.0, 1.0)),
             ("beta", (1.0, [0.0])),
             ("pareto", (0.0,)),
@@ -599,11 +1050,61 @@ class TestChecks:
             ("triangular", (0.0, 2.0, 1.0)),
             ("triangular", (1.0, 1.0, 1.0)),
             ("triangular", ([0.0, 1.0], 1.0, [2.0, 1.0])),
+            ("poisson", (-1.0,)),
+            ("poisson", (nan,)),
+            ("poisson", ([nan],)),
+            ("poisson", ([-1.0],)),
+            ("poisson", (1e19,)),
+            ("poisson", ([1e19],)),
+            ("binomial", (-1, 0.5)),
+            ("binomial", ([-1], 0.5)),
+            ("binomial", (5, 1.5)),
+            ("binomial", (5, nan)),
+            ("binomial", ([5], [nan])),
+            ("binomial", (nan, 0.5)),
+            ("binomial", (-1, 2.0)),
+            ("negative_binomial", (0.0, 0.5)),
+            ("negative_binomial", (nan, 0.5)),
+            ("negative_binomial", ([nan], 0.5)),
+            ("negative_binomial", (1.0, 0.0)),
+            ("negative_binomial", (1.0, [1.5])),
+            ("negative_binomial", (-1.0, 2.0)),
+            ("negative_binomial", (1e20, 1e-10)),
+            ("negative_binomial", ([1e20], 1e-10)),
+            ("geometric", (0.0,)),
+            ("geometric", (1.5,)),
+            ("geometric", ([nan],)),
+            ("logseries", (1.0,)),
+            ("logseries", (-0.1,)),
+            ("logseries", (nan,)),
+            ("logseries", ([nan],)),
+            ("zipf", (1.0,)),
+            ("zipf", (nan,)),
+            ("zipf", ([nan],)),
+            ("hypergeometric", (10**9, 1, 1)),
+            ("hypergeometric", (1, [10**9], 1)),
+            ("hypergeometric", (1, 2, 4)),
+            ("hypergeometric", ([1], 2, 4)),
+            ("hypergeometric", (-1, 5, 3)),
+            ("hypergeometric", (5, -1, 3)),
+            ("hypergeometric", (5, 3, -1)),
+            ("hypergeometric", ([5], [3], [-1])),
+            ("hypergeometric", (nan, 3, 1)),
+            ("poisson", (None,)),
             ("normal", (0.0, None)),
             ("uniform", (None, 1.0)),
+            ("binomial", (5, None)),
             ("triangular", (0.0, None, 2.0)),
             ("exponential", ("x",)),
             ("exponential", (1j,)),
+            ("negative_binomial", (None, 0.5)),
+            ("binomial", (None, [0.5])),
+            ("integers", (1, [0, 1])),
+            ("integers", ([0, 1], 1)),
+            ("integers", (0, [0, 1])),
+            ("integers", (0, [-1, 1], None, numpy.int64, True)),
+            ("integers", ([-1], 3, None, numpy.uint8)),
+            ("integers", ([0], 256, None, numpy.uint8, True)),
         ],
     )
     @pytest.mark.filterwarnings("ignore:overflow encountered:RuntimeWarning")
@@ -636,9 +1137,14 @@ def _draw(self, device, bitgen):
         rng = Generator(bitgen(12, device=device))
         return (
             rng.random(1000, dtype=dpnp.float32),
+            rng.integers(0, 1000, size=1000),
             rng.bit_generator.random_raw(1000),
             rng.normal(size=1000, dtype=dpnp.float32),
             rng.gamma(0.7, size=1000, dtype=dpnp.float32),
+            rng.poisson([3.0, 50.0], size=(500, 2)),
+            rng.binomial(100, 0.3, size=1000),
+            rng.zipf(2.0, size=1000),
+            rng.hypergeometric(50, 70, 30, size=1000),
         )
 
     @pytest.mark.parametrize("bitgen", [Philox4x3210])
@@ -647,12 +1153,18 @@ def test_same_values(self, bitgen):
         ref = [dpnp.asnumpy(x) for x in self._draw(devices[0], bitgen)]
         for dev in devices[1:]:
             res = [dpnp.asnumpy(x) for x in self._draw(dev, bitgen)]
-            # uniform bits are exact
-            for a, b in zip(ref[:2], res[:2]):
+            # uniform bits and integers are exact
+            for a, b in zip(ref[:3], res[:3]):
                 assert_array_equal(a, b)
             # transcendental functions may differ in the last bits
-            for a, b in zip(ref[2:], res[2:]):
+            for a, b in zip(ref[3:5], res[3:5]):
                 assert_allclose(a, b, rtol=1e-4, atol=1e-5)
+            # and so may flip a rejection test, which shifts the rest of
+            # the work-item's stream; the discrete samplers compute in
+            # float32 without fp64, so their values differ from fp64 devices
+            if dev.has_aspect_fp64 == devices[0].has_aspect_fp64:
+                for a, b in zip(ref[5:], res[5:]):
+                    assert (a != b).mean() < 0.05
 
     def test_queue(self):
         for dev in _all_devices():
diff --git a/dpnp/tests/test_sycl_queue.py b/dpnp/tests/test_sycl_queue.py
index b7601cf06a5..45c423030ce 100644
--- a/dpnp/tests/test_sycl_queue.py
+++ b/dpnp/tests/test_sycl_queue.py
@@ -927,16 +927,25 @@ def test_random_state(func, args, kwargs, device, usm_type):
 
 _generator_methods = [
     pytest.param("beta", [0.5, 2.0], {"size": (3, 4)}),
+    pytest.param("binomial", [10, 0.3], {"size": (3, 4)}),
     pytest.param("chisquare", [3.0], {"size": (3, 4)}),
     pytest.param("exponential", [2.0], {"size": (3, 4)}),
     pytest.param("f", [2.0, 3.0], {"size": (3, 4)}),
     pytest.param("gamma", [2.0], {"size": (3, 4)}),
+    pytest.param("geometric", [0.3], {"size": (3, 4)}),
     pytest.param("gumbel", [], {"size": (3, 4)}),
+    pytest.param("hypergeometric", [5, 7, 4], {"size": (3, 4)}),
+    pytest.param("integers", [0, 10], {"size": (3, 4)}),
     pytest.param("laplace", [], {"size": (3, 4)}),
     pytest.param("logistic", [], {"size": (3, 4)}),
     pytest.param("lognormal", [], {"size": (3, 4)}),
+    pytest.param("logseries", [0.5], {"size": (3, 4)}),
+    pytest.param("negative_binomial", [3, 0.4], {"size": (3, 4)}),
+    pytest.param("noncentral_chisquare", [3.0, 1.0], {"size": (3, 4)}),
+    pytest.param("noncentral_f", [2.0, 3.0, 1.0], {"size": (3, 4)}),
     pytest.param("normal", [], {"size": (3, 4)}),
     pytest.param("pareto", [2.0], {"size": (3, 4)}),
+    pytest.param("poisson", [3.0], {"size": (3, 4)}),
     pytest.param("power", [2.0], {"size": (3, 4)}),
     pytest.param("random", [], {"size": (3, 4)}),
     pytest.param("rayleigh", [], {"size": (3, 4)}),
@@ -950,6 +959,7 @@ def test_random_state(func, args, kwargs, device, usm_type):
     pytest.param("vonmises", [0.0, 1.0], {"size": (3, 4)}),
     pytest.param("wald", [1.0, 2.0], {"size": (3, 4)}),
     pytest.param("weibull", [2.0], {"size": (3, 4)}),
+    pytest.param("zipf", [2.0], {"size": (3, 4)}),
 ]
 
 
@@ -983,10 +993,17 @@ def test_bit_generator(bit_generator, device):
 # "f", "i" and "p" stand for float, integer and probability arrays on another
 # queue than the generator's
 _diff_queue_methods = [
+    pytest.param("binomial", ["i", 0.5]),
+    pytest.param("geometric", ["p"]),
+    pytest.param("hypergeometric", ["i", "i", 1]),
+    pytest.param("logseries", ["p"]),
+    pytest.param("negative_binomial", ["f", 0.5]),
     pytest.param("normal", ["f"]),
+    pytest.param("poisson", ["f"]),
     pytest.param("standard_gamma", ["f"]),
     pytest.param("triangular", [0.0, "f", 5.0]),
     pytest.param("uniform", ["f", 5.0]),
+    pytest.param("zipf", ["f"]),
 ]
 
 
@@ -1004,6 +1021,7 @@ def _diff_queue_args(args, device):
     "func, args",
     _diff_queue_methods
     + [
+        pytest.param("integers", ["i"]),
         pytest.param("random", [4, None], id="random-out"),
     ],
 )
diff --git a/dpnp/tests/test_usm_type.py b/dpnp/tests/test_usm_type.py
index 8fda5a32b58..90eb7c938e1 100644
--- a/dpnp/tests/test_usm_type.py
+++ b/dpnp/tests/test_usm_type.py
@@ -1837,7 +1837,10 @@ def test_tensorsolve(self, usm_type, usm_type_other):
 @pytest.mark.parametrize(
     "func, args, kwargs",
     [
+        pytest.param("hypergeometric", [5, 7, 4], {"size": 3}),
+        pytest.param("integers", [[3, 5], [7, 9]], {}),
         pytest.param("normal", [[1.0, 2.0]], {}),
+        pytest.param("poisson", [3.0], {"size": 3}),
         pytest.param("random", [3], {}),
         pytest.param("standard_gamma", [2.0], {"size": 3}),
         pytest.param("triangular", [0.0, 1.0, 3.0], {"size": 3}),