diff --git a/.github/actions/setup-randblas-deps-windows/action.yml b/.github/actions/setup-randblas-deps-windows/action.yml index a21c2ff0..a0046c16 100644 --- a/.github/actions/setup-randblas-deps-windows/action.yml +++ b/.github/actions/setup-randblas-deps-windows/action.yml @@ -17,9 +17,6 @@ outputs: blaspp-dir: description: "Directory containing blasppConfig.cmake." value: ${{ steps.setup.outputs.blaspp-dir }} - random123-dir: - description: "Directory containing the installed Random123 headers." - value: ${{ steps.setup.outputs.random123-dir }} googletest-prefix: description: "GoogleTest installation prefix." value: ${{ steps.setup.outputs.googletest-prefix }} @@ -59,12 +56,6 @@ runs: path: ${{ github.workspace }}\..\windows-deps\blaspp-install key: windows-msvc-blaspp-windows-portability-ilp64-sequential-3-${{ hashFiles('.github/actions/setup-randblas-deps-windows/setup.ps1') }} - - name: cache Random123 installation - uses: actions/cache@v4 - with: - path: ${{ github.workspace }}\..\windows-deps\Random123-install - key: windows-random123-2-${{ hashFiles('.github/actions/setup-randblas-deps-windows/setup.ps1') }} - - name: cache LAPACK++ installation if: inputs.install-lapackpp == 'true' uses: actions/cache@v4 diff --git a/.github/actions/setup-randblas-deps-windows/setup.ps1 b/.github/actions/setup-randblas-deps-windows/setup.ps1 index 2a9d1033..9aae2490 100644 --- a/.github/actions/setup-randblas-deps-windows/setup.ps1 +++ b/.github/actions/setup-randblas-deps-windows/setup.ps1 @@ -172,17 +172,6 @@ if (-not (Test-Path -LiteralPath (Join-Path $gtestInstall "lib\cmake\GTest\GTest ) } -$random123Source = Join-Path $DependencyRoot "Random123" -$random123Install = Join-Path $DependencyRoot "Random123-install" -$random123Include = Join-Path $random123Install "include" -if (-not (Test-Path -LiteralPath (Join-Path $random123Include "Random123\philox.h"))) { - Clone-Head -Url "https://github.com/DEShawResearch/Random123.git" ` - -Destination $random123Source - New-Item -ItemType Directory -Force -Path $random123Include | Out-Null - Copy-Item -LiteralPath (Join-Path $random123Source "include\Random123") ` - -Destination $random123Include -Recurse -} - $blasppSource = Join-Path $DependencyRoot "blaspp" $blasppBuild = Join-Path $DependencyRoot "blaspp-build" $blasppInstall = Join-Path $DependencyRoot "blaspp-install" @@ -254,7 +243,6 @@ if ($InstallLapackpp) { $exports = [ordered]@{ "blaspp_DIR" = Convert-ToCMakePath $blasppDir - "Random123_DIR" = Convert-ToCMakePath $random123Include "googletest_PREFIX" = Convert-ToCMakePath $gtestInstall "MKLROOT" = Convert-ToCMakePath $mklRoot } @@ -275,8 +263,6 @@ if ($env:GITHUB_PATH) { if ($env:GITHUB_OUTPUT) { "blaspp-dir=$(Convert-ToCMakePath $blasppDir)" | Out-File -FilePath $env:GITHUB_OUTPUT -Append -Encoding utf8 - "random123-dir=$(Convert-ToCMakePath $random123Include)" | - Out-File -FilePath $env:GITHUB_OUTPUT -Append -Encoding utf8 "googletest-prefix=$(Convert-ToCMakePath $gtestInstall)" | Out-File -FilePath $env:GITHUB_OUTPUT -Append -Encoding utf8 "mkl-root=$(Convert-ToCMakePath $mklRoot)" | diff --git a/.github/actions/setup-randblas-deps/action.yml b/.github/actions/setup-randblas-deps/action.yml index 18182d7d..5ca4fa79 100644 --- a/.github/actions/setup-randblas-deps/action.yml +++ b/.github/actions/setup-randblas-deps/action.yml @@ -1,8 +1,8 @@ name: setup-randblas-deps description: > - Installs RandBLAS build dependencies (GoogleTest, BLAS++, Random123, and - optionally LAPACK++), with caching for the from-source dependencies. Exports - blaspp_DIR, Random123_DIR, and (optionally) lapackpp_DIR via GITHUB_ENV. + Installs RandBLAS build dependencies (GoogleTest, BLAS++, and optionally + LAPACK++), with caching for the from-source dependencies. Exports blaspp_DIR + and (optionally) lapackpp_DIR via GITHUB_ENV. inputs: cc: @@ -188,26 +188,12 @@ runs: fi make -j"${jobs}" install - - name: install Random123 headers - shell: bash - run: | - set -euxo pipefail - cd .. - rm -rf Random123 Random123-install - git clone --depth 1 https://github.com/DEShawResearch/Random123.git - # Random123's `make install-include` target uses GNU `cp -d`, which - # isn't supported by BSD cp on macOS. Copy the headers directly so - # the same step works on both platforms. - mkdir -p Random123-install/include - cp -R Random123/include/Random123 Random123-install/include/ - - name: export dependency paths shell: bash run: | set -euxo pipefail cd .. echo "blaspp_DIR=$(pwd)/blaspp-install/lib/cmake/blaspp" >> "$GITHUB_ENV" - echo "Random123_DIR=$(pwd)/Random123-install/include" >> "$GITHUB_ENV" if [[ "${{ inputs.install-lapackpp }}" == "true" ]]; then echo "lapackpp_DIR=$(pwd)/lapackpp-install/lib/cmake/lapackpp" >> "$GITHUB_ENV" fi diff --git a/.github/scripts/windows/run-ci.ps1 b/.github/scripts/windows/run-ci.ps1 index 9e03ac3d..f0a48f10 100644 --- a/.github/scripts/windows/run-ci.ps1 +++ b/.github/scripts/windows/run-ci.ps1 @@ -78,7 +78,6 @@ if ($SetupDependencies) { } $blasppDir = Require-EnvironmentVariable "blaspp_DIR" -$random123Dir = Require-EnvironmentVariable "Random123_DIR" $mklRoot = Require-EnvironmentVariable "MKLROOT" $mklBin = Join-Path $mklRoot "bin" if (-not (Test-Path -LiteralPath $mklBin)) { @@ -98,10 +97,8 @@ function Install-RandBLAS { $build = Join-Path $WorkRoot "$Name-build" $install = Join-Path $WorkRoot "$Name-install" $configuredBlasppDir = $blasppDir - $configuredRandom123Dir = $random123Dir if ($UseNativeDependencyPaths) { $configuredBlasppDir = $blasppDir.Replace("/", "\") - $configuredRandom123Dir = $random123Dir.Replace("/", "\") } $arguments = @( "-S", $SourceRoot, @@ -110,7 +107,6 @@ function Install-RandBLAS { "-DCMAKE_BUILD_TYPE=Release", "-DCMAKE_INSTALL_PREFIX=$(Convert-ToCMakePath $install)", "-Dblaspp_DIR=$configuredBlasppDir", - "-DRandom123_DIR=$configuredRandom123Dir", "-DBUILD_TESTS=$(if ($BuildTests) { 'ON' } else { 'OFF' })" ) if ($BuildTests) { @@ -165,9 +161,9 @@ switch ($Task) { -UseNativeDependencyPaths $true $build = Join-Path $WorkRoot "downstream-consumer-build" - # Deliberately omit blaspp_DIR and Random123_DIR. This makes the smoke - # test exercise the native-backslash dependency paths recorded by - # RandBLASConfig.cmake and guards their generated-path normalization. + # Deliberately omit blaspp_DIR. This makes the smoke test exercise the + # native-backslash dependency path recorded by RandBLASConfig.cmake + # and guards its generated-path normalization. Invoke-Checked -Program "cmake" -Arguments @( "-S", (Join-Path $SourceRoot "test\downstream"), "-B", $build, @@ -198,7 +194,6 @@ switch ($Task) { "-DCMAKE_BUILD_TYPE=Release", "-DRandBLAS_DIR=$(Convert-ToCMakePath $randblasDir)", "-Dblaspp_DIR=$blasppDir", - "-DRandom123_DIR=$random123Dir", "-Dlapackpp_DIR=$lapackppDir" ) Invoke-Checked -Program "cmake" -Arguments @("--build", $build) diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index 639f707c..b80a787b 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -206,7 +206,6 @@ jobs: -DCMAKE_EXE_LINKER_FLAGS="${extra_link_flags}" \ ${sanitize_address_arg} \ -Dblaspp_DIR="${blaspp_DIR}" \ - -DRandom123_DIR="${Random123_DIR}" \ -DCMAKE_INSTALL_PREFIX="$(pwd)/../RandBLAS-install" \ "$(pwd)/../RandBLAS" diff --git a/.github/workflows/downstream-consumer.yml b/.github/workflows/downstream-consumer.yml index ef1768e1..d7adb61c 100644 --- a/.github/workflows/downstream-consumer.yml +++ b/.github/workflows/downstream-consumer.yml @@ -33,7 +33,6 @@ jobs: cmake \ -DCMAKE_BUILD_TYPE=Release \ -Dblaspp_DIR="${blaspp_DIR}" \ - -DRandom123_DIR="${Random123_DIR}" \ -DCMAKE_INSTALL_PREFIX="$(pwd)/../RandBLAS-install" \ "$(pwd)/../RandBLAS" make -j"$(nproc)" install @@ -49,7 +48,6 @@ jobs: -DCMAKE_BUILD_TYPE=Release \ -DCMAKE_PREFIX_PATH="$(pwd)/../../../../RandBLAS-install" \ -Dblaspp_DIR="${blaspp_DIR}" \ - -DRandom123_DIR="${Random123_DIR}" \ .. make -j"$(nproc)" ./smoke diff --git a/.github/workflows/examples.yml b/.github/workflows/examples.yml index 25663134..e9a732a5 100644 --- a/.github/workflows/examples.yml +++ b/.github/workflows/examples.yml @@ -34,7 +34,6 @@ jobs: cmake \ -DCMAKE_BUILD_TYPE=Release \ -Dblaspp_DIR="${blaspp_DIR}" \ - -DRandom123_DIR="${Random123_DIR}" \ -DCMAKE_INSTALL_PREFIX="$(pwd)/../RandBLAS-install" \ "$(pwd)/../RandBLAS" make -j"$(nproc)" install @@ -51,7 +50,6 @@ jobs: -DCMAKE_PREFIX_PATH="$(pwd)/../../../RandBLAS-install" \ -Dblaspp_DIR="${blaspp_DIR}" \ -Dlapackpp_DIR="${lapackpp_DIR}" \ - -DRandom123_DIR="${Random123_DIR}" \ .. make -j"$(nproc)" diff --git a/.github/workflows/thread-sanitizer.yml b/.github/workflows/thread-sanitizer.yml index 1e0f42ba..23f7cabe 100644 --- a/.github/workflows/thread-sanitizer.yml +++ b/.github/workflows/thread-sanitizer.yml @@ -47,7 +47,6 @@ jobs: -DCMAKE_CXX_FLAGS="-fsanitize=thread -fno-omit-frame-pointer" \ -DCMAKE_EXE_LINKER_FLAGS="-fsanitize=thread" \ -Dblaspp_DIR="${blaspp_DIR}" \ - -DRandom123_DIR="${Random123_DIR}" \ "$(pwd)/../RandBLAS" make -j"$(nproc)" diff --git a/AGENTS.md b/AGENTS.md index af91f65e..e6e46b26 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -177,7 +177,6 @@ ctest -V # Verbose output Key CMake variables: - `blaspp_DIR`: Path to BLAS++ installation (containing `blasppConfig.cmake`) -- `Random123_DIR`: Path to Random123 headers - `CMAKE_BUILD_TYPE`: Release or Debug ### Installation @@ -186,7 +185,6 @@ Key CMake variables: mkdir RandBLAS-build && cd RandBLAS-build cmake -DCMAKE_BUILD_TYPE=Release \ -Dblaspp_DIR=/path/to/blaspp-install/lib/cmake/blaspp/ \ - -DRandom123_DIR=/path/to/random123-install/include/ \ ../RandBLAS/ make -j install ctest diff --git a/CMake/FindRandom123.cmake b/CMake/FindRandom123.cmake deleted file mode 100644 index f46ffe17..00000000 --- a/CMake/FindRandom123.cmake +++ /dev/null @@ -1,34 +0,0 @@ -if (NOT Random123_FOUND) - -# find the header -# first look where the user told us -if(Random123_DIR) - find_path(Random123_INCLUDE_DIR Random123/philox.h - PATHS "${Random123_DIR}" "${Random123_DIR}/include/" - NO_DEFAULT_PATH) -endif() - -# look in typical system locations -find_path(Random123_INCLUDE_DIR Random123/philox.h - PATHS "/usr/include/" "/usr/local/include/" - NO_DEFAULT_PATH) - -# finally let CMake look -find_path(Random123_INCLUDE_DIR Random123/philox.h) - -mark_as_advanced(Random123_INCLUDE_DIR) - -# handle the QUIETLY and REQUIRED arguments and set Random123_FOUND -include(FindPackageHandleStandardArgs) - -find_package_handle_standard_args(Random123 - "Failed to find Random123. Set -DRandom123_DIR=X with X pointing to the directory where the header files \"Random123/*.h\" are located." - Random123_INCLUDE_DIR) - -if (NOT TARGET Random123::Random123) - add_library(Random123::Random123 INTERFACE IMPORTED GLOBAL) - target_include_directories(Random123::Random123 - SYSTEM INTERFACE "${Random123_INCLUDE_DIR}") -endif() - -endif() diff --git a/CMake/RandBLASConfig.cmake.in b/CMake/RandBLASConfig.cmake.in index 3ac3f89a..06cd954e 100644 --- a/CMake/RandBLASConfig.cmake.in +++ b/CMake/RandBLASConfig.cmake.in @@ -14,12 +14,6 @@ if (NOT blaspp_DIR) endif () find_dependency(blaspp) -# Random123 -if (NOT Random123_DIR) - set(Random123_DIR "@RandBLAS_CONFIG_RANDOM123_DIR@") -endif () -find_dependency(Random123) - # OpenMP set(RandBLAS_HAS_OpenMP @RandBLAS_HAS_OpenMP@) if (RandBLAS_HAS_OpenMP) diff --git a/CMake/rb_config.cmake b/CMake/rb_config.cmake index 83e2c903..91ca107c 100644 --- a/CMake/rb_config.cmake +++ b/CMake/rb_config.cmake @@ -2,7 +2,6 @@ # syntax. In particular, native Windows backslashes would otherwise be parsed # as escape sequences when a downstream project loads RandBLASConfig.cmake. file(TO_CMAKE_PATH "${blaspp_DIR}" RandBLAS_CONFIG_BLASPP_DIR) -file(TO_CMAKE_PATH "${Random123_DIR}" RandBLAS_CONFIG_RANDOM123_DIR) configure_file(CMake/RandBLASConfig.cmake.in ${CMAKE_INSTALL_LIBDIR}/cmake/RandBLAS/RandBLASConfig.cmake @ONLY) @@ -10,11 +9,6 @@ configure_file(CMake/RandBLASConfig.cmake.in configure_file(CMake/RandBLASConfigVersion.cmake.in ${CMAKE_INSTALL_LIBDIR}/cmake/RandBLAS/RandBLASConfigVersion.cmake @ONLY) -if (PROJECT_NAME STREQUAL "RandBLAS") - install(FILES CMake/FindRandom123.cmake - DESTINATION "${CMAKE_INSTALL_LIBDIR}/cmake/RandBLAS") -endif() - install(FILES ${CMAKE_CURRENT_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR}/cmake/RandBLAS/RandBLASConfig.cmake ${CMAKE_CURRENT_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR}/cmake/RandBLAS/RandBLASConfigVersion.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index 1d9dfeb4..f0afa6f0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,7 +12,6 @@ include(RuntimeDLLs) # find dependencies find_package(blaspp REQUIRED) -find_package(Random123 REQUIRED) include(OpenMP) include(MKL_sparse) diff --git a/INSTALL.md b/INSTALL.md index 55cd69e5..25d5d438 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -14,7 +14,7 @@ If you want a TL;DR version of this guide, refer to one of the following. * The [examples folder](https://github.com/BallisticLA/RandBLAS/tree/main/examples). -## 1. Required dependencies: a C++20 compatible compiler, BLAS++, and Random123 +## 1. Required dependencies: a C++20 compatible compiler and BLAS++ RandBLAS uses C++20 [concepts](https://en.cppreference.com/w/cpp/language/constraints). Make sure your compiler supports these. We test gcc ≥13 on Linux, and both Apple Clang @@ -26,10 +26,11 @@ It can be installed with GNU make or CMake. If you want to use RandBLAS' CMake build system, then it will be necessary to have built and installed BLAS++ via CMake. -Random123 is a header-only library of counter-based random number generators. +RandBLAS includes its header-only counter-based random-number generators. There +is no separate random-number package to install or configure. -We give recipes for installing BLAS++ and Random123 below. -Later on, we'll assume these recipes were executed from a directory +We give a recipe for installing BLAS++ below. +Later on, we'll assume this recipe was executed from a directory that contains (or will contain) the ``RandBLAS`` project directory as a subdirectory. One can compile and install BLAS++ from @@ -47,14 +48,6 @@ cmake -DCMAKE_BUILD_TYPE=Release \ make -j install ``` -One can install Random123 from -[source](https://github.com/DEShawResearch/random123) by running -```shell -git clone git@github.com:DEShawResearch/random123.git -cd random123/ -make prefix=`pwd`/../random123-install install-include -``` - ## 2. Optional dependencies: GTest and OpenMP GoogleTest (aka *GTest*) is Google’s C++ testing and mocking framework. It is an optional @@ -93,7 +86,6 @@ The following CMake variables influence the RandBLAS build. |------------------|-------------------------------------------| | CMAKE_BUILD_TYPE | Release or Debug. The default is Release. | | blaspp_DIR | The path to your local BLAS++ install | -| Random123_DIR | The path to your local random123 install | Assuming you used the recipes from Section 1 to get RandBLAS' dependencies, you can download, build, and install RandBLAS as follows: @@ -104,7 +96,6 @@ mkdir RandBLAS-build cd RandBLAS-build cmake -DCMAKE_BUILD_TYPE=Release \ -Dblaspp_DIR=`pwd`/../blaspp-install/lib/cmake/blaspp/ \ - -DRandom123_DIR=`pwd`/../random123-install/include/ \ -DCMAKE_BINARY_DIR=`pwd` \ -DCMAKE_INSTALL_PREFIX=`pwd`/../RandBLAS-install \ ../RandBLAS/ @@ -116,9 +107,6 @@ Here are the conceptual meanings of the recipe's other build flags: * `-Dblaspp_DIR=X` means `X` is the directory containing the file `blasppConfig.cmake`. -* `-DRandom123_DIR=Y` means `Y` is the directory containing the Random123 - header files. - * `-DCMAKE_INSTALL_PREFIX=Z` means subdirectories within `Z` will contain the RandBLAS binaries, header files, and CMake configuration files needed for using RandBLAS in other projects. The CMake configuration files are @@ -153,7 +141,6 @@ find_package(lapackpp REQUIRED) set(myproject_cxx_source my_project.cc) add_executable(my_project ${myproject_cxx_source}) -target_include_directories(myproject PUBLIC ${Random123_DIR}) target_link_libraries(myproject PUBLIC RandBLAS blaspp lapackpp) ``` @@ -169,7 +156,7 @@ Run the commands from an **x64 Native Tools Command Prompt for Visual Studio**. The recipe uses the NMake generator, so `cmake --build` is serial and does not need `--parallel`. -### A.1. Required dependencies: MSVC, oneMKL, BLAS++, and Random123 +### A.1. Required dependencies: MSVC, oneMKL, and BLAS++ Install the following tools first: @@ -223,20 +210,6 @@ cmake --fresh ^ cmake --build C:/randblas-work/build/blaspp --target install ``` -Random123 is header-only. Copy its public headers into a stable installation -prefix so the installed RandBLAS package does not depend on retaining the -Random123 source checkout: - -```bat -git clone https://github.com/DEShawResearch/Random123.git ^ - C:\randblas-work\src\Random123 - -cmake -E make_directory C:/randblas-work/install/Random123/include -cmake -E copy_directory ^ - C:/randblas-work/src/Random123/include/Random123 ^ - C:/randblas-work/install/Random123/include/Random123 -``` - ### A.2. Optional dependencies: GoogleTest and OpenMP GoogleTest is needed only to build and run the RandBLAS test suite. A minimal @@ -269,8 +242,8 @@ command in the next section. ### A.3. Building, installing, and testing RandBLAS -Clone RandBLAS, configure it against the installed BLAS++ and the Random123 -headers, and enable the test suite: +Clone RandBLAS, configure it against the installed BLAS++, and enable the test +suite: ```bat git clone https://github.com/BallisticLA/RandBLAS.git ^ @@ -283,7 +256,6 @@ cmake --fresh ^ -DCMAKE_BUILD_TYPE=Release ^ -DCMAKE_INSTALL_PREFIX=C:/randblas-work/install/RandBLAS ^ -Dblaspp_DIR=C:/randblas-work/install/blaspp/blaspp ^ - -DRandom123_DIR=C:/randblas-work/install/Random123/include ^ -DCMAKE_PREFIX_PATH=C:/randblas-work/install/googletest ^ -DBUILD_TESTS=ON @@ -339,9 +311,8 @@ cmake --fresh ^ cmake --build C:/path/to/my_randblas_project-build ``` -The installed `RandBLASConfig.cmake` records the BLAS++ and Random123 locations -used to build RandBLAS, so ordinary consumers should not need to supply those -paths again. +The installed `RandBLASConfig.cmake` records the BLAS++ location used to build +RandBLAS, so ordinary consumers should not need to supply that path again. Repository-owned tests copy imported dependency DLLs beside their executables. An arbitrary downstream application is responsible for its own deployment. diff --git a/RandBLAS/CMakeLists.txt b/RandBLAS/CMakeLists.txt index b5fd2967..c586da38 100644 --- a/RandBLAS/CMakeLists.txt +++ b/RandBLAS/CMakeLists.txt @@ -1,5 +1,5 @@ -set(RandBLAS_libs blaspp Random123::Random123) +set(RandBLAS_libs blaspp) if (RandBLAS_HAS_OpenMP) list(APPEND RandBLAS_libs OpenMP::OpenMP_CXX) endif() @@ -19,19 +19,12 @@ target_compile_features(RandBLAS INTERFACE cxx_std_20) # RandBLAS is header-only, so these MSVC requirements must propagate to every # translation unit that includes its headers. /Zc:__cplusplus makes MSVC report -# the selected language standard correctly (which Random123 inspects), while -# /EHsc enables standard C++ exception-unwind semantics for code using the API. +# the selected language standard correctly, while /EHsc enables standard C++ +# exception-unwind semantics for code using the API. target_compile_options(RandBLAS INTERFACE $<$:/EHsc> $<$:/Zc:__cplusplus>) -# Random123's boxmuller.hpp otherwise calls the nonstandard sincos/sincosf -# functions, which are unavailable in the MSVC runtime. Select Random123's -# portable separate-sine-and-cosine fallback for native MSVC consumers. -target_compile_definitions(RandBLAS INTERFACE - $<$:R123_NO_SINCOS=1>) - - target_include_directories(RandBLAS INTERFACE $ $ diff --git a/RandBLAS/DevNotes.md b/RandBLAS/DevNotes.md index 30b02218..08e67a3d 100644 --- a/RandBLAS/DevNotes.md +++ b/RandBLAS/DevNotes.md @@ -4,8 +4,8 @@ This file reviews aspects of RandBLAS' implementation that aren't (currently) su for our user guide. - * Our basic random number generation is handled by [Random123](https://github.com/DEShawResearch/random123). - We have small wrappers around Random123 code in ``RandBLAS/base.hh`` and ``RandBLAS/random_gen.hh``. + * The native counter-based random-number subsystem is documented in + [``RandBLAS/rng/DevNotes.md``](rng/DevNotes.md). * ``RandBLAS/dense_skops.hh`` has code for representing and sampling dense sketching operators. The sampling code is complicated because it supports multi-threaded random (sub)matrix generation, and yet the generated (sub)matrices are the same no matter how many threads diff --git a/RandBLAS/base.hh b/RandBLAS/base.hh index f7f60c21..eaf606c4 100644 --- a/RandBLAS/base.hh +++ b/RandBLAS/base.hh @@ -36,7 +36,6 @@ #include #include -#include #include #include @@ -50,141 +49,32 @@ /// code common across the project namespace RandBLAS { -typedef r123::Philox4x32 DefaultRNG; using std::uint64_t; - -/// ------------------------------------------------------------------- -/// This is a stateful version of a -/// *counter-based random number generator* (CBRNG) from Random123. -/// It packages a CBRNG together with two arrays, called "counter" and "key," -/// which are interpreted as extended-width unsigned integers. -/// -/// RNGStates are used in every RandBLAS function that involves random sampling. -/// -template -struct RNGState { - - /// ------------------------------------------------------------------- - /// Type of the underlying Random123 CBRNG. Must be based on - /// Philox or Threefry. We've found that Philox works best for our - /// purposes, and we default to Philox4x32. - using generator = RNG; - - using ctr_type = typename RNG::ctr_type; - // ^ An array type defined in Random123. - using key_type = typename RNG::key_type; - // ^ An array type defined in Random123. - using ctr_uint = typename RNG::ctr_type::value_type; - // ^ The unsigned integer type used in this RNGState's counter array. - using key_uint = typename RNG::key_type::value_type; - // ^ The unsigned integer type used in this RNGState's key array. - - /// ------------------------------------------------------------------ - /// This is a Random123-defined statically-sized array of unsigned integers. - /// With RandBLAS' default, it contains four 32-bit unsigned ints - /// and is interpreted as one 128-bit unsigned int. - /// - /// This member specifies a "location" in the random stream - /// defined by RNGState::generator and RNGState::key. - /// Random sampling functions in RandBLAS effectively consume elements - /// of the random stream starting from this location. - /// - /// **RandBLAS functions do not mutate input RNGStates.** Free-functions - /// return new RNGStates with suitably updated counters. Constructors - /// for SketchingOperator objects store updated RNGStates in the - /// object's next_state member. - typename RNG::ctr_type counter; - - /// ------------------------------------------------------------------ - /// This is a Random123-defined statically-sized array of unsigned integers. - /// With RandBLAS' default, it contains two 32-bit unsigned ints - /// and is interpreted as one 64-bit unsigned int. - /// - /// This member specifices a sequece of pseudo-random numbers - /// that RNGState::generator can produce. Any fixed sequence has - /// fairly large period (\math{2^{132},} with RandBLAS' default) and - /// is statistically independent from sequences induced by different keys. - /// - /// To increment the key by "step," call \math{\ttt{key.incr(step)}}. - typename RNG::key_type key; - - const static int len_c = RNG::ctr_type::static_size; - static_assert(len_c >= 2); - const static int len_k = RNG::key_type::static_size; - - /// Initialize the counter and key to zero. - RNGState() : counter{}, key{} {} - - /// Initialize the counter and key to zero, then increment the key by k. - RNGState(uint64_t k) : counter{}, key{} { key.incr(k); } - - // construct from a key - RNGState(key_type const &k) : counter{}, key(k) {} - - // Initialize counter and key arrays at the given values. - RNGState(ctr_type const &c, key_type const &k) : counter(c), key(k) {} - - // move construct from an initial counter and key - RNGState(ctr_type &&c, key_type &&k) : counter(std::move(c)), key(std::move(k)) {} - - // move constructor. - RNGState(RNGState &&s) : RNGState(std::move(s.counter), std::move(s.key)) {}; - - ~RNGState() {}; - - /// Copy constructor. - RNGState(const RNGState &s) : RNGState(s.counter, s.key) {}; - - // A copy-assignment operator. - RNGState &operator=(const RNGState &s) { - std::memcpy(this->counter.v, s.counter.v, this->len_c * sizeof(ctr_uint)); - std::memcpy(this->key.v, s.key.v, this->len_k * sizeof(key_uint)); - return *this; - }; - - // - // Comparators (for now, these are just for testing and debugging) - // - - bool operator==(const RNGState &s) const { - // the compiler should only allow comparisons between RNGStates of the same type. - for (int i = 0; i < len_c; ++i) { - if (counter.v[i] != s.counter.v[i]) { return false; } - } - for (int i = 0; i < len_k; ++i) { - if (key.v[i] != s.key.v[i]) { return false; } - } - return true; - }; - - bool operator!=(const RNGState &s) const { - return !(*this == s); - }; - -}; - -template -const int RandBLAS::RNGState::len_c; - -template -const int RandBLAS::RNGState::len_k; - -template +template + requires requires(RNGState const& state, std::ostream& stream) { + state.counter.size(); + state.key.size(); + state.counter[0]; + state.key[0]; + stream << state.counter[0]; + stream << state.key[0]; + } std::ostream &operator<<( std::ostream &out, - const RNGState &s + const RNGState &s ) { - int i; + auto const& counter = s.counter; + auto const& key = s.key; out << "counter : {"; - for (i = 0; i < s.len_c - 1; ++i) { - out << s.counter[i] << ", "; + for (std::size_t i = 0; i + 1 < counter.size(); ++i) { + out << counter[i] << ", "; } - out << s.counter[i] << "}\n"; + out << counter[counter.size() - 1] << "}\n"; out << "key : {"; - for (i = 0; i < s.len_k - 1; ++i) { - out << s.key[i] << ", "; + for (std::size_t i = 0; i + 1 < key.size(); ++i) { + out << key[i] << ", "; } - out << s.key[i] << "}"; + out << key[key.size() - 1] << "}"; return out; } @@ -488,4 +378,3 @@ concept SketchingOperator = requires { #endif } // end namespace RandBLAS::base - diff --git a/RandBLAS/dense_skops.hh b/RandBLAS/dense_skops.hh index b89fa91e..82044bf8 100644 --- a/RandBLAS/dense_skops.hh +++ b/RandBLAS/dense_skops.hh @@ -62,9 +62,9 @@ inline void copy_promote(int n, const T_IN &a, T_OUT* b) { * "ptr" is the pointer offset for the desired submatrix in the imagined buffer of the parent matrix. * * @tparam T the data type of the matrix - * @tparam RNG a random123 CBRNG type + * @tparam state_t a counter-based RNG state type * @tparam OP an operator that transforms raw random values into matrix - * elements. See r123ext::uneg11 and r123ext::boxmul. + * elements. See rng::uneg11 and rng::boxmul. * * @param[in] n_cols * The number of columns in the implicitly defined parent matrix. @@ -93,18 +93,25 @@ inline void copy_promote(int n, const T_IN &a, T_OUT* b) { * using OMP_NUM_THREADS. The sequence of values generated does not depend on the number of threads. * */ -template -static RNGState fill_dense_submat_impl(int64_t n_cols, T* smat, int64_t n_srows, int64_t n_scols, int64_t ptr, const RNGState &seed, int64_t lda = 0) { +template +static state_t fill_dense_submat_impl(int64_t n_cols, T* smat, + int64_t n_srows, int64_t n_scols, + int64_t ptr, state_t const& seed, + int64_t lda = 0) { if (lda <= 0) { lda = n_scols; } else { randblas_require(lda >= n_scols); } randblas_require(n_cols >= n_scols); - RNG rng; - using CTR_t = typename RNG::ctr_type; - using KEY_t = typename RNG::key_type; - const int64_t ctr_size = CTR_t::static_size; + using res_t = typename state_t::res_t; + using word_t = typename res_t::value_type; + constexpr int64_t ctr_size = std::tuple_size_v; + static_assert(ctr_size % 2 == 0, + "dense sampling requires an even RNG result length"); + static_assert(std::numeric_limits::digits == 32 || + std::numeric_limits::digits == 64, + "dense sampling requires 32- or 64-bit result words"); int64_t pad = 0; // ^ computed such that n_cols+pad is divisible by ctr_size @@ -125,19 +132,17 @@ static RNGState fill_dense_submat_impl(int64_t n_cols, T* smat, int64_t n_s const bool one_block_per_row = ctr_mat_start == ctr_mat_row_end; const int64_t first_block_len = ((one_block_per_row) ? last_block_stop : ctr_size) - first_block_start; - CTR_t temp_c = seed.counter; - temp_c.incr(ctr_mat_start); - const CTR_t c = temp_c; - const KEY_t k = seed.key; + state_t first_state = seed; + first_state.advance(ctr_mat_start); #pragma omp parallel for schedule(static) for (int64_t row = 0; row < n_srows; row++) { int64_t incr_from_c = safe_int_product(ctr_inter_row_stride, row); - auto c_row = c; - c_row.incr(incr_from_c); - auto rv = OP::generate(rng, c_row, k); + state_t row_state = first_state; + row_state.advance(incr_from_c); + auto rv = OP::generate(row_state); T* smat_row = smat + row*lda; for (int i = 0; i < first_block_len; i++) { @@ -149,35 +154,34 @@ static RNGState fill_dense_submat_impl(int64_t n_cols, T* smat, int64_t n_s // middle blocks int64_t ind = first_block_len; for (int i = 0; i < (ctr_mat_row_end - ctr_mat_start - 1); ++i) { - c_row.incr(); - rv = OP::generate(rng, c_row, k); + row_state.advance(1); + rv = OP::generate(row_state); copy_promote(ctr_size, rv, smat_row + ind); ind = ind + ctr_size; } // last block - c_row.incr(); - rv = OP::generate(rng, c_row, k); + row_state.advance(1); + rv = OP::generate(row_state); copy_promote(last_block_stop, rv, smat_row + ind); } - // find the largest counter in the counter array - CTR_t max_c = c; - max_c.incr(n_srows * ctr_inter_row_stride); - return RNGState {max_c, k}; + state_t next_state = seed; + next_state.advance(ctr_mat_start + n_srows * ctr_inter_row_stride); + return next_state; } -template -RNGState compute_next_state(DD dist, RNGState state) { +template +state_t compute_next_state(DD dist, state_t state) { int64_t major_len = dist.dim_major; int64_t minor_len = dist.dim_minor; - int64_t ctr_size = RNG::ctr_type::static_size; + constexpr int64_t ctr_size = std::tuple_size_v; int64_t pad = 0; if (major_len % ctr_size != 0) { pad = ctr_size - major_len % ctr_size; } int64_t ctr_major_axis_stride = (major_len + pad) / ctr_size; int64_t full_incr = safe_int_product(ctr_major_axis_stride, minor_len); - state.counter.incr(full_incr); + state.advance(full_incr); return state; } @@ -203,7 +207,7 @@ namespace RandBLAS { // Forward declaration of DenseSkOp. It's returnable by // DenseDist.sample(), but its definition involves DenseDist. -template +template struct DenseSkOp; @@ -327,8 +331,8 @@ struct DenseDist { // ------------------------------------------------------------------------------------- /// Construct a DenseSkOp with this distribution and the provided seed_state. - template - DenseSkOp sample(RNGState &seed_state) { + template + DenseSkOp sample(state_t &seed_state) { return {*this, seed_state}; } @@ -351,7 +355,7 @@ struct DenseDist { /// A sample from a distribution over matrices whose entries are iid /// mean-zero variance-one random variables. /// This type conforms to the SketchingOperator concept. -template +template struct DenseSkOp { // --------------------------------------------------------------------------- @@ -360,7 +364,7 @@ struct DenseSkOp { // --------------------------------------------------------------------------- /// Type alias. - using state_t = RNGState; + using state_t = generator_state_t; // --------------------------------------------------------------------------- /// Real scalar type used in matrix representations of this operator. @@ -451,7 +455,7 @@ struct DenseSkOp { // Move constructor DenseSkOp( - DenseSkOp &&S + DenseSkOp &&S ) : // Initializations dist(S.dist), seed_state(S.seed_state), @@ -557,8 +561,10 @@ static_assert(SketchingOperator>); /// - Used to define :math:`\mtxS` as a sample from :math:`\D.` /// /// @endverbatim -template -RNGState fill_dense_unpacked(blas::Layout layout, const DenseDist &D, int64_t n_rows, int64_t n_cols, int64_t ro_s, int64_t co_s, T* buff, const RNGState &seed) { +template +state_t fill_dense_unpacked(blas::Layout layout, const DenseDist &D, + int64_t n_rows, int64_t n_cols, int64_t ro_s, + int64_t co_s, T* buff, state_t const& seed) { using RandBLAS::dense::fill_dense_submat_impl; randblas_require(D.n_rows >= n_rows + ro_s); randblas_require(D.n_cols >= n_cols + co_s); @@ -575,14 +581,16 @@ RNGState fill_dense_unpacked(blas::Layout layout, const DenseDist &D, int64 n_cols_ = n_cols; ptr = safe_int_product(ro_s, ma_len) + co_s; } - RNGState next_state{}; + state_t next_state{}; switch (D.family) { case ScalarDist::Gaussian: { - next_state = fill_dense_submat_impl(ma_len, buff, n_rows_, n_cols_, ptr, seed); + next_state = fill_dense_submat_impl( + ma_len, buff, n_rows_, n_cols_, ptr, seed); break; } case ScalarDist::Uniform: { - next_state = fill_dense_submat_impl(ma_len, buff, n_rows_, n_cols_, ptr, seed); + next_state = fill_dense_submat_impl( + ma_len, buff, n_rows_, n_cols_, ptr, seed); blas::scal(n_rows_ * n_cols_, (T)std::sqrt(3), buff, 1); break; } @@ -617,8 +625,8 @@ RNGState fill_dense_unpacked(blas::Layout layout, const DenseDist &D, int64 /// A CBRNG state /// - Used to define \math{\mat(\buff)} as a sample from \math{\D}. /// -template -RNGState fill_dense(const DenseDist &D, T *buff, const RNGState &seed) { +template +state_t fill_dense(const DenseDist &D, T *buff, state_t const& seed) { return fill_dense_unpacked(D.natural_layout, D, D.n_rows, D.n_cols, 0, 0, buff, seed); } diff --git a/RandBLAS/random_gen.hh b/RandBLAS/random_gen.hh index 92980703..807f951a 100644 --- a/RandBLAS/random_gen.hh +++ b/RandBLAS/random_gen.hh @@ -32,138 +32,57 @@ /// @file #include "compilers.hh" -#include -#include -#include -#include -// NOTE: we do not support Random123's AES or ARS generators. - -// Host-side shim for sincospi / sincospif. Random123/boxmuller.hpp calls these -// unqualified; its own fallback definitions are gated behind -// `!defined(CUDART_VERSION) || CUDART_VERSION < 5000` -// and are therefore skipped on host compiles where any dependency leaks -// into the include chain (notably blaspp's CMake config -// transitively exports CUDA_INCLUDE_DIRS when blaspp was built with CUDA -// support, which then poisons RandBLAS's host test compiles). Without these -// definitions the host-compiled code fails to link. -// -// An earlier version of this shim was removed in PR #156 (9b73a25, 2026-03-05) -// on the theory that CI green = unused. RandBLAS CI does not currently -// exercise the Linux + CUDA-aware blaspp configuration where the failure -// manifests; restoring the shim so host builds succeed on that config. -#if !defined(__CUDACC__) -#include -static inline void sincospif(float x, float *s, float *c) { - const float PIf = 3.1415926535897932f; - *s = std::sin(PIf * x); - *c = std::cos(PIf * x); -} -static inline void sincospi(double x, double *s, double *c) { - const double PI = 3.1415926535897932; - *s = std::sin(PI * x); - *c = std::cos(PI * x); -} -#endif - -#include -#include - -/// our extensions to random123 -namespace r123ext -{ -/** Apply boxmuller transform to all elements of ri. The number of elements of r - * must be evenly divisible by 2. See also r123::uneg11all. - * - * @tparam CTR a random123 CBRNG ctr_type - * @tparam T the return element type. The default return type is dictated by - * the RNG's ctr_type's value_type : float for 32 bit counter elements - * and double for 64. - * - * @param[in] ri a sequence of N random values generated using random123 CBRNG - * type RNG. The transform is applied pair wise to the sequence. - * - * @returns a std::array of transformed floating point values. - */ -template ::type> -auto boxmulall( - CTR const &ri -) { - std::array ro; - int nit = CTR::static_size / 2; - for (int i = 0; i < nit; ++i) - { - auto [v0, v1] = r123::boxmuller(ri[2*i], ri[2*i + 1]); - ro[2*i ] = v0; - ro[2*i + 1] = v1; - } - return ro; -} - -/** @defgroup generators - * Generators take CBRNG, counter,and key instances and return a sequence of - * random floating point numbers in a std::array. The length of the squence is - * the length of the counter and the precision is float for 32 bit counters and - * double for 64. - */ -/// @{ - -/// Generate a sequence of random values and apply a Box-Muller transform. -struct boxmul -{ - /** Generate a sequence of random values and apply a Box-Muller transform. - * - * @tparam RNG a random123 CBRNG type - * - * @param[in] rng: a random123 CBRNG instance used to generate the sequence - * @param[in] c: the CBRNG counter - * @param[in] k: the CBRNG key - * - * @returns a std::array where N is the CBRNG's ctr_type::static_size - * and T is deduced from the RNG's counter element type : float - * for 32 bit counter elements and double for 64. For example when - * RNG is Philox4x32 the return is a std::array. - */ - template - static - auto generate( - RNG &rng, - typename RNG::ctr_type const &c, - typename RNG::key_type const &k - ) { - return boxmulall(rng(c,k)); +#include "rng/concepts.hh" +#include "rng/distributions.hh" +#include "rng/philox.hh" +#include "rng/repacked_output.hh" +#include "rng/word_array.hh" + +#include +#include +#include + +namespace RandBLAS { + +using DefaultRNG = rng::Philox<4, 32, 10>; + +/// Copyable state that binds an engine to one counter and one key. +template +struct RNGState { + using engine_t = Engine; + using ctr_t = typename Engine::ctr_t; + using key_t = typename Engine::key_t; + using res_t = typename Engine::res_t; + + ctr_t counter{}; + key_t key{}; + [[no_unique_address]] Engine engine{}; + + constexpr RNGState() = default; + + constexpr RNGState(std::uint64_t seed) noexcept( + noexcept(Engine::make_key(seed))) + requires rng::SeedMappableEngine + : key(Engine::make_key(seed)) {} + + explicit constexpr RNGState(key_t const& input_key) : key(input_key) {} + + constexpr RNGState(ctr_t const& input_counter, key_t const& input_key) + : counter(input_counter), key(input_key) {} + + constexpr void generate(res_t& output) const noexcept( + noexcept(engine.generate(counter, key, output))) { + engine.generate(counter, key, output); } -}; -/// Generate a sequence of random values and transform to -1.0 to 1.0. -struct uneg11 -{ - /** Generate a sequence of random values and transform to -1.0 to 1.0. - * - * @tparam RNG a random123 CBRNG type - * - * @param[in] rng: a random123 CBRNG instance used to generate the sequence - * @param[in] c: CBRNG counter - * @param[in] k: CBRNG key - * - * @returns a std::array where N is the CBRNG's ctr_type::static_size - * and T is deduced from the RNG's counter element type : float - * for 32 bit counter elements and double for 64. For example when - * RNG is Philox4x32 the return is a std::array. - */ - template ::type> - static - auto generate( - RNG &rng, - typename RNG::ctr_type const &c, - typename RNG::key_type const &k - ) { - return r123::uneg11all(rng(c,k)); + constexpr void advance(std::uint64_t blocks) noexcept( + noexcept(counter.advance(blocks))) { + counter.advance(blocks); } -}; -/// @} + friend constexpr bool operator==(RNGState const&, RNGState const&) = default; +}; -} // end of namespace r123ext +using DefaultRNGState = RNGState; +} // namespace RandBLAS diff --git a/RandBLAS/rng/DevNotes.md b/RandBLAS/rng/DevNotes.md new file mode 100644 index 00000000..11772109 --- /dev/null +++ b/RandBLAS/rng/DevNotes.md @@ -0,0 +1,168 @@ +# Random-number generation developer notes + +RandBLAS provides native, header-only counter-based random-number generation. +This document records its public contracts, reproducibility guarantees, and +validation requirements. + +## Public engine and state contracts + +A counter-based engine is a stateless block function. It owns the types and +meaning of its counter, key, and result: + +```cpp +void generate(ctr_t const& counter, + key_t const& key, + res_t& output) const; +``` + +The third argument is output-only. The engine writes every lane and does not +mutate the counter or key. The engine contract is structural; engines do not +inherit from a RandBLAS base class. + +RandBLAS algorithms consume state-like objects instead of engines directly. +The structural `GeneratorState` concept requires a fixed-size unsigned `res_t` +and these operations: + +```cpp +void generate(res_t& output) const; +void advance(std::uint64_t blocks); +``` + +`RNGState` is the provided transparent adapter. Its public `counter`, +`key`, and `engine` values represent the complete state. Generic algorithms +depend only on the `GeneratorState` operations and do not require or inspect +those public members. + +## Relationship to the C++ standard random facilities + +| RandBLAS abstraction | Nearest standard-library abstraction | Deliberate difference | +|---|---|---| +| `rng::WordArray` | `std::array` | Adds little-endian, extended-width modular `advance(uint64_t)`; it is not a generator. | +| `rng::Philox` | C++26 `std::philox_engine` | RandBLAS is a stateless `(counter, key) -> block` function. The standard engine owns state, caches a block position, and returns one scalar per mutating `operator()`. | +| `RNGState` | State stored inside a standard random-number engine | Provides transparent counter, key, and engine values together with nonmutating block generation and explicit block advancement. It has no scalar `operator()`, serialization, seed sequence, or cached lane index. | +| `rng::RepackedOutput` | `std::independent_bits_engine` | Re-expresses every bit of one existing block in fixed LSB-first chunks. It does not draw a variable number of scalar values or define a new stream position. | +| `rng::u01`, `rng::boxmuller`, `rng::uneg11`, `rng::boxmul` | `std::uniform_real_distribution` and `std::normal_distribution` | Transforms explicit words or result blocks without owning or mutating generator state. | +| `GeneratorState` | `std::uniform_random_bit_generator` | Produces a fixed result block without mutation; a URBG produces one scalar by mutating itself. Neither native RandBLAS concept models URBG. | + +Neither a RandBLAS engine nor state is a standard uniform random bit generator. +A scalar standard-engine adaptor can be designed independently if one is ever +needed. + +## Counter and seed semantics + +Each engine chooses its `ctr_t`, including the counter's period and the meaning +of `advance(1)`. The counter type implements modular `advance(uint64_t)`; +generic state and sampler code delegates to that operation. + +Native Philox uses `WordArray` for `ctr_t`. Lane zero is the least +significant word of an `N * W`-bit unsigned integer, and `advance(k)` adds `k` +to that integer modulo `2^(N * W)`. One call to `generate` produces exactly one +`N`-word block at the current counter without mutating the counter or key. + +An engine may provide `static key_t make_key(uint64_t)`. Only engines with that +hook support `RNGState(uint64_t)`. Explicit key and counter/key construction +remains available for known-answer tests and expert use. Native Philox preserves +the existing scalar-seed interpretation: start with a zero key and add the +64-bit seed using the key's extended-width unsigned representation. + +## Output blocks and repacking order + +`res_t` is a fixed-size array of unsigned words. `RepackedOutput` preserves all +bits and the wrapped engine's block boundary while exposing narrower unsigned +lanes. Source-word order is preserved, and chunks within a word are emitted +least-significant first. Thus `0xAABBCCDD` becomes `{0xCCDD, 0xAABB}` for +16-bit output and `{0xDD, 0xCC, 0xBB, 0xAA}` for 8-bit output, independently of +host byte order. + +Current samplers continue to consume native 32- or 64-bit lanes. This migration +does not define bit assembly for 8- or 16-bit sampler inputs. + +## Coordinate-addressed sampling + +Sampling is indexed by coordinates and reserved counter blocks, not by the +execution order of a shared stream. A worker copies the input state, advances +that copy by the block offset for its region, generates local blocks, and leaves +the caller's state unchanged. The returned state is a copy advanced by the +total reservation. + +This mapping is what makes full/submatrix generation consistent and makes +generated operators independent of OpenMP thread count. Dense row padding and +sparse block reservations are compatibility constraints during the migration. + +`RandBLAS::testing::detail::RNGStream` is test infrastructure only. It adapts +result blocks to sequential scalar draws for random sparse test-matrix +generation. See `test/DevNotes.md` for its consumption details. + +## Floating-point reproducibility + +Native transforms retain the integer-to-floating formulas, constants, +endpoints, word assignment, and output order used by RandBLAS through +Random123. They use `std::sin`, `std::cos`, `std::log`, and `std::sqrt` on the +host. Dense Gaussian values may therefore differ in the last bits across math +libraries, compilers, and architectures. The integer Philox stream and default +sparse operator output remain bitwise compatibility requirements. + +The supported transform names are `u01`, `boxmuller`, `uneg11`, and `boxmul`. + +## Algorithm provenance and licensing + +The Philox algorithm is described by Salmon, Moraes, Dror, and Shaw in +[Parallel Random Numbers: As Easy as 1, 2, 3](https://doi.org/10.1145/2063384.2063405). +The implementation, floating-point transformations, and known-answer material +are adapted from D. E. Shaw Research's Random123 project. Files containing +adapted implementation or test material retain the applicable D. E. Shaw +Research BSD-3-Clause notice. Static vectors were generated from Random123 +commit [`9545ff6413f258be2f04c1d319d99aaef7521150`](https://github.com/DEShawResearch/random123/commit/9545ff6413f258be2f04c1d319d99aaef7521150). + +Native Philox is a statistical counter-based generator, not a cryptographic +random-number generator. + +The default engine is `rng::Philox<4, 32, 10>`. Its integer blocks are bitwise +identical to Random123 Philox4x32-10 for the same counter and key. This guarantee +also covers the default sparse sketch stream. Dense Gaussian sampling preserves +the same formulas and block assignment, subject to the host-math limitation +described above. + +## Known-answer and statistical testing + +Static known-answer vectors cover each supported Philox word count, word width, +and round count. Those vectors are generated once from the pinned Random123 +commit `9545ff6413f258be2f04c1d319d99aaef7521150`; normal builds and tests never +locate Random123. Separate tests cover +counter carries and wraparound, engine/state concepts, seed mapping, output +repacking, floating-point endpoints, Box--Muller reference values, statistical +behavior, sampler state advancement, full/submatrix agreement, and OpenMP +thread-count independence. + +Characterization fixtures captured before the migration protect the default +dense and sparse streams. Installed-package and example builds protect the +absence of a transitive Random123 dependency. + +Dense samplers require 32- or 64-bit result words and an even number of lanes. +`sample_indices_iid_uniform` requires at least two 32-bit lanes (three when +also producing Rademacher signs), or at least one 64-bit lane (two with signs). +Sparse Fisher--Yates sampling requires at least three 32-bit lanes or two +64-bit lanes. These constraints are sampler contracts, not requirements of the +base engine and state concepts. In particular, current samplers do not accept +8- or 16-bit `RepackedOutput` results. + +## Adding another engine + +A new engine supplies value-semantic `ctr_t`, `key_t`, and fixed unsigned +`res_t` types plus the output-only `generate` function. Its counter supplies +`advance(uint64_t)`. It may supply `make_key(uint64_t)` if scalar seeding has a +stable, documented mapping. No inheritance or Philox-specific representation is +required. + +The engine should have direct known-answer tests, counter-period tests, seed +mapping tests, and appropriate statistical validation before it is used by a +sampler. Sampler-specific result width and lane-count requirements remain +separate from the base engine concept. + +## Performance validation + +The migration records pre-change and native timings for the basic dense RNG +benchmark and sparse sketch sampling benchmark under the same compiler, build +type, dimensions, and thread count. A regression outside ordinary run-to-run +variation must be investigated before merge. Optimization beyond parity is not +part of the dependency-removal change. diff --git a/RandBLAS/rng/concepts.hh b/RandBLAS/rng/concepts.hh new file mode 100644 index 00000000..dc4fa81a --- /dev/null +++ b/RandBLAS/rng/concepts.hh @@ -0,0 +1,87 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#pragma once + +#include +#include +#include + +namespace RandBLAS::rng { + +namespace detail { + +template +concept FixedUnsignedBlock = requires { + typename block_t::value_type; + requires std::unsigned_integral; + requires(std::tuple_size_v > 0); +}; + +} // namespace detail + +/// Stateless counter-based engine producing one fixed-size result block. +template +concept CounterBasedEngine = + std::semiregular && requires { + typename engine_t::ctr_t; + typename engine_t::key_t; + typename engine_t::res_t; + requires std::regular; + requires std::regular; + requires detail::FixedUnsignedBlock; + } && requires(engine_t const& engine, typename engine_t::ctr_t& counter, + typename engine_t::ctr_t const& const_counter, + typename engine_t::key_t const& key, + typename engine_t::res_t& output, std::uint64_t blocks) { + { counter.advance(blocks) } -> std::same_as; + { engine.generate(const_counter, key, output) } -> std::same_as; + }; + +template +concept SeedMappableEngine = + CounterBasedEngine && requires(std::uint64_t seed) { + { engine_t::make_key(seed) } -> std::same_as; + }; + +} // namespace RandBLAS::rng + +namespace RandBLAS { + +/// Copyable generator state that produces and advances fixed-size blocks. +template +concept GeneratorState = std::copyable && requires { + typename state_t::res_t; + requires rng::detail::FixedUnsignedBlock; +} && requires(state_t& state, state_t const& const_state, + typename state_t::res_t& output, std::uint64_t blocks) { + { const_state.generate(output) } -> std::same_as; + { state.advance(blocks) } -> std::same_as; +}; + +} // namespace RandBLAS diff --git a/RandBLAS/rng/distributions.hh b/RandBLAS/rng/distributions.hh new file mode 100644 index 00000000..628e84a5 --- /dev/null +++ b/RandBLAS/rng/distributions.hh @@ -0,0 +1,137 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +/* +Copyright 2010-2011, D. E. Shaw Research. +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 D. E. Shaw Research nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. +*/ + +// The integer conversions and Box--Muller mapping in this file are adapted +// from Random123's uniform.hpp and boxmuller.hpp. + +#pragma once + +#include "concepts.hh" + +#include +#include +#include +#include +#include +#include + +namespace RandBLAS::rng { + +namespace detail { + +template +using default_real_t = std::conditional_t; + +} // namespace detail + +/// Convert a random unsigned word to a floating-point value in (0, 1]. +template +[[nodiscard]] constexpr real_t u01(word_t input) noexcept { + static_assert(std::is_unsigned_v); + static_assert(sizeof(word_t) == 4 || sizeof(word_t) == 8); + static_assert(std::is_same_v || std::is_same_v); + constexpr real_t factor = + real_t{1} / + (static_cast(std::numeric_limits::max()) + real_t{1}); + constexpr real_t half_factor = real_t{0.5} * factor; + return static_cast(input) * factor + half_factor; +} + +/// Symmetric-uniform conversion and dense-sampling transform policy. +struct uneg11 { + /// Convert a random unsigned word to a floating-point value in [-1, 1]. + template + [[nodiscard]] static constexpr real_t convert(word_t input) noexcept { + static_assert(std::is_unsigned_v); + static_assert(sizeof(word_t) == 4 || sizeof(word_t) == 8); + static_assert(std::is_same_v || std::is_same_v); + using signed_word_t = std::make_signed_t; + constexpr real_t factor = + real_t{1} / + (static_cast( + std::numeric_limits::max()) + real_t{1}); + constexpr real_t half_factor = real_t{0.5} * factor; + return static_cast(static_cast(input)) * factor + + half_factor; + } + + template + [[nodiscard]] static auto generate(state_t const& state) { + using bits_t = typename state_t::res_t; + using word_t = typename bits_t::value_type; + using real_t = detail::default_real_t; + constexpr std::size_t count = std::tuple_size_v; + bits_t bits{}; + std::array output{}; + state.generate(bits); + for (std::size_t i = 0; i < count; ++i) { + output[i] = convert(bits[i]); + } + return output; + } +}; + +/// Transform an angle word and a radius word into sine-then-cosine normals. +template +[[nodiscard]] inline auto boxmuller(word_t angle_word, word_t radius_word) { + static_assert(std::is_unsigned_v); + static_assert(sizeof(word_t) == 4 || sizeof(word_t) == 8); + using real_t = detail::default_real_t; + constexpr real_t pi = real_t{3.1415926535897932}; + auto angle = pi * uneg11::convert(angle_word); + auto radius = std::sqrt(real_t{-2} * std::log(u01(radius_word))); + return std::array{std::sin(angle) * radius, + std::cos(angle) * radius}; +} + +/// Box--Muller dense-sampling transform policy. +struct boxmul { + template + requires(std::tuple_size_v % 2 == 0) + [[nodiscard]] static auto generate(state_t const& state) { + using bits_t = typename state_t::res_t; + using word_t = typename bits_t::value_type; + using real_t = detail::default_real_t; + constexpr std::size_t count = std::tuple_size_v; + bits_t bits{}; + std::array output{}; + state.generate(bits); + for (std::size_t i = 0; i < count; i += 2) { + auto pair = boxmuller(bits[i], bits[i + 1]); + output[i] = pair[0]; + output[i + 1] = pair[1]; + } + return output; + } +}; + +} // namespace RandBLAS::rng diff --git a/RandBLAS/rng/philox.hh b/RandBLAS/rng/philox.hh new file mode 100644 index 00000000..ca84bed7 --- /dev/null +++ b/RandBLAS/rng/philox.hh @@ -0,0 +1,207 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +/* +Copyright 2010-2011, D. E. Shaw Research. +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 D. E. Shaw Research nor the names of its contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. +*/ + +// Adapted from Random123's Philox implementation. The public interface and +// storage types are native RandBLAS facilities; constants and round operations +// intentionally preserve Random123's bit stream. + +#pragma once + +#include "word_array.hh" + +#include +#include +#include +#include +#include + +#if defined(_MSC_VER) && defined(_M_X64) +#include +#endif + +namespace RandBLAS::rng { + +namespace detail { + +template +[[nodiscard]] constexpr Word mulhilo(Word left, Word right, Word* high) noexcept { + static_assert(sizeof(Word) == 4 || sizeof(Word) == 8); + + if constexpr (sizeof(Word) == 4) { + auto product = static_cast(left) * + static_cast(right); + *high = static_cast(product >> 32); + return static_cast(product); + } else { +#if defined(_MSC_VER) && defined(_M_X64) + unsigned __int64 native_high; + auto low = _umul128(static_cast(left), + static_cast(right), + &native_high); + *high = static_cast(native_high); + return static_cast(low); +#elif defined(__SIZEOF_INT128__) + using double_word_t = unsigned __int128; + auto product = static_cast(left) * + static_cast(right); + *high = static_cast(product >> 64); + return static_cast(product); +#else + static_assert(sizeof(Word) != 8, + "64-bit Philox requires unsigned __int128 or _umul128"); +#endif + } +} + +template +struct PhiloxConstants { + using word_t = std::conditional_t; + + static constexpr word_t multiplier_0 = [] { + if constexpr (W == 32 && N == 2) { + return UINT32_C(0xd256d193); + } else if constexpr (W == 32) { + return UINT32_C(0xd2511f53); + } else if constexpr (N == 2) { + return UINT64_C(0xd2b74407b1ce6e93); + } else { + return UINT64_C(0xd2e7470ee14c6c93); + } + }(); + + static constexpr word_t multiplier_1 = [] { + if constexpr (W == 32) { + return UINT32_C(0xcd9e8d57); + } else { + return UINT64_C(0xca5a826395121157); + } + }(); + + static constexpr word_t weyl_0 = [] { + if constexpr (W == 32) { + return UINT32_C(0x9e3779b9); + } else { + return UINT64_C(0x9e3779b97f4a7c15); + } + }(); + + static constexpr word_t weyl_1 = [] { + if constexpr (W == 32) { + return UINT32_C(0xbb67ae85); + } else { + return UINT64_C(0xbb67ae8584caa73b); + } + }(); +}; + +template +constexpr void apply_philox_round(std::array& block, + key_t const& key) noexcept { + using constants_t = PhiloxConstants; + + auto input_0 = block[0]; + auto input_1 = block[1]; + word_t high_0; + auto low_0 = mulhilo(constants_t::multiplier_0, input_0, &high_0); + + if constexpr (N == 2) { + block[0] = static_cast(high_0 ^ key[0] ^ input_1); + block[1] = low_0; + } else { + auto input_2 = block[2]; + auto input_3 = block[3]; + word_t high_1; + auto low_1 = mulhilo(constants_t::multiplier_1, input_2, &high_1); + block[0] = static_cast(high_1 ^ input_1 ^ key[0]); + block[1] = low_1; + block[2] = static_cast(high_0 ^ input_3 ^ key[1]); + block[3] = low_0; + } +} + +template +constexpr void bump_philox_key(key_t& key) noexcept { + using constants_t = PhiloxConstants; + using word_t = typename key_t::value_type; + + key[0] = static_cast(key[0] + constants_t::weyl_0); + if constexpr (N == 4) { + key[1] = static_cast(key[1] + constants_t::weyl_1); + } +} + +} // namespace detail + +/// Stateless Philox counter-based random-number engine. +/// +/// Word zero is the least-significant word of counters and keys. `generate` +/// maps one counter/key pair to one result block without modifying its inputs. +template +struct Philox { + static_assert(N == 2 || N == 4, "Philox supports two or four words"); + static_assert(W == 32 || W == 64, "Philox supports 32- or 64-bit words"); + static_assert(R <= 16, "Philox supports at most 16 rounds"); + + using word_t = std::conditional_t; + using ctr_t = WordArray; + using key_t = WordArray; + using res_t = std::array; + + /// Map a scalar seed to a key by adding it to an all-zero key. + [[nodiscard]] static constexpr key_t make_key(std::uint64_t seed) noexcept { + key_t key{}; + key.advance(seed); + return key; + } + + /// Generate a complete result block, overwriting every output lane. + constexpr void generate(ctr_t const& counter, key_t const& key, + res_t& output) const noexcept { + res_t block{}; + for (std::size_t i = 0; i < N; ++i) { + block[i] = counter[i]; + } + + key_t round_key = key; + for (std::size_t round = 0; round < R; ++round) { + detail::apply_philox_round(block, round_key); + if (round + 1 < R) { + detail::bump_philox_key(round_key); + } + } + output = block; + } + + friend constexpr bool operator==(Philox const&, Philox const&) = default; +}; + +} // namespace RandBLAS::rng diff --git a/RandBLAS/rng/repacked_output.hh b/RandBLAS/rng/repacked_output.hh new file mode 100644 index 00000000..849b05db --- /dev/null +++ b/RandBLAS/rng/repacked_output.hh @@ -0,0 +1,135 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#pragma once + +#include "concepts.hh" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace RandBLAS::rng { + +namespace detail { + +template +inline constexpr bool valid_repacking_widths = + OutputBits > 0 && OutputBits <= SourceBits && + SourceBits % OutputBits == 0 && + std::has_single_bit(SourceBits / OutputBits); + +} // namespace detail + +template +concept ValidRepacking = + std::unsigned_integral && + std::unsigned_integral && + detail::valid_repacking_widths::digits, + std::numeric_limits::digits>; + +/// Express each result word of an engine as fixed-width, LSB-first chunks. +/// +/// The adaptor preserves the wrapped engine's counter, key, seed mapping, and +/// total number of bits per block. Equal-width adaptation is an identity. +template + requires CounterBasedEngine && + ValidRepacking +struct RepackedOutput { + using source_res_t = typename Engine::res_t; + using source_word_t = typename source_res_t::value_type; + + static constexpr std::size_t source_word_bits = + std::numeric_limits::digits; + static constexpr std::size_t output_word_bits = + std::numeric_limits::digits; + static constexpr std::size_t chunks_per_source_word = + source_word_bits / output_word_bits; + static constexpr std::size_t source_word_count = + std::tuple_size_v; + + using word_t = OutputWord; + using ctr_t = typename Engine::ctr_t; + using key_t = typename Engine::key_t; + static constexpr std::size_t repacked_word_count = + source_word_count * chunks_per_source_word; + using res_t = std::array; + + constexpr RepackedOutput() requires std::default_initializable = default; + + constexpr explicit RepackedOutput(Engine engine) noexcept( + std::is_nothrow_move_constructible_v) + : engine(std::move(engine)) {} + + [[nodiscard]] static constexpr key_t make_key(std::uint64_t seed) noexcept( + noexcept(Engine::make_key(seed))) + requires requires { + { Engine::make_key(seed) } -> std::same_as; + } + { + return Engine::make_key(seed); + } + + constexpr void generate(ctr_t const& counter, key_t const& key, + res_t& output) const noexcept( + noexcept(engine.generate(counter, key, + std::declval()))) { + source_res_t source{}; + engine.generate(counter, key, source); + + constexpr source_word_t mask = [] { + if constexpr (output_word_bits == source_word_bits) { + return std::numeric_limits::max(); + } else { + return static_cast( + (source_word_t{1} << output_word_bits) - 1); + } + }(); + + std::size_t output_index = 0; + for (source_word_t source_word : source) { + for (std::size_t chunk = 0; chunk < chunks_per_source_word; ++chunk) { + auto shifted = static_cast( + source_word >> (chunk * output_word_bits)); + output[output_index++] = static_cast(shifted & mask); + } + } + } + + friend constexpr bool operator==(RepackedOutput const&, RepackedOutput const&) = default; + + [[no_unique_address]] Engine engine{}; +}; + +} // namespace RandBLAS::rng diff --git a/RandBLAS/rng/word_array.hh b/RandBLAS/rng/word_array.hh new file mode 100644 index 00000000..9b894362 --- /dev/null +++ b/RandBLAS/rng/word_array.hh @@ -0,0 +1,95 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#pragma once + +#include +#include +#include +#include +#include + +namespace RandBLAS::rng { + +/// Fixed-size little-endian words interpreted as one modular unsigned value. +template +struct WordArray { + static_assert(WordCount > 0, "WordArray requires at least one word"); + + using value_type = Word; + static constexpr std::size_t static_size = WordCount; + + std::array words{}; + + constexpr Word& operator[](std::size_t i) noexcept { + return words[i]; + } + + constexpr Word const& operator[](std::size_t i) const noexcept { + return words[i]; + } + + [[nodiscard]] static constexpr std::size_t size() noexcept { + return WordCount; + } + + /// Add an unsigned 64-bit amount, discarding overflow past the last word. + constexpr void advance(std::uint64_t amount) noexcept { + constexpr auto word_bits = std::numeric_limits::digits; + std::uint64_t remaining = amount; + Word carry = 0; + + for (std::size_t i = 0; i < WordCount; ++i) { + Word addend; + if constexpr (word_bits < 64) { + constexpr std::uint64_t mask = + (std::uint64_t{1} << word_bits) - 1; + addend = static_cast(remaining & mask); + remaining >>= word_bits; + } else { + addend = static_cast(remaining); + remaining = 0; + } + + Word after_addend = static_cast(words[i] + addend); + bool addend_overflow = after_addend < words[i]; + Word after_carry = static_cast(after_addend + carry); + bool carry_overflow = after_carry < after_addend; + words[i] = after_carry; + carry = static_cast(addend_overflow || carry_overflow); + + if (remaining == 0 && carry == 0) { + break; + } + } + } + + friend constexpr bool operator==(WordArray const&, WordArray const&) = default; +}; + +} // namespace RandBLAS::rng diff --git a/RandBLAS/skge.hh b/RandBLAS/skge.hh index aab9f58c..75fccef4 100644 --- a/RandBLAS/skge.hh +++ b/RandBLAS/skge.hh @@ -535,7 +535,7 @@ void _rskges_compress_and_apply_coo( /// - Leading dimension of \math{\mat(B)} when reading from \math{B}. /// - Refer to documentation for \math{\lda} for details. /// -template +template void lskges( blas::Layout layout, blas::Op opS, @@ -544,7 +544,7 @@ void lskges( int64_t n, // \op(A) is m-by-n int64_t m, // \op(submat(S)) is d-by-m T alpha, - const SparseSkOp &S, + const SparseSkOp &S, int64_t ro_s, int64_t co_s, const T *A, @@ -674,7 +674,7 @@ void lskges( /// - Leading dimension of \math{\mat(B)} when reading from \math{B}. /// - Refer to documentation for \math{\lda} for details. /// -template +template inline void rskges( blas::Layout layout, blas::Op opA, @@ -685,7 +685,7 @@ inline void rskges( T alpha, const T *A, int64_t lda, - const SparseSkOp &S, + const SparseSkOp &S, int64_t ro_s, int64_t co_s, T beta, @@ -856,7 +856,7 @@ inline void sketch_general( int64_t ldb ); -template +template inline void sketch_general( blas::Layout layout, blas::Op opS, @@ -865,7 +865,7 @@ inline void sketch_general( int64_t n, // op(A) is m-by-n int64_t m, // op(submat(\mtxS)) is d-by-m T alpha, - const SparseSkOp &S, + const SparseSkOp &S, int64_t ro_s, int64_t co_s, const T *A, @@ -880,7 +880,7 @@ inline void sketch_general( ); } -template +template inline void sketch_general( blas::Layout layout, blas::Op opS, @@ -889,7 +889,7 @@ inline void sketch_general( int64_t n, // op(A) is m-by-n int64_t m, // op(submat(\mtxS)) is d-by-m T alpha, - const DenseSkOp &S, + const DenseSkOp &S, int64_t ro_s, int64_t co_s, const T *A, @@ -1028,7 +1028,7 @@ inline void sketch_general( int64_t ldb ); -template +template inline void sketch_general( blas::Layout layout, blas::Op opA, @@ -1039,7 +1039,7 @@ inline void sketch_general( T alpha, const T *A, int64_t lda, - const DenseSkOp &S, + const DenseSkOp &S, int64_t ro_s, int64_t co_s, T beta, @@ -1052,7 +1052,7 @@ inline void sketch_general( } -template +template inline void sketch_general( blas::Layout layout, blas::Op opA, @@ -1063,7 +1063,7 @@ inline void sketch_general( T alpha, const T *A, int64_t lda, - const SparseSkOp &S, + const SparseSkOp &S, int64_t ro_s, int64_t co_s, T beta, diff --git a/RandBLAS/sparse_data/base.hh b/RandBLAS/sparse_data/base.hh index ecf98471..f2c37b71 100644 --- a/RandBLAS/sparse_data/base.hh +++ b/RandBLAS/sparse_data/base.hh @@ -31,6 +31,7 @@ #include "RandBLAS/config.h" #include "RandBLAS/base.hh" #include +#include #include #ifdef __cpp_concepts @@ -58,7 +59,7 @@ int64_t nnz_in_dense(int64_t n_rows, int64_t n_cols, int64_t stride_row, int64_t int64_t nnz = 0; for (int64_t i = 0; i < n_rows; ++i) { for (int64_t j = 0; j < n_cols; ++j) { - if (abs(MAT(i, j)) > abs_tol) + if (std::abs(MAT(i, j)) > abs_tol) nnz += 1; } } diff --git a/RandBLAS/sparse_data/coo_matrix.hh b/RandBLAS/sparse_data/coo_matrix.hh index 33b15f88..bd6739cf 100644 --- a/RandBLAS/sparse_data/coo_matrix.hh +++ b/RandBLAS/sparse_data/coo_matrix.hh @@ -36,6 +36,7 @@ #include "RandBLAS/util.hh" +#include #include #include #include @@ -360,7 +361,7 @@ void dense_to_coo(int64_t stride_row, int64_t stride_col, T *mat, T abs_tol, COO for (int64_t i = 0; i < n_rows; ++i) { for (int64_t j = 0; j < n_cols; ++j) { T val = MAT(i, j); - if (abs(val) > abs_tol) { + if (std::abs(val) > abs_tol) { spmat.vals[nnz] = val; spmat.rows[nnz] = i; spmat.cols[nnz] = j; @@ -414,4 +415,3 @@ void coo_to_dense(const COOMatrix &spmat, Layout layout, T *mat) { } } // end namespace RandBLAS::sparse_data::coo - diff --git a/RandBLAS/sparse_data/csc_matrix.hh b/RandBLAS/sparse_data/csc_matrix.hh index 7789b52c..d206735d 100644 --- a/RandBLAS/sparse_data/csc_matrix.hh +++ b/RandBLAS/sparse_data/csc_matrix.hh @@ -34,6 +34,7 @@ #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/conversions.hh" #include +#include namespace RandBLAS::sparse_data { @@ -307,7 +308,7 @@ void dense_to_csc(int64_t stride_row, int64_t stride_col, T *mat, T abs_tol, CSC for (int64_t j = 0; j < n_cols; ++j) { for (int64_t i = 0; i < n_rows; ++i) { T val = MAT(i, j); - if (abs(val) > abs_tol) { + if (std::abs(val) > abs_tol) { spmat.vals[nnz] = val; spmat.rowidxs[nnz] = i; nnz += 1; diff --git a/RandBLAS/sparse_data/csr_matrix.hh b/RandBLAS/sparse_data/csr_matrix.hh index 026e1b66..d868a72c 100644 --- a/RandBLAS/sparse_data/csr_matrix.hh +++ b/RandBLAS/sparse_data/csr_matrix.hh @@ -34,6 +34,7 @@ #include "RandBLAS/sparse_data/base.hh" #include "RandBLAS/sparse_data/conversions.hh" #include +#include namespace RandBLAS::sparse_data { @@ -311,7 +312,7 @@ void dense_to_csr(int64_t stride_row, int64_t stride_col, T *mat, T abs_tol, CSR for (int64_t i = 0; i < n_rows; ++i) { for (int64_t j = 0; j < n_cols; ++j) { T val = MAT(i, j); - if (abs(val) > abs_tol) { + if (std::abs(val) > abs_tol) { spmat.vals[nnz] = val; spmat.colidxs[nnz] = j; nnz += 1; diff --git a/RandBLAS/sparse_data/sksp.hh b/RandBLAS/sparse_data/sksp.hh index f3eefdd8..d0f0eab9 100644 --- a/RandBLAS/sparse_data/sksp.hh +++ b/RandBLAS/sparse_data/sksp.hh @@ -40,7 +40,7 @@ namespace RandBLAS::sparse_data { // ============================================================================= /// \fn lsksp3(blas::Layout layout, blas::Op opS, blas::Op opA, int64_t d, -/// int64_t n, int64_t m, T alpha, const DenseSkOp &S, int64_t ro_s, int64_t co_s, +/// int64_t n, int64_t m, T alpha, const DenseSkOp &S, int64_t ro_s, int64_t co_s, /// SpMat &A, int64_t ro_a, int64_t co_a, T beta, T *B, int64_t ldb /// ) /// @verbatim embed:rst:leading-slashes @@ -186,7 +186,7 @@ void lsksp3( // ============================================================================= /// \fn rsksp3(blas::Layout layout, blas::Op opA, blas::Op opS, int64_t m, /// int64_t d, int64_t n, T alpha, const SpMat &A, int64_t ro_a, int64_t co_a, -/// const DenseSkOp &S, int64_t ro_s, int64_t co_s, T beta, T *B, int64_t ldb +/// const DenseSkOp &S, int64_t ro_s, int64_t co_s, T beta, T *B, int64_t ldb /// ) /// @verbatim embed:rst:leading-slashes /// Sketch from the right in an SpMM-like operation diff --git a/RandBLAS/sparse_skops.hh b/RandBLAS/sparse_skops.hh index 70e4c55f..ccb8900d 100644 --- a/RandBLAS/sparse_skops.hh +++ b/RandBLAS/sparse_skops.hh @@ -40,15 +40,19 @@ #include #include #include +#include +#include #include #include +#include #define MAX(a, b) (((a) < (b)) ? (b) : (a)) #define MIN(a, b) (((a) < (b)) ? (a) : (b)) namespace RandBLAS::sparse { -template > +template void _considerate_fisher_yates( const state_t &state, int64_t k, @@ -63,25 +67,37 @@ void _considerate_fisher_yates( // indices = {0, 1, 2, ..., n - 1}; input-output; not const. // work_piv = buffer of length k; output-only. randblas_require( k <= n ); - if (vals != nullptr) { - randblas_require(state.len_c >= 4); - } - typename state_t::generator gen; - auto ctr = state.counter; + using res_t = typename state_t::res_t; + using word_t = typename res_t::value_type; + constexpr std::size_t block_size = std::tuple_size_v; + constexpr auto word_bits = std::numeric_limits::digits; + static_assert(word_bits == 32 || word_bits == 64, + "sparse Fisher-Yates sampling requires 32- or 64-bit result words"); + static_assert((word_bits == 32 && block_size >= 3) || + (word_bits == 64 && block_size >= 2), + "sparse Fisher-Yates sampling requires enough result words for an index and sign"); + state_t work = state; for (sint_t j = 0; j < k; ++j) { - auto rv = gen(ctr, state.key); - // ^ Array of uint32's, sampled uniformly at random. - ctr.incr(); + res_t rv{}; + work.generate(rv); + // ^ Array of unsigned words, sampled uniformly at random. + work.advance(1); // ^ The counter is incremented every loop, even if this means // we aren't being super efficient in terms of samples. - auto s = promote_uint_pair(rv[0], rv[1]); + std::uint64_t s; + if constexpr (word_bits == 32) { + s = promote_uint_pair(rv[0], rv[1]); + } else { + s = rv[0]; + } sint_t p = j + static_cast(s % (n - j)); // ^ sample from {j, j+1, ...., n - 1} work_piv[j] = p; std::swap(indices[p], indices[j]); samples[j] = indices[j]; if (vals != nullptr) { - vals[j] = (rv[2] % 2 == 0) ? 1.0 : -1.0; + constexpr std::size_t sign_lane = word_bits == 32 ? 2 : 1; + vals[j] = (rv[sign_lane] % 2 == 0) ? 1.0 : -1.0; } } for (sint_t j = 1; j <= k; ++j) { @@ -94,7 +110,8 @@ void _considerate_fisher_yates( return; } -template > +template static state_t repeated_fisher_yates( const state_t &state, int64_t vec_nnz, @@ -121,14 +138,13 @@ static state_t repeated_fisher_yates( std::vector vec_work(dim_major); std::iota(vec_work.begin(), vec_work.end(), 0); std::vector pivots(vec_nnz); - auto [ctr, key] = state; + state_t work = state; for (sint_t i = 0; i < dim_minor; ++i) { - state_t state_work{ctr, state.key}; _considerate_fisher_yates( - state_work, vec_nnz, dim_major, + work, vec_nnz, dim_major, idxs_major, vec_work.data(), pivots.data(), vals ); - ctr.incr(vec_nnz); + work.advance(vec_nnz); idxs_major += vec_nnz; if (idxs_minor != nullptr) { std::fill(idxs_minor, idxs_minor + vec_nnz, i); @@ -138,7 +154,7 @@ static state_t repeated_fisher_yates( vals += vec_nnz; } } - return state_t {ctr, key}; + return work; } inline double isometry_scale(Axis major_axis, int64_t vec_nnz, int64_t dim_major, int64_t dim_minor) { @@ -155,7 +171,8 @@ namespace RandBLAS { // Forward declaration of SparseSkOp. It's returnable by // SparseDist.sample(), but its definition involves SparseDist. -template +template struct SparseSkOp; // ============================================================================= @@ -261,8 +278,9 @@ struct SparseDist { // ------------------------------------------------------------------------------------- /// Construct a SparseSkOp with this distribution and the provided seed_state. - template - SparseSkOp sample(RNGState &seed_state) { + template + SparseSkOp sample(state_t &seed_state) { return {*this, seed_state}; } @@ -292,22 +310,23 @@ struct SparseDist { /// be used for the next call to a random sampling function whose output should be statistically /// independent from \math{\ttt{samples}.} /// -template > +template inline state_t repeated_fisher_yates( int64_t k, int64_t n, int64_t r, sint_t *samples, const state_t &state ) { return sparse::repeated_fisher_yates(state, k, n, r, samples, (sint_t*) nullptr, (double*) nullptr); } -template -RNGState compute_next_state(SparseDist dist, RNGState state) { +template +state_t compute_next_state(SparseDist dist, state_t state) { // Both _considerate_fisher_yates (SASO with vec_nnz > 1) and // sample_indices_iid_uniform (SASO with vec_nnz == 1, and LASO) consume // exactly one CBRNG counter increment per nonzero. int64_t num_major_axis_vec = (dist.major_axis == Axis::Short) ? std::max(dist.n_rows, dist.n_cols) : std::min(dist.n_rows, dist.n_cols); - state.counter.incr(num_major_axis_vec * dist.vec_nnz); + state.advance(num_major_axis_vec * dist.vec_nnz); return state; } @@ -315,7 +334,7 @@ RNGState compute_next_state(SparseDist dist, RNGState state) { /// A sample from a distribution over structured sparse matrices with either /// independent rows or independent columns. This type conforms to the /// SketchingOperator concept. -template +template struct SparseSkOp { // --------------------------------------------------------------------------- @@ -324,7 +343,7 @@ struct SparseSkOp { // --------------------------------------------------------------------------- /// Type alias. - using state_t = RNGState; + using state_t = generator_state_t; // --------------------------------------------------------------------------- /// Real scalar type used for nonzeros in matrix representations of this operator. @@ -457,7 +476,7 @@ struct SparseSkOp { nnz(nnz), vals(vals), rows(rows), cols(cols){ }; // Move constructor - SparseSkOp(SparseSkOp &&S + SparseSkOp(SparseSkOp &&S ) : dist(S.dist), seed_state(S.seed_state), next_state(S.next_state), n_rows(dist.n_rows), n_cols(dist.n_cols), own_memory(S.own_memory), nnz(S.nnz), rows(S.rows), cols(S.cols), vals(S.vals) @@ -565,7 +584,7 @@ void laso_merge_long_axis_vector_coo_data( /// - A CBRNG state used to define :math:`\mtxS.` /// /// @endverbatim -template +template state_t fill_sparse_unpacked( const SparseDist &D, int64_t n_rows_sub, int64_t n_cols_sub, @@ -625,7 +644,7 @@ state_t fill_sparse_unpacked( // and LASO) consume exactly vec_nnz counter increments per major-axis vector, so the // skip amount is uniform. state_t work_state = seed_state; - work_state.counter.incr(num_major_off * vec_nnz); + work_state.advance(num_major_off * vec_nnz); // Identify which output array holds the major-axis coordinate and which holds the // minor-axis coordinate (the index of the major-axis vector). We sample directly @@ -709,7 +728,7 @@ state_t fill_sparse_unpacked( // ro_s = co_s = 0 and the full operator dimensions instead. It writes the COO data for // the operator (D, seed_state) into the first nnz entries of (vals, rows, cols), which // must have length at least D.full_nnz. -template +template state_t fill_sparse_unpacked_nosub( const SparseDist &D, int64_t &nnz, T* vals, sint_t* rows, sint_t *cols, diff --git a/RandBLAS/testing/comparison.hh b/RandBLAS/testing/comparison.hh index b52fb35f..44e00679 100644 --- a/RandBLAS/testing/comparison.hh +++ b/RandBLAS/testing/comparison.hh @@ -70,11 +70,11 @@ bool approx_equal(T A, T B, std::ostream &str, { // Check if the numbers are really close -- needed // when comparing numbers near zero. - T diff_ab = abs(A - B); + T diff_ab = std::abs(A - B); if (diff_ab <= atol) return true; - T max_ab = std::max(abs(B), abs(A)); + T max_ab = std::max(std::abs(B), std::abs(A)); if (diff_ab <= max_ab * rtol) return true; diff --git a/RandBLAS/testing/lapack_like.hh b/RandBLAS/testing/lapack_like.hh index dd07ff98..e625eedc 100644 --- a/RandBLAS/testing/lapack_like.hh +++ b/RandBLAS/testing/lapack_like.hh @@ -220,8 +220,9 @@ inline int64_t required_powermethod_iters(int64_t n, T p_fail, T tol) { return num_iters; } -template -std::pair> power_method(int64_t n, FUNC &A, T* v, T tol, T failure_prob, const RNGState &state) { +template +std::pair power_method(int64_t n, FUNC &A, T* v, T tol, + T failure_prob, const state_t &state) { auto next_state = RandBLAS::fill_dense_unpacked(blas::Layout::ColMajor, {n, 1}, n, 1, 0, 0, v, state); std::vector work(n, 0.0); T* u = work.data(); @@ -241,8 +242,11 @@ std::pair> power_method(int64_t n, FUNC &A, T* v, T tol, T fail } -template -std::tuple> exeigs_powermethod(int64_t n, const T* A, T* eigvecs, T tol, T failure_prob, const RNGState &state, std::vector work) { +template +std::tuple exeigs_powermethod(int64_t n, const T* A, + T* eigvecs, T tol, T failure_prob, + const state_t &state, + std::vector work) { auto layout = blas::Layout::ColMajor; RandBLAS::util::require_symmetric(layout, A, n, n, (T) 0.0); diff --git a/RandBLAS/testing/linops.hh b/RandBLAS/testing/linops.hh index 6082ca9d..86339f96 100644 --- a/RandBLAS/testing/linops.hh +++ b/RandBLAS/testing/linops.hh @@ -36,6 +36,7 @@ #include "RandBLAS/skge.hh" #include "RandBLAS/sparse_data/spmm_dispatch.hh" #include "RandBLAS/util.hh" +#include #include #include #include @@ -74,12 +75,12 @@ std::vector eye(int64_t n) { return A; } -template -auto random_matrix(int64_t m, int64_t n, RNGState s) { +template +auto random_matrix(int64_t m, int64_t n, state_t s) { std::vector A(m * n); DenseDist DA(m, n); auto next_state = RandBLAS::fill_dense(DA, A.data(), s); - std::tuple, Layout, RNGState> t{A, DA.natural_layout, next_state}; + std::tuple, Layout, state_t> t{A, DA.natural_layout, next_state}; return t; } @@ -237,7 +238,7 @@ void reference_left_apply( for (int64_t i = 0; i < rows_S; ++i) { for (int64_t j = 0; j < cols_S; ++j) { auto ell = i * s_row_stride + j * s_col_stride; - S_dense_abs[ell] = abs(S_dense[ell]); + S_dense_abs[ell] = std::abs(S_dense[ell]); } } @@ -251,14 +252,14 @@ void reference_left_apply( std::vector A_abs_vec(size_A); T* A_abs = A_abs_vec.data(); for (int64_t i = 0; i < size_A; ++i) - A_abs[i] = abs(A[i]); + A_abs[i] = std::abs(A[i]); if (beta != 0.0) { for (int64_t i = 0; i < size_B; ++i) - E[i] = abs(B[i]); + E[i] = std::abs(B[i]); } T eps = std::numeric_limits::epsilon(); - T err_alpha = (abs(alpha) * m) * (2 * eps); - T err_beta = abs(beta) * eps; + T err_alpha = (std::abs(alpha) * m) * (2 * eps); + T err_beta = std::abs(beta) * eps; T* S_abs_ptr = S_dense_abs.data(); blas::gemm(layout, transS, transA, d, n, m, err_alpha, &S_abs_ptr[pos], lds, A_abs, lda, err_beta, E, ldb diff --git a/RandBLAS/testing/rng.hh b/RandBLAS/testing/rng.hh new file mode 100644 index 00000000..504e260e --- /dev/null +++ b/RandBLAS/testing/rng.hh @@ -0,0 +1,95 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#pragma once + +#include "RandBLAS/random_gen.hh" + +#include +#include +#include +#include + +namespace RandBLAS::testing::detail { + +/// Sequential word stream used by RandBLAS test-data generators. +template +struct RNGStream { + using res_t = typename state_t::res_t; + using word_t = typename res_t::value_type; + static constexpr std::size_t block_size = std::tuple_size_v; + + state_t state; + res_t buffer{}; + std::size_t pos = block_size; + double spare = 0.0; + bool has_spare = false; + + explicit RNGStream(state_t const& initial_state) : state(initial_state) {} + + word_t next_word() { + if (pos >= block_size) { + state.generate(buffer); + state.advance(1); + pos = 0; + } + return buffer[pos++]; + } + + /// Return a uniform value in (0, 1]. + double uniform_01() { + return rng::u01(next_word()); + } + + /// Return one normal value, caching the second Box--Muller result. + template + value_t gaussian() { + if (has_spare) { + has_spare = false; + return static_cast(spare); + } + word_t angle_word = next_word(); + word_t radius_word = next_word(); + auto [first, second] = rng::boxmuller(angle_word, radius_word); + spare = second; + has_spare = true; + return static_cast(first); + } + + /// Return the number of failures before the first Bernoulli success. + std::int64_t geometric(double log_1_minus_p) { + double u = uniform_01(); + return static_cast( + std::floor(std::log(1.0 - u) / log_1_minus_p)); + } + + /// Report state after every result block already loaded into the buffer. + state_t get_state() const { return state; } +}; + +} // namespace RandBLAS::testing::detail diff --git a/RandBLAS/testing/sparse_data.hh b/RandBLAS/testing/sparse_data.hh index 43671953..a5515c3a 100644 --- a/RandBLAS/testing/sparse_data.hh +++ b/RandBLAS/testing/sparse_data.hh @@ -47,6 +47,7 @@ #include "RandBLAS/sparse_data/csr_matrix.hh" #include "RandBLAS/sparse_data/csc_matrix.hh" #include "RandBLAS/sparse_data/conversions.hh" +#include "RandBLAS/testing/rng.hh" namespace RandBLAS::testing { @@ -64,77 +65,8 @@ using RandBLAS::SignedInteger; #endif -namespace detail { - -// Sequential wrapper around a Random123 CBRNG. Philox4x32 produces 4 uint32_t -// values per counter increment; this helper dispenses them one at a time and -// provides uniform, Gaussian, and geometric draws. -// -// Subclasses RNGState so that counter and key are inherited directly. -template -struct PhiloxStream : public RandBLAS::RNGState { - using state_t = RandBLAS::RNGState; - using ctr_t = typename state_t::ctr_type; - static constexpr int ctr_size = state_t::len_c; - - RNG rng; - ctr_t buffer; - int pos; - double spare; - bool has_spare; - - PhiloxStream(const state_t &state) - : state_t(state), pos(ctr_size), spare(0.0), has_spare(false) {} - - uint32_t next_u32() { - if (pos >= ctr_size) { - buffer = rng(this->counter, this->key); - this->counter.incr(); - pos = 0; - } - return buffer.v[pos++]; - } - - // Uniform in (0, 1], never 0.0 (safe for log). - double uniform_01() { - return r123::u01(next_u32()); - } - - // Box-Muller Gaussian. Each call to boxmuller produces two independent values; - // we cache the second and return it on the next invocation. - template - T gaussian() { - if (has_spare) { - has_spare = false; - return static_cast(spare); - } - uint32_t u1 = next_u32(); - uint32_t u2 = next_u32(); - auto [g1, g2] = r123::boxmuller(u1, u2); - spare = g2; - has_spare = true; - return static_cast(g1); - } - - // Geometric distribution: number of failures before first success in Bernoulli(p). - // Uses inverse CDF: floor(log(1 - u) / log(1 - p)) with u ~ Uniform(0, 1]. - // Since u01 returns (0, 1], we have 1 - u in [0, 1). The only problematic value - // is 1 - u = 0, i.e., u = 1.0 exactly. With u01(uint32_t), the maximum - // is 1.0 - 2^-33, so this never happens. - int64_t geometric(double log_1_minus_p) { - double u = uniform_01(); - return static_cast(std::floor(std::log(1.0 - u) / log_1_minus_p)); - } - - state_t get_state() const { - return state_t{this->counter, this->key}; - } -}; - -} // end namespace detail - - -template +template void iid_sparsify_random_dense( int64_t n_rows, int64_t n_cols, @@ -142,7 +74,7 @@ void iid_sparsify_random_dense( int64_t stride_col, T* mat, T prob_of_zero, - RandBLAS::RNGState state + state_t state ) { auto spar = new T[n_rows * n_cols]; auto dist = RandBLAS::DenseDist(n_rows, n_cols, RandBLAS::ScalarDist::Uniform); @@ -173,14 +105,15 @@ void iid_sparsify_random_dense( } -template +template void iid_sparsify_random_dense( int64_t n_rows, int64_t n_cols, Layout layout, T* mat, T prob_of_zero, - RandBLAS::RNGState state + state_t state ) { if (layout == Layout::ColMajor) { iid_sparsify_random_dense(n_rows, n_cols, 1, n_rows, mat, prob_of_zero, state); @@ -261,17 +194,18 @@ int64_t trianglize_coo( // Returns {CSRMatrix, next_state}. Use with structured bindings: // auto [A, next_state] = random_csr(m, n, density, state); // ============================================================================ -template -std::pair, RandBLAS::RNGState> random_csr( +template +std::pair, state_t> random_csr( int64_t m, int64_t n, double density, - const RandBLAS::RNGState &state + const state_t &state ) { randblas_require(density >= 0.0 && density <= 1.0); CSRMatrix A(m, n); - detail::PhiloxStream stream(state); + detail::RNGStream stream(state); if (density == 0.0 || m == 0 || n == 0) { if (m > 0) { @@ -343,17 +277,18 @@ std::pair, RandBLAS::RNGState> random_csr( // Returns {CSCMatrix, next_state}. Use with structured bindings: // auto [A, next_state] = random_csc(m, n, density, state); // ============================================================================ -template -std::pair, RandBLAS::RNGState> random_csc( +template +std::pair, state_t> random_csc( int64_t m, int64_t n, double density, - const RandBLAS::RNGState &state + const state_t &state ) { randblas_require(density >= 0.0 && density <= 1.0); CSCMatrix A(m, n); - detail::PhiloxStream stream(state); + detail::RNGStream stream(state); if (density == 0.0 || m == 0 || n == 0) { if (n > 0) { @@ -422,17 +357,18 @@ std::pair, RandBLAS::RNGState> random_csc( // Returns {COOMatrix, next_state}. Use with structured bindings: // auto [A, next_state] = random_coo(m, n, density, state); // ============================================================================ -template -std::pair, RandBLAS::RNGState> random_coo( +template +std::pair, state_t> random_coo( int64_t m, int64_t n, double density, - const RandBLAS::RNGState &state + const state_t &state ) { randblas_require(density >= 0.0 && density <= 1.0); COOMatrix A(m, n); - detail::PhiloxStream stream(state); + detail::RNGStream stream(state); if (density == 0.0 || m == 0 || n == 0) { return {std::move(A), stream.get_state()}; diff --git a/RandBLAS/util.hh b/RandBLAS/util.hh index dcff6b5a..970b0894 100644 --- a/RandBLAS/util.hh +++ b/RandBLAS/util.hh @@ -33,8 +33,6 @@ #include #include #include -#include -#include #include #include @@ -44,10 +42,13 @@ # include #endif #include +#include #include #include #include #include +#include +#include namespace RandBLAS::util { @@ -134,8 +135,8 @@ void require_symmetric(blas::Layout layout, const T* A, int64_t n, int64_t lda, for (int64_t j = i+1; j < n; ++j) { T Aij = matA(i,j); T Aji = matA(j,i); - T viol = abs(Aij - Aji); - T rel_tol = (abs(Aij) + abs(Aji) + 1)*tol; + T viol = std::abs(Aij - Aji); + T rel_tol = (std::abs(Aij) + std::abs(Aji) + 1)*tol; if (viol > rel_tol) { std::string message = "Symmetry check failed. |A(%i,%i) - A(%i,%i)| was %e, which exceeds tolerance of %e."; auto _message = message.c_str(); @@ -490,13 +491,12 @@ static inline TO uneg11_to_u01(TI in) { /// be used for the next call to a random sampling function whose output should be statistically /// independent from :math:`\ttt{samples}.` /// @endverbatim -template > +template state_t sample_indices_iid(int64_t n, const T* cdf, int64_t k, sint_t* samples, const state_t &state) { - auto [ctr, key] = state; - using RNG = typename state_t::generator; - RNG gen; - auto rv_array = r123ext::uneg11::generate(gen, ctr, key); - int64_t len_c = (int64_t) state.len_c; + state_t work = state; + auto rv_array = rng::uneg11::generate(work); + constexpr int64_t len_c = std::tuple_size_v; int64_t rv_index = 0; for (int64_t i = 0; i < k; ++i) { auto random_unif01 = uneg11_to_u01(rv_array[rv_index]); @@ -504,41 +504,63 @@ state_t sample_indices_iid(int64_t n, const T* cdf, int64_t k, sint_t* samples, samples[i] = sample_index; rv_index += 1; if (rv_index == len_c) { - ctr.incr(1); - rv_array = r123ext::uneg11::generate(gen, ctr, key); + work.advance(1); + if (i + 1 < k) { + rv_array = rng::uneg11::generate(work); + } rv_index = 0; } } - if (0 < rv_index) ctr.incr(1); - return state_t(ctr, key); + if (0 < rv_index) { + work.advance(1); + } + return work; } inline std::uint64_t promote_uint_pair(std::uint32_t a, std::uint32_t b) { return static_cast(a) + (static_cast(b) << 32); } -template > +template state_t sample_indices_iid_uniform(int64_t n, int64_t k, sint_t* samples, T* rademachers, const state_t &state) { - if constexpr (WriteRademachers) { - randblas_require(state.len_c >= 4); - } else { - randblas_require(state.len_c >= 2); + using res_t = typename state_t::res_t; + using word_t = typename res_t::value_type; + constexpr std::size_t block_size = std::tuple_size_v; + constexpr auto word_bits = std::numeric_limits::digits; + static_assert(word_bits == 32 || word_bits == 64, + "uniform index sampling requires 32- or 64-bit result words"); + if constexpr (word_bits == 32) { + static_assert(block_size >= 2, + "32-bit index sampling requires at least two result words"); + if constexpr (WriteRademachers) { + static_assert(block_size >= 3, + "32-bit Rademacher sampling requires a third result word"); + } + } else if constexpr (WriteRademachers) { + static_assert(block_size >= 2, + "64-bit Rademacher sampling requires two result words"); } - using RNG = typename state_t::generator; - RNG gen; - auto ctr = state.counter; - auto key = state.key; + + state_t work = state; std::uint64_t n_64 = static_cast(n); for (int64_t i = 0; i < k; ++i) { - auto rv = gen(ctr, key); - ctr.incr(); - std::uint64_t s = promote_uint_pair(rv[0], rv[1]); + res_t rv{}; + work.generate(rv); + work.advance(1); + std::uint64_t s; + if constexpr (word_bits == 32) { + s = promote_uint_pair(rv[0], rv[1]); + } else { + s = rv[0]; + } samples[i] = static_cast(s % n_64); if constexpr (WriteRademachers) { - rademachers[i] = (rv[2] % 2 == 0) ? (T) 1 : (T) -1; + constexpr std::size_t sign_lane = word_bits == 32 ? 2 : 1; + rademachers[i] = (rv[sign_lane] % 2 == 0) ? (T) 1 : (T) -1; } } - return state_t(ctr, key); + return work; } @@ -550,11 +572,11 @@ state_t sample_indices_iid_uniform(int64_t n, int64_t k, sint_t* samples, T* rad /// independent from :math:`\ttt{samples}.` /// /// @endverbatim -template > +template state_t sample_indices_iid_uniform(int64_t n, int64_t k, sint_t* samples, const state_t &state) { return sample_indices_iid_uniform(n, k, samples, (float*) nullptr, state); } } // end namespace RandBLAS - diff --git a/docker/tsan/Dockerfile b/docker/tsan/Dockerfile index 0acccac7..cb25a761 100644 --- a/docker/tsan/Dockerfile +++ b/docker/tsan/Dockerfile @@ -6,9 +6,9 @@ # Build: docker build -t randblas-tsan docker/tsan/ # Run: ./docker/tsan/run.sh # -# The image bakes blaspp + Random123 into /opt so that iterating on RandBLAS -# itself is fast — only RandBLAS gets rebuilt per session. Bump the upstream -# refs in this Dockerfile to refresh those dependencies. +# The image bakes blaspp into /opt so that iterating on RandBLAS itself is +# fast — only RandBLAS gets rebuilt per session. Bump the upstream ref in this +# Dockerfile to refresh that dependency. FROM ubuntu:24.04 ENV DEBIAN_FRONTEND=noninteractive @@ -38,15 +38,7 @@ RUN git clone --depth 1 https://github.com/icl-utk-edu/blaspp.git /tmp/blaspp \ && cmake --build /tmp/blaspp-build -j"$(nproc)" --target install \ && rm -rf /tmp/blaspp /tmp/blaspp-build -# Random123: header-only. cp -R is portable (BSD/GNU) so we don't need the -# Random123 Makefile's `install-include` target (which uses GNU `cp -d`). -RUN git clone --depth 1 https://github.com/DEShawResearch/Random123.git /tmp/Random123 \ - && mkdir -p /opt/Random123-install/include \ - && cp -R /tmp/Random123/include/Random123 /opt/Random123-install/include/ \ - && rm -rf /tmp/Random123 - ENV blaspp_DIR=/opt/blaspp-install/lib/cmake/blaspp -ENV Random123_DIR=/opt/Random123-install/include WORKDIR /work CMD ["/bin/bash"] diff --git a/docker/tsan/run.sh b/docker/tsan/run.sh index c38e9244..0f489c0e 100755 --- a/docker/tsan/run.sh +++ b/docker/tsan/run.sh @@ -53,7 +53,6 @@ case "${subcommand}" in -DCMAKE_CXX_FLAGS="-fsanitize=thread -fno-omit-frame-pointer" \ -DCMAKE_EXE_LINKER_FLAGS="-fsanitize=thread" \ -Dblaspp_DIR="${blaspp_DIR}" \ - -DRandom123_DIR="${Random123_DIR}" \ /work/RandBLAS make -j"$(nproc)" # ignore_noninstrumented_modules=1 silences spurious reports diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 2595d374..368f3315 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -32,9 +32,6 @@ set( add_executable( tls_dense_skop ${tls_dense_skop_cxx} ) -target_include_directories( - tls_dense_skop PUBLIC ${Random123_DIR} -) target_link_libraries( tls_dense_skop PUBLIC RandBLAS blaspp lapackpp ) @@ -45,9 +42,6 @@ set( add_executable( tls_sparse_skop ${tls_sparse_skop_cxx} ) -target_include_directories( - tls_sparse_skop PUBLIC ${Random123_DIR} -) target_link_libraries( tls_sparse_skop PUBLIC RandBLAS blaspp lapackpp ) @@ -55,9 +49,6 @@ target_link_libraries( add_executable( slra_svd_synthetic sparse-low-rank-approx/svd_rank1_plus_noise.cc ) -target_include_directories( - slra_svd_synthetic PUBLIC ${Random123_DIR} -) target_link_libraries( slra_svd_synthetic PUBLIC RandBLAS blaspp lapackpp ) @@ -77,9 +68,6 @@ FetchContent_MakeAvailable( add_executable( slra_svd_fmm sparse-low-rank-approx/svd_matrixmarket.cc ) -target_include_directories( - slra_svd_fmm PUBLIC ${Random123_DIR} -) target_link_libraries( slra_svd_fmm PUBLIC RandBLAS blaspp lapackpp fast_matrix_market::fast_matrix_market ) @@ -87,9 +75,6 @@ target_link_libraries( add_executable( slra_qrcp sparse-low-rank-approx/qrcp_matrixmarket.cc ) -target_include_directories( - slra_qrcp PUBLIC ${Random123_DIR} -) target_link_libraries( slra_qrcp PUBLIC RandBLAS blaspp lapackpp fast_matrix_market::fast_matrix_market ) @@ -97,9 +82,6 @@ target_link_libraries( add_executable( spmm_performance simple-kernel-benchmarks/spmm_performance.cc ) -target_include_directories( - spmm_performance PUBLIC ${Random123_DIR} -) target_link_libraries( spmm_performance PUBLIC RandBLAS blaspp lapackpp ) @@ -107,9 +89,6 @@ target_link_libraries( add_executable( sketch_general_performance simple-kernel-benchmarks/sketch_general_performance.cc ) -target_include_directories( - sketch_general_performance PUBLIC ${Random123_DIR} -) target_link_libraries( sketch_general_performance PUBLIC RandBLAS blaspp lapackpp ) diff --git a/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc b/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc index 70e8dc0f..c158087f 100644 --- a/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc +++ b/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc @@ -387,7 +387,7 @@ int run(SpMat &A, int64_t k, int64_t power_iteration_steps, StabilizationMethod T *Q = new T[m*k]{}; T *R = new T[k*n]{}; int64_t *piv = new int64_t[n]{}; - RandBLAS::RNGState state(0); + RandBLAS::DefaultRNGState state(0); auto start_timer = std_clock::now(); TIMED_LINE( diff --git a/examples/sparse-low-rank-approx/svd_matrixmarket.cc b/examples/sparse-low-rank-approx/svd_matrixmarket.cc index 243cde03..1673ca68 100644 --- a/examples/sparse-low-rank-approx/svd_matrixmarket.cc +++ b/examples/sparse-low-rank-approx/svd_matrixmarket.cc @@ -220,7 +220,7 @@ int main(int argc, char** argv) { double *U = new double[m*k]{}; double *VT = new double[k*n]{}; double *qb_work = new double[std::max(m, n)]; - RandBLAS::RNGState state(0); + RandBLAS::DefaultRNGState state(0); /* Effect of various parameters on performance: It's EXTREMELY important to use -O3 if you want reasonably diff --git a/examples/sparse-low-rank-approx/svd_rank1_plus_noise.cc b/examples/sparse-low-rank-approx/svd_rank1_plus_noise.cc index 23c28676..0250cdb6 100644 --- a/examples/sparse-low-rank-approx/svd_rank1_plus_noise.cc +++ b/examples/sparse-low-rank-approx/svd_rank1_plus_noise.cc @@ -68,9 +68,11 @@ auto parse_dimension_args(int argc, char** argv) { return std::make_tuple(m, n, vec_nnz); } -template +template void iid_sparsify_random_dense( - int64_t n_rows, int64_t n_cols, int64_t stride_row, int64_t stride_col, T* mat, T prob_of_zero, RandBLAS::RNGState state + int64_t n_rows, int64_t n_cols, int64_t stride_row, + int64_t stride_col, T* mat, T prob_of_zero, State state ) { auto spar = new T[n_rows * n_cols]; auto dist = RandBLAS::DenseDist(n_rows, n_cols, RandBLAS::ScalarDist::Uniform); diff --git a/examples/total-least-squares/tls_dense_skop.cc b/examples/total-least-squares/tls_dense_skop.cc index 543add01..2e28a1a5 100644 --- a/examples/total-least-squares/tls_dense_skop.cc +++ b/examples/total-least-squares/tls_dense_skop.cc @@ -38,6 +38,7 @@ #include #include #include +#include #include using std::chrono::high_resolution_clock; @@ -138,8 +139,8 @@ int main(int argc, char* argv[]){ // Sample the sketching operator auto time_constructsketch1 = high_resolution_clock::now(); RandBLAS::DenseDist Dist{ sk_dim, m }; - uint32_t seed = 1997; - RandBLAS::DenseSkOp S(Dist, seed); + std::uint64_t seed = 1997; + RandBLAS::DenseSkOp S(Dist, seed); RandBLAS::fill_dense(S); auto time_constructsketch2 = high_resolution_clock::now(); double sampling_time = (double) duration_cast(time_constructsketch2 - time_constructsketch1).count()/1000; diff --git a/examples/total-least-squares/tls_sparse_skop.cc b/examples/total-least-squares/tls_sparse_skop.cc index 054fff86..6a3ba39d 100644 --- a/examples/total-least-squares/tls_sparse_skop.cc +++ b/examples/total-least-squares/tls_sparse_skop.cc @@ -38,6 +38,7 @@ #include #include #include +#include #include using std::chrono::high_resolution_clock; @@ -145,8 +146,8 @@ int main(int argc, char* argv[]){ 8, // Number of non-zero entires per column, RandBLAS::Axis::Short // A "SASO" (aka SJLT, aka OSNAP, aka generalized CountSketch) ); - uint32_t seed = 1997; - RandBLAS::SparseSkOp S(Dist, seed); + std::uint64_t seed = 1997; + RandBLAS::SparseSkOp S(Dist, seed); RandBLAS::fill_sparse(S); auto time_constructsketch2 = high_resolution_clock::now(); double sampling_time = (double) duration_cast(time_constructsketch2 - time_constructsketch1).count()/1000; diff --git a/rtd/source/FAQ.rst b/rtd/source/FAQ.rst index ea9a0e1c..c0ecc2d4 100644 --- a/rtd/source/FAQ.rst +++ b/rtd/source/FAQ.rst @@ -124,8 +124,9 @@ C++ idioms and features we do use Things that affect our API: * Templates. We template for floating point precision just about everywhere. - We also template for stateful random number generators (see :cpp:any:`RandBLAS::RNGState`) - and arrays of 32-bit versus 64-bit signed integers. + Sampling functions and sketching operators also template on random-number + state types satisfying :cpp:any:`RandBLAS::GeneratorState`, and on arrays + of 32-bit versus 64-bit signed integers. * Standard constructors. We use these for any nontrivial struct type in RandBLAS. They're important because many of our datatypes have const members that need to be initialized as functions (albeit simple functions) of other members. @@ -188,4 +189,3 @@ Some discussion We have no plans for consistent naming of overload-free sparse BLAS functions. The most we do in this regard is offer functions called [left/right]_spmm for SpMM where the sparse matrix operand appears on the left or on the right. - diff --git a/rtd/source/api_reference/skops_and_dists.rst b/rtd/source/api_reference/skops_and_dists.rst index 8a7c35e5..a54e8141 100644 --- a/rtd/source/api_reference/skops_and_dists.rst +++ b/rtd/source/api_reference/skops_and_dists.rst @@ -41,6 +41,13 @@ Preliminaries .. doxygenenum:: RandBLAS::Axis :project: RandBLAS +.. dropdown:: GeneratorState + :animate: fade-in-slide-down + :color: light + + .. doxygenconcept:: RandBLAS::GeneratorState + :project: RandBLAS + .. dropdown:: RNGState :animate: fade-in-slide-down :color: light @@ -77,7 +84,7 @@ Dense sketching, with Gaussians *et al.* .. doxygenfunction:: RandBLAS::fill_dense(DenseSkOp &S) :project: RandBLAS - .. doxygenfunction:: RandBLAS::fill_dense_unpacked(blas::Layout layout, const DenseDist &D, int64_t n_rows, int64_t n_cols, int64_t S_ro, int64_t S_co, T *buff, const RNGState &seed) + .. doxygenfunction:: RandBLAS::fill_dense_unpacked(blas::Layout layout, const DenseDist &D, int64_t n_rows, int64_t n_cols, int64_t S_ro, int64_t S_co, T *buff, const State &seed) :project: RandBLAS @@ -127,4 +134,3 @@ The unifying (C++20) concepts .. doxygenconcept:: RandBLAS::SketchingOperator :project: RandBLAS - diff --git a/rtd/source/installation/index.rst b/rtd/source/installation/index.rst index 4f1b47ee..7dd7e84b 100644 --- a/rtd/source/installation/index.rst +++ b/rtd/source/installation/index.rst @@ -2,9 +2,10 @@ Installation ============ -RandBLAS is a header-only C++20 library with two required dependencies. One of these -dependencies (`Random123 `_) is header-only, -while the other (`BLAS++ `_) needs to be compiled. +RandBLAS is a header-only C++20 library with one required library dependency: +`BLAS++ `_, which needs to be compiled. +The counter-based random-number engines are headers included with RandBLAS; +there is no separate random-number package to install or configure. Having a compiled dependency makes setting up RandBLAS a little more complicated than setting up other header-only libraries. RandBLAS also has OpenMP and GoogleTest as @@ -63,14 +64,12 @@ MKL sparse features, pass ``-DRandBLAS_USE_MKL_SPARSE=OFF`` at configure time:: Everyone else ------------- -Strictly speaking, we only need three things to use RandBLAS in other projects. +Strictly speaking, we only need two things to use RandBLAS in other projects. 1. ``RandBLAS/config.h``, filled according to the instructions in ``RandBLAS/config.h.in``. -2. The locations of Random123 header files. - -3. The locations of the header files and compiled binary for BLAS++ (which will - referred to as "blaspp" when installed on your system). +2. The locations of the header files and compiled binary for BLAS++ (which will + be referred to as "blaspp" when installed on your system). If you have these things at hand, then compiling a RandBLAS-dependent program is just a matter of specifying standard compiler flags. diff --git a/rtd/source/tutorial/distributions.rst b/rtd/source/tutorial/distributions.rst index 1511039d..2805d4b0 100644 --- a/rtd/source/tutorial/distributions.rst +++ b/rtd/source/tutorial/distributions.rst @@ -132,7 +132,7 @@ narrow circumstances where one of these might be preferred in practice. We'll ex // Assume previous code defined integers (d1, d2, n) where 0 < d1 < d2 < n, // and "family" variable equal to ScalarDist::Gaussian or ScalarDist::Uniform, - // and a "state" variable of type RNGState. + // and a "state" variable satisfying GeneratorState. DenseDist D1(d1, n, family, Axis::Long); DenseDist D2(d2, n, family, Axis::Long); DenseSkOp S1(D1, state); diff --git a/rtd/source/tutorial/index.rst b/rtd/source/tutorial/index.rst index c7aea8af..d6cd52a6 100644 --- a/rtd/source/tutorial/index.rst +++ b/rtd/source/tutorial/index.rst @@ -31,12 +31,12 @@ RandBLAS, at a glance .. code:: c++ // step 1 - RandBLAS::RNGState state(); + RandBLAS::DefaultRNGState state{}; // step 2 RandBLAS::DenseDist D(10000, 50); RandBLAS::DenseSkOp S(D, state); // step 3 - double B* = new double[20000 * 50]; + double* B = new double[20000 * 50]; RandBLAS::sketch_general( blas::Layout::ColMajor, blas::Op::NoTrans, blas::Op::NoTrans, 20000, 50, 10000, diff --git a/rtd/source/tutorial/sampling_skops.rst b/rtd/source/tutorial/sampling_skops.rst index 1b21a64d..1cd97424 100644 --- a/rtd/source/tutorial/sampling_skops.rst +++ b/rtd/source/tutorial/sampling_skops.rst @@ -6,112 +6,129 @@ Sampling a sketching operator ****************************************************************************** -RandBLAS relies on counter-based random number generators (CBRNGs) from Random123. -A CBRNG returns a random number upon being called with two integer parameters: the *counter* and the *key*. -The time required for the CBRNG to return does not depend on either of these parameters. -A serial application can set the key at the outset of the program and never change it, while -parallel applications should use different keys across different threads. -Sequential calls to the CBRNG with a fixed key should use different values for the counter. +RandBLAS includes native, header-only counter-based random-number generators +(CBRNGs). A CBRNG is a stateless block function of a counter and a key. RandBLAS +sampling functions consume a state abstraction that binds such an engine to one +counter and one key. +This organization lets RandBLAS assign counter blocks to matrix coordinates +instead of relying on the order in which threads happen to run. Sampling a full +operator or one of its submatrices is therefore reproducible and independent of +the OpenMP thread count. -RandBLAS doesn't expose CBRNGs directly. Instead, it exposes an abstraction of -a CBRNG's state as defined in the :cpp:struct:`RandBLAS::RNGState` type. -RNGState objects are needed to construct sketching operators. .. _constructing_rng_states_tut: -Constructing RNGStates -====================== +Constructing RNG states +======================= -There are two ways to construct an RNGState from scratch: +Most applications should use :cpp:type:`RandBLAS::DefaultRNGState`, whose +engine is ``RandBLAS::rng::Philox<4, 32, 10>``: .. code:: c++ - RandBLAS::RNGState s1(); // key and counter are initialized to 0. - RandBLAS::RNGState s2(42); // key set to 42, counter set to 0. + RandBLAS::DefaultRNGState s1{}; // zero counter and default key + RandBLAS::DefaultRNGState s2{42}; // zero counter and key mapped from seed 42 -Note that in both cases the counter is initialized to zero. -This is important: you should never set the counter yourself! -If you want statistically independent runs of the same program, then you can start with different values for the key. +The scalar seed is mapped to a key by the engine's stable ``make_key`` function; +it should be treated as a seed rather than as the key's representation. Using +different seeds is the ordinary way to obtain independent program runs. -You can also construct an RNGState with a copy operation: +An RNG state is copyable: .. code:: c++ - RandBLAS::RNGState s3(s1); // s3 is a copy of s1. + RandBLAS::DefaultRNGState s3{s1}; Constructing your first sketching operator ========================================== -RandBLAS provides several constructors for the DenseSkOp and SparseSkOp classes. -However, the *recommended* constructors for these classes just accept two parameters: -a representation of a distribution (i.e., a DenseDist or a SparseDist) and an RNGState. +The recommended DenseSkOp and SparseSkOp constructors accept a distribution and +an RNG state. For example, the following code defines a +:math:`10000 \times 50` dense sketching operator whose entries are independent +standard-normal samples: -For example, the following code produces a :math:`10000 \times 50` dense sketching operator -whose entries are iid samples from the standard normal distribution. +.. code:: c++ - .. code:: c++ + RandBLAS::DefaultRNGState state{}; + RandBLAS::DenseDist dist(10000, 50); + RandBLAS::DenseSkOp S(dist, state); + // state is copied into S.seed_state. Sampling happens only when needed. - RandBLAS::RNGState my_state(); - RandBLAS::DenseDist my_dist(10000, 50); - RandBLAS::DenseSkOp S(my_dist, my_state); - // my_state is stored as a constant value S.seed_state. - // S.seed_state will be accessed by RandBLAS' random sampling - // functions behind the scenes, only when needed. - -We note that the numerical precision of the sketching operator must be specified with a template parameter; -the entries of the sketching operator are defined by sampling in single precision and then -casting the sample to double if needed. +The numerical precision of the sketching operator is a template parameter. +Entries are sampled in single precision and cast to double when needed. -Formal API docs for the recommended constructors can be found :ref:`here ` and :ref:`here `. +Formal API docs for the recommended constructors can be found +:ref:`here ` and +:ref:`here `. Constructing your :math:`N^{\text{th}}` sketching operator, for :math:`N > 1` ============================================================================== -Suppose you have an application that requires two statistically independent dense -sketching operators, :math:`\texttt{S1}` and :math:`\texttt{S2}`, each of size -:math:`10000 \times 50`. How should you get your hands on these objects? +Constructing two operators from the same distribution and state defines the +same mathematical operator, not two independent ones: + +.. code:: c++ -.. warning:: - If you try to construct those sketching operators as follows ... + RandBLAS::DefaultRNGState state{}; + RandBLAS::DenseDist dist(10000, 50); + RandBLAS::DenseSkOp S1(dist, state); + RandBLAS::DenseSkOp S2(dist, state); // S2 is the same as S1 - .. code:: c++ +Use the first operator's ``next_state`` for the second operator: - RandBLAS::RNGState my_state(); - RandBLAS::DenseDist my_dist(10000, 50); - RandBLAS::DenseSkOp S1(my_dist, my_state); - RandBLAS::DenseSkOp S2(my_dist, my_state); +.. code:: c++ - *then your results would be invalid! Far from being independent,* :math:`\texttt{S1}` - *and* :math:`\texttt{S2}` *would be equal from a mathematical perspective.* + RandBLAS::DefaultRNGState state{}; + RandBLAS::DenseDist dist(10000, 50); + RandBLAS::DenseSkOp S1(dist, state); + RandBLAS::DenseSkOp S2(dist, S1.next_state); -One correct approach is to then call the constructor for :math:`\texttt{S2}` -using :math:`\texttt{S1.next_state}` as its RNGState argument: +Alternatively, start with two states constructed from different seeds: - .. code:: c++ +.. code:: c++ - RandBLAS::RNGState my_state(); - RandBLAS::DenseDist my_dist(10000, 50); - RandBLAS::DenseSkOp S1(my_dist, my_state); - // ^ Defines S1 from a mathematical perspective. Computes S1.next_state, - // but otherwise performs no work. - RandBLAS::DenseSkOp S2(my_dist, S1.next_state); + RandBLAS::DefaultRNGState state1{19}; + RandBLAS::DefaultRNGState state2{93}; + RandBLAS::DenseDist dist(10000, 50); + RandBLAS::DenseSkOp S1(dist, state1); + RandBLAS::DenseSkOp S2(dist, state2); -Another valid approach is to declare two RNGState objects from the beginning using -different keys, as in the following code: - .. code:: c++ +Engine and state details +======================== - RandBLAS::RNGState my_state1(19); - // ^ An RNGState with zero'd counter and key initialized to 19. - RandBLAS::RNGState my_state2(93); - // ^ An RNGState with zero'd counter and key initialized to 93. - RandBLAS::DenseDist my_dist(10000, 50); - RandBLAS::DenseSkOp S1(my_dist, my_state1); - RandBLAS::DenseSkOp S2(my_dist, my_state2); - // ^ S1 and S2 are defined only from a mathematical perspective. - // No real work is performed here. +Users who need direct block access can name the engine and state explicitly: +.. code:: c++ + using Engine = RandBLAS::rng::Philox<4, 32, 10>; + using State = RandBLAS::RNGState; + + State state{1234}; + Engine::res_t block{}; + state.generate(block); // writes every element of the output-only array + state.advance(1); // add one to the engine's multiword counter + +The engine provides ``ctr_t``, ``key_t``, and ``res_t`` aliases. The concrete +``RNGState`` adapter exposes the same aliases and public ``counter``, ``key``, +and ``engine`` data. Generic samplers require only the ``GeneratorState`` +operations, so custom states need not use that representation. Generation does +not mutate the state: advancing by one block is always explicit. The default +Philox engine produces exactly the same integer blocks as Philox4x32-10 in +Random123 for the same counter and key. Philox is a statistical generator, not +a cryptographic random-number generator. + +``RandBLAS::rng::RepackedOutput`` can expose each result word as narrower chunks +without changing the block boundary. Chunks are least-significant first within +each source word: ``0xAABBCCDD`` becomes ``{0xCCDD, 0xAABB}`` with 16-bit words +and ``{0xDD, 0xCC, 0xBB, 0xAA}`` with 8-bit words, regardless of host byte order. +Current RandBLAS samplers accept native 32- or 64-bit result words; repacked +8- and 16-bit results are provided for direct use and future sampler work. + +The dense Gaussian transform calls the host ``sin``, ``cos``, ``log``, and +``sqrt`` functions. Its final bits can therefore vary across math libraries, +compilers, or architectures even though the underlying integer stream and block +assignment are fixed. diff --git a/rtd/source/tutorial/sketch_updates.rst b/rtd/source/tutorial/sketch_updates.rst index fa05bba5..a79dab9e 100644 --- a/rtd/source/tutorial/sketch_updates.rst +++ b/rtd/source/tutorial/sketch_updates.rst @@ -207,7 +207,7 @@ Implementation // Since d < m and we're short-axis major, the columns of matrices sampled from // D1 or D1 will be sampled i.i.d. from some distribution on R^d. - auto S1 = D1.sample( seed_state ); // seed_state is some RNGState. + auto S1 = D1.sample( seed_state ); // seed_state satisfies GeneratorState. auto S = D.sample( seed_state ); // With these definitions, S1 is *always* equal to the first m columns of S. // We recover S2 by working implicitly with the trailing k columns of S. @@ -271,8 +271,7 @@ Implementation // Since n > d and we're short-axis major, the rows of matrices sampled from // D1 or D1 will be sampled i.i.d. from some distribution on R^d. - auto S1 = D1.sample( seed_state ); // seed_state is some RNGState. + auto S1 = D1.sample( seed_state ); // seed_state satisfies GeneratorState. auto S = D.sample( seed_state ); // With these definitions, S1 is *always* equal to the first m rows of S. // We recover S2 by working implicitly with the last k rows of S. - diff --git a/rtd/source/updates/index.rst b/rtd/source/updates/index.rst index d66218aa..e308e056 100644 --- a/rtd/source/updates/index.rst +++ b/rtd/source/updates/index.rst @@ -122,7 +122,7 @@ This makes it possible to sample from a templated .. code:: c++ - RNGState seed_state(8675309); + DefaultRNGState seed_state(8675309); auto S = D.sample(seed_state); In RandBLAS 1.0 it was necessary to construct a sketching operator by calling diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b60dc84a..45781520 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -60,16 +60,21 @@ if (GTest_FOUND) ##################################################################### set(STAT_SOURCES - basic_rng/test_r123.cc + basic_rng/test_philox.cc + basic_rng/test_repacked_output.cc + basic_rng/test_distributions.cc + basic_rng/test_rng_state.cc basic_rng/test_discrete.cc basic_rng/test_continuous.cc basic_rng/test_distortion.cc + basic_rng/test_sampler_regression.cc + basic_rng/test_word_array.cc ) add_executable(stat_tests ${STAT_SOURCES}) target_link_libraries(stat_tests RandBLAS GTest::GTest GTest::Main) # Bake the KAT vector path in at compile time so the test is cwd-independent. target_compile_definitions(stat_tests PRIVATE - KAT_VECTORS_PATH="${CMAKE_CURRENT_SOURCE_DIR}/basic_rng/r123_kat_vectors.txt") + PHILOX_KAT_VECTORS_PATH="${CMAKE_CURRENT_SOURCE_DIR}/basic_rng/philox_kat_vectors.txt") randblas_stage_runtime_dlls(stat_tests) gtest_discover_tests(stat_tests) @@ -79,7 +84,12 @@ if (GTest_FOUND) # ##################################################################### - set(META_SOURCES meta/test_lapack_like.cc meta/test_sparse_data_generators.cc meta/test_comparison.cc) + set(META_SOURCES + meta/test_lapack_like.cc + meta/test_sparse_data_generators.cc + meta/test_comparison.cc + meta/test_rng_stream.cc + ) add_executable(meta_tests ${META_SOURCES}) target_link_libraries(meta_tests RandBLAS GTest::GTest GTest::Main) randblas_stage_runtime_dlls(meta_tests) diff --git a/test/DevNotes.md b/test/DevNotes.md index b2c121d6..994eb3be 100644 --- a/test/DevNotes.md +++ b/test/DevNotes.md @@ -29,18 +29,41 @@ Tests of wrapper functions. * sketch_symmetric. It reduces to the same sketch_general no matter the type of the sketching operator. -### test_basic_rng +### basic_rng Relies on RandBLAS/testing/stats.hh. - * test_r123.cc has deterministic tests for Random123. The tests compare generated values - to reference values computed ahead of time. The tests are __extremely__ messy, since they're - adapted from tests in the official Random123 repository, and Random123 needs to handle a far wider - range of compilers and languages than we assume for RandBLAS. + * `test_philox.cc` validates the native Philox implementation against 204 + static vectors generated offline from a pinned Random123 revision. The test + suite never locates or executes Random123. + * `test_word_array.cc` covers modular carry propagation and wraparound. + * `test_rng_state.cc` checks engine/state concepts, scalar seed mapping, + output-only block generation, explicit advancement, public state data, + and compatibility with an engine whose counter representation is opaque. + * `test_repacked_output.cc` checks direct and nested repacking, including + least-significant-chunk-first ordering. + * `test_distributions.cc` checks the native integer-to-floating transforms, + Box--Muller reference values, endpoints, policy loops, and word assignment. + * `test_sampler_regression.cc` protects the pre-migration dense and sparse + streams. Sparse outputs are bitwise exact; dense comparisons use the narrow + floating-point tolerance required for host math-library differences. * test_discrete.cc includes statistical tests for sampling from an index set with or without replacement. +Sampler tests elsewhere cover state advancement, full/submatrix agreement, and +OpenMP thread-count independence. Downstream-package and example builds verify +that no external random-number package is required. + +### RNG stream + +`RandBLAS/testing/rng.hh` contains the test-only `detail::RNGStream` adapter. +It turns fixed result blocks into a sequential word stream for random sparse +test-matrix generation and supplies the uniform, Gaussian, and geometric draws +needed there. Loading a block advances its held state immediately; unread lanes +remain in its local buffer. Production RandBLAS sampling remains +coordinate-addressed and does not use this sequential adapter. + # OLD diff --git a/test/basic_rng/benchmark_speed.cc b/test/basic_rng/benchmark_speed.cc index 5f08396d..5588d11d 100644 --- a/test/basic_rng/benchmark_speed.cc +++ b/test/basic_rng/benchmark_speed.cc @@ -58,12 +58,12 @@ std::ostream &operator<<(std::ostream &os, std::vector &v) -template +template auto run_test(RandBLAS::DenseDist D, T *mat) { auto t0 = std::chrono::high_resolution_clock::now(); - RNGState seed; - RandBLAS::dense::fill_dense_submat_impl(D.n_cols, mat, D.n_rows, D.n_cols, 0, seed); + State seed; + RandBLAS::dense::fill_dense_submat_impl(D.n_cols, mat, D.n_rows, D.n_cols, 0, seed); auto t1 = std::chrono::high_resolution_clock::now(); return (t1 - t0).count(); } @@ -74,8 +74,8 @@ int main(int argc, char **argv) (void) argc; using T = float; - using RNG = r123::Philox4x32; - using OP = r123ext::uneg11; + using State = RandBLAS::DefaultRNGState; + using OP = RandBLAS::rng::uneg11; int64_t m = atoi(argv[1]); int64_t n = atoi(argv[2]); @@ -84,9 +84,9 @@ int main(int argc, char **argv) std::vector mat(d); - auto dt = run_test(dist, mat.data()); + auto dt = run_test(dist, mat.data()); - std::cerr << "[" << typeid(RNG).name() << ", " + std::cerr << "[" << typeid(State).name() << ", " << typeid(OP).name() << "] dt = " << dt << std::endl; if (d < 100) @@ -94,4 +94,3 @@ int main(int argc, char **argv) return 0; } - diff --git a/test/basic_rng/philox_kat_vectors.txt b/test/basic_rng/philox_kat_vectors.txt new file mode 100644 index 00000000..eff1443e --- /dev/null +++ b/test/basic_rng/philox_kat_vectors.txt @@ -0,0 +1,240 @@ +# Copyright, 2026. See LICENSE for copyright holder information. +# Copyright 2010-2011, D. E. Shaw Research. +# 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 D. E. Shaw Research nor the names of its contributors +# may be used to endorse or promote products derived from this software +# without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# Generated offline from Random123 commit +# 9545ff6413f258be2f04c1d319d99aaef7521150. +# Each family/round has zero, published nonzero, and carry-heavy inputs. +# The carry-heavy case repeats {maximum word, 1, high bit} over counter then key. +# +# family round counter... key... result... +philox2x32 0 00000000 00000000 00000000 00000000 00000000 +philox2x32 0 243f6a88 85a308d3 13198a2e 243f6a88 85a308d3 +philox2x32 0 ffffffff 00000001 80000000 ffffffff 00000001 +philox2x32 1 00000000 00000000 00000000 00000000 00000000 +philox2x32 1 243f6a88 85a308d3 13198a2e 8b72d2a0 f0133418 +philox2x32 1 ffffffff 00000001 80000000 5256d193 2da92e6d +philox2x32 2 00000000 00000000 00000000 9e3779b9 00000000 +philox2x32 2 243f6a88 85a308d3 13198a2e 33d14c7d d2a391e0 +philox2x32 2 ffffffff 00000001 80000000 703973a9 26555a69 +philox2x32 3 00000000 00000000 00000000 bd91d970 ca60ee3b +philox2x32 3 243f6a88 85a308d3 13198a2e b7b8af39 dc9ef8c7 +philox2x32 3 ffffffff 00000001 80000000 c60e9917 d530630b +philox2x32 4 00000000 00000000 00000000 8b076d9d d8e44b50 +philox2x32 4 243f6a88 85a308d3 13198a2e a7d2c1b0 e54026bb +philox2x32 4 ffffffff 00000001 80000000 2d2d2abf f117af35 +philox2x32 5 00000000 00000000 00000000 d202938e 827f1e27 +philox2x32 5 243f6a88 85a308d3 13198a2e e7549abb df45e810 +philox2x32 5 ffffffff 00000001 80000000 2cd416e9 daff7aad +philox2x32 6 00000000 00000000 00000000 39e72b6e bfa5a88a +philox2x32 6 243f6a88 85a308d3 13198a2e 4b7ad861 48b98461 +philox2x32 6 ffffffff 00000001 80000000 693f348d b8c360cb +philox2x32 7 00000000 00000000 00000000 257a3673 cd26be2a +philox2x32 7 243f6a88 85a308d3 13198a2e bedbbe6b e4c770b3 +philox2x32 7 ffffffff 00000001 80000000 dbf6301d 5a9049f7 +philox2x32 8 00000000 00000000 00000000 80681c43 ec432709 +philox2x32 8 243f6a88 85a308d3 13198a2e 1e8a52ca bc95b271 +philox2x32 8 ffffffff 00000001 80000000 3dae91c6 3f634da7 +philox2x32 9 00000000 00000000 00000000 74781986 2f5ced79 +philox2x32 9 243f6a88 85a308d3 13198a2e a15736e6 fae073fe +philox2x32 9 ffffffff 00000001 80000000 7c769ed7 54c45ab2 +philox2x32 10 00000000 00000000 00000000 ff1dae59 6cd10df2 +philox2x32 10 243f6a88 85a308d3 13198a2e dd7ce038 f62a4c12 +philox2x32 10 ffffffff 00000001 80000000 3d749939 0a04bc75 +philox2x32 11 00000000 00000000 00000000 9367111f 4d47c61b +philox2x32 11 243f6a88 85a308d3 13198a2e 02958a49 64927828 +philox2x32 11 ffffffff 00000001 80000000 96500064 b63184bb +philox2x32 12 00000000 00000000 00000000 f83975be f79323cd +philox2x32 12 243f6a88 85a308d3 13198a2e b9f636c0 e84a00eb +philox2x32 12 ffffffff 00000001 80000000 81d32477 c9d9dd6c +philox2x32 13 00000000 00000000 00000000 56f9c679 dbf2ba1a +philox2x32 13 243f6a88 85a308d3 13198a2e 0d3237ed af943040 +philox2x32 13 ffffffff 00000001 80000000 49eb52df cafd1755 +philox2x32 14 00000000 00000000 00000000 9455f794 cb1bc07b +philox2x32 14 243f6a88 85a308d3 13198a2e bea9235f c71c9a17 +philox2x32 14 ffffffff 00000001 80000000 7e9022bc 55b2a50d +philox2x32 15 00000000 00000000 00000000 15f3bb02 f034fdfc +philox2x32 15 243f6a88 85a308d3 13198a2e e199cfb0 60e9de8d +philox2x32 15 ffffffff 00000001 80000000 15472595 3c476df4 +philox2x32 16 00000000 00000000 00000000 a77dbd6a 0d4d0426 +philox2x32 16 243f6a88 85a308d3 13198a2e 81ecdc32 4300f210 +philox2x32 16 ffffffff 00000001 80000000 e87cd805 2797398f +philox4x32 0 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 +philox4x32 0 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 243f6a88 85a308d3 13198a2e 03707344 +philox4x32 0 ffffffff 00000001 80000000 ffffffff 00000001 80000000 ffffffff 00000001 80000000 ffffffff +philox4x32 1 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 +philox4x32 1 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 2efd7704 ad2d4ba2 3728c377 b37e0218 +philox4x32 1 ffffffff 00000001 80000000 ffffffff 00000001 80000000 66cf46ab 80000000 adaee0ad 2daee0ad +philox4x32 2 00000000 00000000 00000000 00000000 00000000 00000000 9e3779b9 00000000 bb67ae85 00000000 +philox4x32 2 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 c320345a 20f4f871 70e22dde dddb124c +philox4x32 2 ffffffff 00000001 80000000 ffffffff 00000001 80000000 95b7e207 8bf3a3cb 42bfd20a 10dd9e71 +philox4x32 3 00000000 00000000 00000000 00000000 00000000 00000000 aae8eb44 02719033 f734f9c9 4942ddfb +philox4x32 3 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 9a27db43 0524dc72 ddfbdc94 8e41df2e +philox4x32 3 ffffffff 00000001 80000000 ffffffff 00000001 80000000 8201a2d7 fa0be366 9d128b91 11342145 +philox4x32 4 00000000 00000000 00000000 00000000 00000000 00000000 1e597a2c 4fa2984f f71cefed d685830c +philox4x32 4 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 c9c7529c d06d7a4c ab326927 24ac33b9 +philox4x32 4 ffffffff 00000001 80000000 ffffffff 00000001 80000000 5e84afa7 d8aa4b47 c9cda0cb 3d46d4b5 +philox4x32 5 00000000 00000000 00000000 00000000 00000000 00000000 f10446c0 c841128b 23f43d26 6cb9f044 +philox4x32 5 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 450b009c 131e3741 9654aad8 a9fcac94 +philox4x32 5 ffffffff 00000001 80000000 ffffffff 00000001 80000000 02610bb4 a2bb73fd 1d7ea502 48212c25 +philox4x32 6 00000000 00000000 00000000 00000000 00000000 00000000 c3b4ab98 e922b5ea 03be793c 8bb43040 +philox4x32 6 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 d0be47ee 222f0768 43e1d190 e1001694 +philox4x32 6 ffffffff 00000001 80000000 ffffffff 00000001 80000000 a21ebdf6 38282dae 60d31824 e1d5975c +philox4x32 7 00000000 00000000 00000000 00000000 00000000 00000000 5f6fb709 0d893f64 4f121f81 4f730a48 +philox4x32 7 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 4dfccaba 190a87f0 c47362ba b6b5242a +philox4x32 7 ffffffff 00000001 80000000 ffffffff 00000001 80000000 c0a1f5eb e641083c 808b1726 2fce60c2 +philox4x32 8 00000000 00000000 00000000 00000000 00000000 00000000 618f177a 9920c1d7 1ec12dc0 c43b6eeb +philox4x32 8 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 734d1a7d 3162ff36 bfd3c78a 9e5c404e +philox4x32 8 ffffffff 00000001 80000000 ffffffff 00000001 80000000 d2fba276 0278cbea 2e5a43ea a6a53031 +philox4x32 9 00000000 00000000 00000000 00000000 00000000 00000000 7028b1c8 d6594c40 4f2051bb 76d6628e +philox4x32 9 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 3eb49fbb 7f0bd1e6 c4392c8d a3c1b987 +philox4x32 9 ffffffff 00000001 80000000 ffffffff 00000001 80000000 d6f9fb8c c480f686 50cd02ad 1397f642 +philox4x32 10 00000000 00000000 00000000 00000000 00000000 00000000 6627e8d5 e169c58d bc57ac4c 9b00dbd8 +philox4x32 10 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 d16cfe09 94fdcceb 5001e420 24126ea1 +philox4x32 10 ffffffff 00000001 80000000 ffffffff 00000001 80000000 0b95874c 8feb31cb b5affb50 67ce8264 +philox4x32 11 00000000 00000000 00000000 00000000 00000000 00000000 5805dfa6 7e9969d4 9ae116f5 5987480f +philox4x32 11 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 068d24b3 a70a26e0 f3b3a8dc dbf273eb +philox4x32 11 ffffffff 00000001 80000000 ffffffff 00000001 80000000 302f83db d7997830 bc4607e1 03e911a4 +philox4x32 12 00000000 00000000 00000000 00000000 00000000 00000000 ce9d78f7 f859be43 1ca381af 2f829cd2 +philox4x32 12 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 14dc9852 29d78ec4 e98019fc a9d79309 +philox4x32 12 ffffffff 00000001 80000000 ffffffff 00000001 80000000 8cc3edbe 98019a77 a90b5104 beaa4501 +philox4x32 13 00000000 00000000 00000000 00000000 00000000 00000000 85c0a17d f7007579 4ee014bb f2db2115 +philox4x32 13 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 9cf950c3 be60a0a4 4a8f6249 efe95096 +philox4x32 13 ffffffff 00000001 80000000 ffffffff 00000001 80000000 755eecea e2efbc5c 85d33f00 726e169a +philox4x32 14 00000000 00000000 00000000 00000000 00000000 00000000 c08b116e c3fc0a8d 1b7a9b1d 538f7e87 +philox4x32 14 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 2958397e 6dea9bcf c2fc1c50 994fcc39 +philox4x32 14 ffffffff 00000001 80000000 ffffffff 00000001 80000000 8143998c b67d6900 1640c5e6 768125de +philox4x32 15 00000000 00000000 00000000 00000000 00000000 00000000 72e6919e d2ffafdb f20be21d 0400f8aa +philox4x32 15 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 ba67de63 a8a5af30 d1f2080d a36ee5da +philox4x32 15 ffffffff 00000001 80000000 ffffffff 00000001 80000000 00aa67a3 8ef6ef2a a318c4b1 85cebc64 +philox4x32 16 00000000 00000000 00000000 00000000 00000000 00000000 55d6e305 9479d0db a1764d17 db61583a +philox4x32 16 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 e94c2174 08b5e56b 1ef8c858 96ee1719 +philox4x32 16 ffffffff 00000001 80000000 ffffffff 00000001 80000000 49491625 33fe5527 fe567aac 5a5f56d9 +philox2x64 0 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 +philox2x64 0 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 243f6a8885a308d3 13198a2e03707344 +philox2x64 0 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 +philox2x64 1 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 +philox2x64 1 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 aac54a573e83c7cf 891461e0c732bb29 +philox2x64 1 ffffffffffffffff 0000000000000001 8000000000000000 52b74407b1ce6e93 2d48bbf84e31916d +philox2x64 2 0000000000000000 0000000000000000 0000000000000000 9e3779b97f4a7c15 0000000000000000 +philox2x64 2 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 47c4c75e50bd0ac8 534f0acf6d18addd +philox2x64 2 ffffffffffffffff 0000000000000001 8000000000000000 776a4e8a2a982c8a a1b13f676a56a869 +philox2x64 3 0000000000000000 0000000000000000 0000000000000000 be5436ff55565a24 1571808376fc460f +philox2x64 3 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 8825ed120afabc16 986119e48c1f20d8 +philox2x64 3 ffffffffffffffff 0000000000000001 8000000000000000 7f957d32dba8bb96 9b0fd4a91990df3e +philox2x64 4 0000000000000000 0000000000000000 0000000000000000 537eb84bd082ec62 827eb21488493aac +philox2x64 4 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 96de2f69441856fb 1e4c107c9a7f74a2 +philox2x64 4 ffffffffffffffff 0000000000000001 8000000000000000 a8aa4ff2a6cf3033 5f1a232a4b322b22 +philox2x64 5 0000000000000000 0000000000000000 0000000000000000 be1aec056fc65891 f58d3498f59bd846 +philox2x64 5 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 7e854cbb77267bf5 6ff0214b1853cc21 +philox2x64 5 ffffffffffffffff 0000000000000001 8000000000000000 2d13b1e9344e41b2 05aa8473dbb89749 +philox2x64 6 0000000000000000 0000000000000000 0000000000000000 7ee2796782e4de12 6921e1f4eea12943 +philox2x64 6 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 bccd525069e680f7 5adc9137188273af +philox2x64 6 ffffffffffffffff 0000000000000001 8000000000000000 b7a589ecb77ba74d f168a3a598663536 +philox2x64 7 0000000000000000 0000000000000000 0000000000000000 b41da69fbfefc666 511e9ce1a5534056 +philox2x64 7 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 98ed1534392bf372 67528b1568882fd5 +philox2x64 7 ffffffffffffffff 0000000000000001 8000000000000000 530d3aba27a914f3 db87150c5ada2737 +philox2x64 8 0000000000000000 0000000000000000 0000000000000000 96db8b4daf8a8498 220c121c6a02c092 +philox2x64 8 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 ed00f2a35daf258f a228f25b6e93c676 +philox2x64 8 ffffffffffffffff 0000000000000001 8000000000000000 4c5f0fad93937b68 a5f1a61f4ba17189 +philox2x64 9 0000000000000000 0000000000000000 0000000000000000 af9bc89ab801df97 d57866551bd37348 +philox2x64 9 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 f4f9860798d8a54f f00604e125c8031d +philox2x64 9 ffffffffffffffff 0000000000000001 8000000000000000 ea96ccc89200c0cc 11304eed52668cb8 +philox2x64 10 0000000000000000 0000000000000000 0000000000000000 ca00a0459843d731 66c24222c9a845b5 +philox2x64 10 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 0a5e742c2997341c b0f883d38000de5d +philox2x64 10 ffffffffffffffff 0000000000000001 8000000000000000 dfd4ad482d619b3d 51fd3d64596e5d24 +philox2x64 11 0000000000000000 0000000000000000 0000000000000000 eead9c85d4cb3288 8e8c441a9fd99f23 +philox2x64 11 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 6a439d200edc821a 13708eb224bef414 +philox2x64 11 ffffffffffffffff 0000000000000001 8000000000000000 47eb4ef2eda67133 6e8bf77d17d65a07 +philox2x64 12 0000000000000000 0000000000000000 0000000000000000 869b390f11b9eb23 9cf69d7331d47418 +philox2x64 12 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 346c74de4013b0bd fb95e1e5f371e0ee +philox2x64 12 ffffffffffffffff 0000000000000001 8000000000000000 19db4e900f10c47b 5ce9ac345540ea49 +philox2x64 13 0000000000000000 0000000000000000 0000000000000000 98a4984bab9d3883 fd97c66879f50f19 +philox2x64 13 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 de117265efb790fd c63ad15d2655b287 +philox2x64 13 ffffffffffffffff 0000000000000001 8000000000000000 a3387193ff1451a5 1e99c5e60107aca1 +philox2x64 14 0000000000000000 0000000000000000 0000000000000000 88e2a32430ea8f1b 8e741b5c25f9bd39 +philox2x64 14 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 dc29fe8332af8bca 0b54796a1f4af747 +philox2x64 14 ffffffffffffffff 0000000000000001 8000000000000000 1011ce6909d1f00f c8dc2de5fb85c7bf +philox2x64 15 0000000000000000 0000000000000000 0000000000000000 59d75b4e4c71fc48 5fe8da3ab2e7c681 +philox2x64 15 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 f57d966355f390a0 a06c288caf6a10fe +philox2x64 15 ffffffffffffffff 0000000000000001 8000000000000000 e2ee9990e00902ec 8a7542db63c54a9d +philox2x64 16 0000000000000000 0000000000000000 0000000000000000 535a025ca08249a1 815a6beda9cacd58 +philox2x64 16 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 8335b1ec8bf1e0ab 4e8c0f5603c0cbe0 +philox2x64 16 ffffffffffffffff 0000000000000001 8000000000000000 f5ff4dd3c7de2026 c632b52e0a561584 +philox4x64 0 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 +philox4x64 0 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 +philox4x64 0 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff +philox4x64 1 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 +philox4x64 1 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c d798eedd15d15f3f c378d0ff4808bdb0 aba658fb081df253 c219bc7795fb1529 +philox4x64 1 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 652d4131ca8908ab 8000000000000000 ad18b8f11eb3936d 2d18b8f11eb3936d +philox4x64 2 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 9e3779b97f4a7c15 0000000000000000 bb67ae8584caa73b 0000000000000000 +philox4x64 2 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c a78ab97e6b327c30 c6a1c518181add35 0a3b93ebfbf63174 6e1c9b78811c452d +philox4x64 2 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 96e5cd70550e7370 b622e7912a7b570b 4525911280d2b11b 7b7ae8affa1c1e31 +philox4x64 3 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 a84ce2eebc058cfe 2a891897af1fc00d f4972cd89c434076 3c970488cf5a1c0f +philox4x64 3 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c 4f2070a15e9b9156 8f69fae1f81b826c d138707af41f5597 352d2593e8a18f90 +philox4x64 3 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 bce4020647726707 57966b2df042fb2d f1e56061f4f25407 0399af42a53f8950 +philox4x64 4 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 317ab6a74bdd23de bdc29adec86ebe1a 84071127f7ef4e38 ad9de0dd68131dda +philox4x64 4 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c 35c7f4dbd5151c76 e10c9268f2f31d51 84964f6cb2368a5a 659e5e32b02cbc62 +philox4x64 4 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 320485e44ca5640d 13a69335516d0561 2a330b071aea889d 1e9625a5073c1d05 +philox4x64 5 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 ad433d7fb14414ec e59be30e3e754d08 68c007808d235e2d 6639e0629f05407a +philox4x64 5 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c 37c7e1af5f896f3d 28b6b0f04c0cfe96 e223e3e0246aaadf 212f70dc6c291fc2 +philox4x64 5 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 4a2047f37de2e0bc 0f832cebb4d0da5b 5a3c7a7f5551e0a3 09bf642cea09ef77 +philox4x64 6 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 a0420d5917a1ff76 db43ddc3206ffe4b 418238cffd2a6360 64be3e9d24fb9384 +philox4x64 6 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c c64b77eed398428e ffb44ac84446e0c9 6b81e44c8fcc0ee5 3f6dbbe979f49c07 +philox4x64 6 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 5fc5d395929e8275 9257627a19342a65 1da87e0e911300c0 5d323f14f3e65bf4 +philox4x64 7 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 5dc8ee6268ec62cd 139bc570b6c125a0 84d6deb4fb65f49e aff7583376d378c2 +philox4x64 7 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c 513a366704edf755 f05d9924c07044d3 bef2cb9cbea74c6c 8db948de4caa1f8a +philox4x64 7 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 306ad7ec0b838fd7 2052db9b5e020140 f7bae59de25f68a5 668d05db9fca452f +philox4x64 8 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 291f138f2e410d0d 130d8dd85e009fb2 fd611071d5db6999 7fb091ff7b4737b7 +philox4x64 8 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c fe1efb9aa5af16c9 557ec7f0228624b4 11784618ec1b1da9 119a69e95b38e1cf +philox4x64 8 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 b0078099b1487f64 9daa391ca8f98513 debba285d7862664 b7bb1ca2c80e4c75 +philox4x64 9 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 2afe6c34d161a4b9 f0d16bf1b9560bff 856d415341667984 6229dc394db7fa77 +philox4x64 9 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c 6e523c3c3ffdb642 36aa26ac7d114d6f 59514df1419547b8 abbbc8b26dd2e16b +philox4x64 9 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 dc1f3971f9423b32 5af0fb3f0e2baffc 7d83479180404bca 4988bc00180f566c +philox4x64 10 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 16554d9eca36314c db20fe9d672d0fdc d7e772cee186176b 7e68b68aec7ba23b +philox4x64 10 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c a528f45403e61d95 38c72dbd566e9788 a5a1610e72fd18b5 57bd43b5e52b7fe6 +philox4x64 10 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 b6364a3702395472 5f9c4fb6c1162ba6 ea75e4d47bbd90b2 3a3b620696d915b6 +philox4x64 11 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 5fa2ca4e80f3a9d0 32740cb878a6105d 3e0241658290f26a 54ceff0f687a5ea4 +philox4x64 11 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c c97e2b358122e443 bd49ceb7f8616a83 cfccb58ad9a5d4d4 7233175a0ed9d88f +philox4x64 11 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 c8e565e51757575b b59641cf6a8bfe7e 7e2a9540046cebb8 a6311affba639576 +philox4x64 12 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 cf1585ba5a73898d dda91571bccf6c06 177370d77b04c2c1 2293b371114e4270 +philox4x64 12 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c 08836d67628ae74b 46cdb3285565680c 1c0479e9722b8726 73bb73d886395679 +philox4x64 12 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 1a4e4fc8bd01c025 89f8460f679b5388 8ec4c91b25afd3aa b55ade142a058d41 +philox4x64 13 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 a5b9c0513ade557d bc1166b9261f0097 40d520aaa2efb5aa 404b6f12713b77f7 +philox4x64 13 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c ff2915ba59344e21 044ddb554f7073ea f388667a9174279f d258398fc49a7411 +philox4x64 13 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 93b8322b0642ab70 f15251074fc338c6 e82af13ffbc7ef9a d01ceb6cd50cf13f +philox4x64 14 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 87ff5f6bc386183a 8b25fd24027b06c6 4c80ab8a9e9e4df2 934fee74d7d7d2c7 +philox4x64 14 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c 89cb31de5d9e2e75 fd5fdd0779490609 42f642a1858173c1 09c9f8a37ccac8f3 +philox4x64 14 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 4e0088d9527207c8 9aca2b8852afa756 adedbaf650b58d8e a476d9f10ddbb150 +philox4x64 15 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 1055c68e3d267d82 5f10d3b9bffd8f3e dcee3c4caf7d62a7 311f91a40370614e +philox4x64 15 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c 2582ea4dad9b3729 27dedfde70402797 85b3213ece9ca2da fa1f415f1d2a092f +philox4x64 15 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 b4b9901cd3d2c662 9f84172316158942 5b9fbad83e92e6e1 dde1db338a22d7d8 +philox4x64 16 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 b4f2e5bb6a75ef84 c4cbdb8819e79dc1 c779b6b9510afe6f 7035ec51dea4e9a6 +philox4x64 16 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c c4180154e673d2e3 efaad82af35fd216 5d9fa640404a49a6 8c72a8a59491f88b +philox4x64 16 ffffffffffffffff 0000000000000001 8000000000000000 ffffffffffffffff 0000000000000001 8000000000000000 92a858bc712fa617 3f784d5d36136777 32119f33ea0dee0b 9f0f143594d14246 diff --git a/test/basic_rng/r123_kat_vectors.txt b/test/basic_rng/r123_kat_vectors.txt deleted file mode 100644 index 9398f32c..00000000 --- a/test/basic_rng/r123_kat_vectors.txt +++ /dev/null @@ -1,75 +0,0 @@ -# This file is copied from https://github.com/DEShawResearch/random123/blob/main/tests/kat_vectors -# -# For each generator, we test: gen(0, 0), gen(fff, fff) and gen(ctr=digits_of_pi, key=more_digits_of_pi). -# Ignoring endianness, these are the first few hexdigits of pi: -# 243F6A88 85A308D3 13198A2E 03707344 A4093822 299F31D0 082EFA98 EC4E6C89 452821E6 38D01377 BE5466CF 34E90C6C C0AC29B7 C97C50DD 3F84D5B5 B5470917 9216D5D9 8979FB1BD -# -#nameNxW R CTR KEY EXPECTED -# -philox2x32 7 00000000 00000000 00000000 257a3673 cd26be2a -philox2x32 7 ffffffff ffffffff ffffffff ab302c4d 3dc9d239 -philox2x32 7 243f6a88 85a308d3 13198a2e bedbbe6b e4c770b3 -philox2x32 10 00000000 00000000 00000000 ff1dae59 6cd10df2 -philox2x32 10 ffffffff ffffffff ffffffff 2c3f628b ab4fd7ad -philox2x32 10 243f6a88 85a308d3 13198a2e dd7ce038 f62a4c12 -# -philox4x32 7 00000000 00000000 00000000 00000000 00000000 00000000 5f6fb709 0d893f64 4f121f81 4f730a48 -philox4x32 7 ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff 5207ddc2 45165e59 4d8ee751 8c52f662 -philox4x32 7 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 4dfccaba 190a87f0 c47362ba b6b5242a -philox4x32 10 00000000 00000000 00000000 00000000 00000000 00000000 6627e8d5 e169c58d bc57ac4c 9b00dbd8 -philox4x32 10 ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff 408f276d 41c83b0e a20bc7c6 6d5451fd -philox4x32 10 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 d16cfe09 94fdcceb 5001e420 24126ea1 -# -philox2x64 7 0000000000000000 0000000000000000 0000000000000000 b41da69fbfefc666 511e9ce1a5534056 -philox2x64 7 ffffffffffffffff ffffffffffffffff ffffffffffffffff a4696cc04462015d 724782dae17169e9 -philox2x64 7 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 98ed1534392bf372 67528b1568882fd5 -philox2x64 10 0000000000000000 0000000000000000 0000000000000000 ca00a0459843d731 66c24222c9a845b5 -philox2x64 10 ffffffffffffffff ffffffffffffffff ffffffffffffffff 65b021d60cd8310f 4d02f3222f86df20 -philox2x64 10 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 0a5e742c2997341c b0f883d38000de5d -# -philox4x64 7 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 5dc8ee6268ec62cd 139bc570b6c125a0 84d6deb4fb65f49e aff7583376d378c2 -philox4x64 7 ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff 071dd84367903154 48e2bbdc722b37d1 6afa9890bb89f76c 9194c8d8ada56ac7 -philox4x64 7 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c 513a366704edf755 f05d9924c07044d3 bef2cb9cbea74c6c 8db948de4caa1f8a -philox4x64 10 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 16554d9eca36314c db20fe9d672d0fdc d7e772cee186176b 7e68b68aec7ba23b -philox4x64 10 ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff 87b092c3013fe90b 438c3c67be8d0224 9cc7d7c69cd777b6 a09caebf594f0ba0 -philox4x64 10 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c a528f45403e61d95 38c72dbd566e9788 a5a1610e72fd18b5 57bd43b5e52b7fe6 -# -threefry2x32 13 00000000 00000000 00000000 00000000 9d1c5ec6 8bd50731 -threefry2x32 13 ffffffff ffffffff ffffffff ffffffff fd36d048 2d17272c -threefry2x32 13 243f6a88 85a308d3 13198a2e 03707344 ba3e4725 f27d669e -threefry2x32 20 00000000 00000000 00000000 00000000 6b200159 99ba4efe -threefry2x32 20 ffffffff ffffffff ffffffff ffffffff 1cb996fc bb002be7 -threefry2x32 20 243f6a88 85a308d3 13198a2e 03707344 c4923a9c 483df7a0 -threefry2x32 32 00000000 00000000 00000000 00000000 cee3d47e a23dfd5c -threefry2x32 32 ffffffff ffffffff ffffffff ffffffff 6e2fe0d0 b1b76f82 -threefry2x32 32 243f6a88 85a308d3 13198a2e 03707344 e2827716 c3c05cdf -# -threefry4x32 13 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 531c7e4f 39491ee5 2c855a92 3d6abf9a -threefry4x32 13 ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff c4189358 1c9cc83a d5881c67 6a0a89e0 -threefry4x32 13 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 082efa98 ec4e6c89 4aa71d8f 734738c2 431fc6a8 ae6debf1 -threefry4x32 20 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 9c6ca96a e17eae66 fc10ecd4 5256a7d8 -threefry4x32 20 ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff 2a881696 57012287 f6c7446e a16a6732 -threefry4x32 20 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 082efa98 ec4e6c89 59cd1dbb b8879579 86b5d00c ac8b6d84 -threefry4x32 72 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 93171da6 9220326d b392b7b1 ff58a002 -threefry4x32 72 ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff 60743f3d 9961e684 aab21c34 8c65fb7d -threefry4x32 72 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0 082efa98 ec4e6c89 09930adf 7f27bd55 9ed68ce1 97f803f6 -# -threefry2x64 13 0000000000000000 0000000000000000 0000000000000000 0000000000000000 f167b032c3b480bd e91f9fee4b7a6fb5 -threefry2x64 13 ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ccdec5c917a874b1 4df53abca26ceb01 -threefry2x64 13 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 c3aac71561042993 3fe7ae8801aff316 -threefry2x64 20 0000000000000000 0000000000000000 0000000000000000 0000000000000000 c2b6e3a8c2c69865 6f81ed42f350084d -threefry2x64 20 ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff e02cb7c4d95d277a d06633d0893b8b68 -threefry2x64 20 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 263c7d30bb0f0af1 56be8361d3311526 -threefry2x64 32 0000000000000000 0000000000000000 0000000000000000 0000000000000000 38ba854d7f13cfb3 d02fca729d54fadc -threefry2x64 32 ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff 6b532f4f6e288646 0388f1ec135ee18e -threefry2x64 32 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 dad492f32efbd0c4 b6d7d0cd1f193e84 -# -threefry4x64 13 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 4071fabee1dc8e05 02ed3113695c9c62 397311b5b89f9d49 e21292c3258024bc -threefry4x64 13 ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff 7eaed935479722b5 90994358c429f31c 496381083e07a75b 627ed0d746821121 -threefry4x64 13 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c c0ac29b7c97c50dd 3f84d5b5b5470917 4361288ef9c1900c 8717291521782833 0d19db18c20cf47e a0b41d63ac8581e5 -threefry4x64 20 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 09218ebde6c85537 55941f5266d86105 4bd25e16282434dc ee29ec846bd2e40b -threefry4x64 20 ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff 29c24097942bba1b 0371bbfb0f6f4e11 3c231ffa33f83a1c cd29113fde32d168 -threefry4x64 20 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c be5466cf34e90c6c c0ac29b7c97c50dd a7e8fde591651bd9 baafd0c30138319b 84a5c1a729e685b9 901d406ccebc1ba4 -threefry4x64 72 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 0000000000000000 94eeea8b1f2ada84 adf103313eae6670 952419a1f4b16d53 d83f13e63c9f6b11 -threefry4x64 72 ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff ffffffffffffffff 11518c034bc1ff4c 193f10b8bcdcc9f7 d024229cb58f20d8 563ed6e48e05183f -threefry4x64 72 243f6a8885a308d3 13198a2e03707344 a4093822299f31d0 082efa98ec4e6c89 452821e638d01377 be5466cf34e90c6c be5466cf34e90c6c c0ac29b7c97c50dd acf412ccaa3b2270 c9e99bd53f2e9173 43dad469dc825948 fbb19d06c8a2b4dc \ No newline at end of file diff --git a/test/basic_rng/r123_rngNxW.mm b/test/basic_rng/r123_rngNxW.mm deleted file mode 100644 index 9112b11d..00000000 --- a/test/basic_rng/r123_rngNxW.mm +++ /dev/null @@ -1,53 +0,0 @@ -/* -Copyright 2010-2011, D. E. Shaw Research. -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 D. E. Shaw Research nor the names of its - contributors may be used to endorse or promote products derived from - this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR -A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT -OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - -/* - We need this file of some crazy logic in test_random123.cc. - It very deliberately DOES NOT use the "#pragma once" directive. - - Please just ignore any warning messages from code linters or - integrated development environments like VSCode. They won't be - able to correctly infer how we intend to use this file. -*/ - -RNGNxW_TPL(philox, 2, 32) -RNGNxW_TPL(philox, 4, 32) -RNGNxW_TPL(threefry, 2, 32) -RNGNxW_TPL(threefry, 4, 32) -#if R123_USE_64BIT -#if R123_USE_PHILOX_64BIT -RNGNxW_TPL(philox, 2, 64) -RNGNxW_TPL(philox, 4, 64) -#endif -RNGNxW_TPL(threefry, 2, 64) -RNGNxW_TPL(threefry, 4, 64) -#endif diff --git a/test/basic_rng/test_discrete.cc b/test/basic_rng/test_discrete.cc index 54c504ba..1354995b 100644 --- a/test/basic_rng/test_discrete.cc +++ b/test/basic_rng/test_discrete.cc @@ -169,21 +169,21 @@ class TestSampleIndices : public ::testing::Test static void test_updated_rngstates_iid_uniform() { RNGState seed; int offset = 3456; - seed.counter.incr(offset); + seed.advance(offset); int n = 40; int k = 17; vector unimportant(2*k); auto s1 = sample_indices_iid_uniform(n, k, unimportant.data(), seed); auto s2 = sample_indices_iid_uniform(n, k, unimportant.data(), s1); // check that counter increments are the same for the two samples of k indices. - auto total_2call = s2.counter.v[0]; - EXPECT_EQ(total_2call-offset, 2*(s1.counter.v[0]-offset)); + auto total_2call = s2.counter[0]; + EXPECT_EQ(total_2call-offset, 2*(s1.counter[0]-offset)); // check that the counter increment for a single sample of size 2k is (a) no larger // than the total increment for two samples of size k, and (b) is at most one less // than the total increment for two samples of size k. auto t = sample_indices_iid_uniform(n, 2*k, unimportant.data(), seed); - auto total_1call = t.counter.v[0]; + auto total_1call = t.counter[0]; EXPECT_LE( total_1call, total_2call ); EXPECT_LE( total_2call, total_1call + 1); } @@ -191,7 +191,7 @@ class TestSampleIndices : public ::testing::Test static void test_updated_rngstates_iid() { RNGState seed; int offset = 8675309; - seed.counter.incr(offset); + seed.advance(offset); int n = 29; int k = 13; vector unimportant(2*k); @@ -201,14 +201,14 @@ class TestSampleIndices : public ::testing::Test auto s1 = sample_indices_iid(n, cdf.data(), k, unimportant.data(), seed); auto s2 = sample_indices_iid(n, cdf.data(), k, unimportant.data(), s1); // check that counter increments are the same for the two samples of k indices. - auto total_2call = s2.counter.v[0]; - EXPECT_EQ(total_2call-offset, 2*(s1.counter.v[0]-offset)); + auto total_2call = s2.counter[0]; + EXPECT_EQ(total_2call-offset, 2*(s1.counter[0]-offset)); // check that the counter increment for a single sample of size 2k is (a) no larger // than the total increment for two samples of size k, and (b) is at most one less // than the total increment for two samples of size k. auto t = sample_indices_iid(n, cdf.data(), 2*k, unimportant.data(), seed); - auto total_1call = t.counter.v[0]; + auto total_1call = t.counter[0]; EXPECT_LE( total_1call, total_2call ); EXPECT_LE( total_2call, total_1call + 1); } @@ -304,7 +304,7 @@ class TestSampleIndices : public ::testing::Test static void test_updated_rngstates_fisher_yates() { RNGState seed; int offset = 306; - seed.counter.incr(offset); + seed.advance(offset); int n = 29; int k = 17; int r1 = 1; @@ -315,12 +315,12 @@ class TestSampleIndices : public ::testing::Test auto s1 = repeated_fisher_yates(k, n, r1, twocall.data(), seed); auto s2 = repeated_fisher_yates(k, n, r2, twocall.data() + r1*k, s1); - auto ctr_twocall = (int) s2.counter.v[0]; - auto expect_incr = (int) std::ceil(((float)r_total/r1)*(s1.counter.v[0]-offset)); + auto ctr_twocall = (int) s2.counter[0]; + auto expect_incr = (int) std::ceil(((float)r_total/r1)*(s1.counter[0]-offset)); EXPECT_EQ(ctr_twocall - offset, expect_incr); auto t = repeated_fisher_yates(k, n, r_total, onecall.data(), seed); - auto ctr_onecall = t.counter.v[0]; + auto ctr_onecall = t.counter[0]; EXPECT_EQ( ctr_onecall, ctr_twocall ); auto msg = RandBLAS::testing::buffs_approx_equal(onecall.data(), twocall.data(), r_total*k, __RANDBLAS_PRETTY_FUNCTION__, __FILE__, __LINE__); diff --git a/test/basic_rng/test_distortion.cc b/test/basic_rng/test_distortion.cc index 89a44a6b..31a18a34 100644 --- a/test/basic_rng/test_distortion.cc +++ b/test/basic_rng/test_distortion.cc @@ -53,7 +53,7 @@ class TestSubspaceDistortion : public ::testing::Test { DenseDist D(d, N, name); std::vector S(d*N); std::cout << "(d, N) = ( " << d << ", " << N << " )\n"; - RandBLAS::RNGState state(key); + RandBLAS::DefaultRNGState state(key); auto next_state = RandBLAS::fill_dense(D, S.data(), state); T inv_stddev = (name == ScalarDist::Gaussian) ? (T) 1.0 : (T) 1.0; blas::scal(d*N, inv_stddev / std::sqrt(d), S.data(), 1); diff --git a/test/basic_rng/test_distributions.cc b/test/basic_rng/test_distributions.cc new file mode 100644 index 00000000..38e7f4e9 --- /dev/null +++ b/test/basic_rng/test_distributions.cc @@ -0,0 +1,202 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +template +struct FixedState { + using res_t = std::array; + + res_t values{}; + std::uint64_t blocks{}; + + constexpr void generate(res_t& output) const noexcept { + output = values; + } + + constexpr void advance(std::uint64_t amount) noexcept { + blocks += amount; + } +}; + +static_assert(RandBLAS::GeneratorState>); + +template +concept CanGenerateNormals = requires(State const& state) { + RandBLAS::rng::boxmul::generate(state); +}; + +template +void expect_near_reference(Real actual, Real expected) { + auto scale = std::max(Real{1}, std::abs(expected)); + auto tolerance = Real{8} * std::numeric_limits::epsilon() * scale; + EXPECT_NEAR(actual, expected, tolerance); +} + +TEST(DistributionConversion, U01MatchesRetainedReferenceValues) { + using RandBLAS::rng::u01; + + EXPECT_EQ(u01(UINT32_C(0)), 0x1p-33f); + EXPECT_EQ(u01(UINT32_C(1)), 0x1.8p-32f); + EXPECT_EQ(u01(UINT32_C(0x80000000)), 0x1p-1f); + EXPECT_EQ(u01(UINT32_MAX), 0x1p+0f); + EXPECT_EQ(u01(UINT32_C(0x243f6a88)), 0x1.21fb54p-3f); + + EXPECT_EQ(u01(UINT32_C(0)), 0x1p-33); + EXPECT_EQ(u01(UINT32_C(1)), 0x1.8p-32); + EXPECT_EQ(u01(UINT32_C(0x80000000)), 0x1.00000001p-1); + EXPECT_EQ(u01(UINT32_MAX), 0x1.ffffffffp-1); + EXPECT_EQ(u01(UINT32_C(0x243f6a88)), 0x1.21fb5444p-3); + + EXPECT_EQ(u01(UINT64_C(0)), 0x1p-65); + EXPECT_EQ(u01(UINT64_C(1)), 0x1.8p-64); + EXPECT_EQ(u01(UINT64_C(0x8000000000000000)), 0x1p-1); + EXPECT_EQ(u01(UINT64_MAX), 0x1p+0); + EXPECT_EQ(u01(UINT64_C(0x243f6a8885a308d3)), + 0x1.21fb54442d184p-3); +} + +TEST(DistributionConversion, U01HasTheRetainedEndpointConvention) { + EXPECT_GT(RandBLAS::rng::u01(UINT32_C(0)), 0.0f); + EXPECT_LE(RandBLAS::rng::u01(UINT32_MAX), 1.0f); + EXPECT_EQ(RandBLAS::rng::u01(UINT32_MAX), 1.0f); + + EXPECT_GT(RandBLAS::rng::u01(UINT32_C(0)), 0.0); + EXPECT_LT(RandBLAS::rng::u01(UINT32_MAX), 1.0); + EXPECT_GT(RandBLAS::rng::u01(UINT64_C(0)), 0.0); + EXPECT_EQ(RandBLAS::rng::u01(UINT64_MAX), 1.0); +} + +TEST(DistributionConversion, Uneg11MatchesRetainedReferenceValues) { + using Policy = RandBLAS::rng::uneg11; + + EXPECT_EQ(Policy::convert(UINT32_C(0)), 0x1p-32f); + EXPECT_EQ(Policy::convert(UINT32_C(1)), 0x1.8p-31f); + EXPECT_EQ(Policy::convert(UINT32_C(0x80000000)), -0x1p+0f); + EXPECT_EQ(Policy::convert(UINT32_MAX), -0x1p-32f); + EXPECT_EQ(Policy::convert(UINT32_C(0x243f6a88)), + 0x1.21fb54p-2f); + + EXPECT_EQ(Policy::convert(UINT32_C(0)), 0x1p-32); + EXPECT_EQ(Policy::convert(UINT32_C(1)), 0x1.8p-31); + EXPECT_EQ(Policy::convert(UINT32_C(0x80000000)), + -0x1.fffffffep-1); + EXPECT_EQ(Policy::convert(UINT32_MAX), -0x1p-32); + EXPECT_EQ(Policy::convert(UINT32_C(0x243f6a88)), + 0x1.21fb5444p-2); + + EXPECT_EQ(Policy::convert(UINT64_C(0)), 0x1p-64); + EXPECT_EQ(Policy::convert(UINT64_C(1)), 0x1.8p-63); + EXPECT_EQ(Policy::convert(UINT64_C(0x8000000000000000)), + -0x1p+0); + EXPECT_EQ(Policy::convert(UINT64_MAX), -0x1p-64); + EXPECT_EQ(Policy::convert(UINT64_C(0x243f6a8885a308d3)), + 0x1.21fb54442d184p-2); +} + +TEST(DistributionConversion, Uneg11IsClosedAndNeverZero) { + std::array inputs{ + 0, 1, UINT32_C(0x7fffffff), UINT32_C(0x80000000), UINT32_MAX}; + for (auto input : inputs) { + auto value = RandBLAS::rng::uneg11::convert(input); + EXPECT_GE(value, -1.0f); + EXPECT_LE(value, 1.0f); + EXPECT_NE(value, 0.0f); + } +} + +TEST(DistributionConversion, BoxMullerMatchesRetainedReferences) { + auto result32 = RandBLAS::rng::boxmuller( + UINT32_C(0x243f6a88), UINT32_C(0x85a308d3)); + expect_near_reference(result32[0], 0x1.c5857cp-1f); + expect_near_reference(result32[1], 0x1.6f9a8ep-1f); + + auto result64 = RandBLAS::rng::boxmuller( + UINT64_C(0x243f6a8885a308d3), + UINT64_C(0x13198a2e03707344)); + expect_near_reference(result64[0], 0x1.c51c6804651e6p+0); + expect_near_reference(result64[1], 0x1.6f4563170165cp+0); +} + +TEST(DistributionConversion, BoxMullerUsesAngleThenRadiusAndReturnsSinThenCos) { + constexpr std::uint32_t angle_word = UINT32_C(0x243f6a88); + constexpr std::uint32_t radius_word = UINT32_C(0x85a308d3); + constexpr float pi = 3.1415926535897932f; + auto angle = pi * RandBLAS::rng::uneg11::convert(angle_word); + auto radius = std::sqrt(-2.0f * + std::log(RandBLAS::rng::u01(radius_word))); + auto result = RandBLAS::rng::boxmuller(angle_word, radius_word); + + expect_near_reference(result[0], std::sin(angle) * radius); + expect_near_reference(result[1], std::cos(angle) * radius); +} + +TEST(DistributionPolicy, GeneratesOneFixedBlockWithoutAdvancingState) { + FixedState state{{ + UINT32_C(0x243f6a88), UINT32_C(0x85a308d3), + UINT32_C(0x13198a2e), UINT32_C(0x03707344)}}; + typename decltype(state)::res_t bits{}; + state.generate(bits); + auto state_before = state.values; + auto uniform = RandBLAS::rng::uneg11::generate(state); + auto normal = RandBLAS::rng::boxmul::generate(state); + + EXPECT_EQ(state.values, state_before); + EXPECT_EQ(uniform.size(), bits.size()); + EXPECT_EQ(normal.size(), bits.size()); + for (std::size_t i = 0; i < bits.size(); ++i) { + EXPECT_EQ(uniform[i], + RandBLAS::rng::uneg11::convert(bits[i])); + } + for (std::size_t i = 0; i < bits.size(); i += 2) { + auto pair = RandBLAS::rng::boxmuller(bits[i], bits[i + 1]); + EXPECT_EQ(normal[i], pair[0]); + EXPECT_EQ(normal[i + 1], pair[1]); + } +} + +TEST(DistributionPolicy, RejectsOddNormalBlockLengths) { + static_assert(CanGenerateNormals>); + static_assert(!CanGenerateNormals>); + SUCCEED(); +} + +} // namespace diff --git a/test/basic_rng/test_philox.cc b/test/basic_rng/test_philox.cc new file mode 100644 index 00000000..22e828db --- /dev/null +++ b/test/basic_rng/test_philox.cc @@ -0,0 +1,209 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +struct KatRecord { + std::string family; + std::size_t rounds{}; + std::vector words; + std::size_t line{}; +}; + +std::vector read_vectors() { + std::ifstream input(PHILOX_KAT_VECTORS_PATH); + EXPECT_TRUE(input.is_open()) << PHILOX_KAT_VECTORS_PATH; + + std::vector records; + std::string text; + for (std::size_t line = 1; std::getline(input, text); ++line) { + if (text.empty() || text.front() == '#') { + continue; + } + + KatRecord record; + record.line = line; + std::istringstream fields(text); + fields >> record.family >> record.rounds; + std::uint64_t word; + while (fields >> std::hex >> word) { + record.words.push_back(word); + } + EXPECT_TRUE(fields.eof()) << "malformed vector at line " << line; + records.push_back(std::move(record)); + } + return records; +} + +template +void check_vector(KatRecord const& record) { + using word_t = typename Engine::word_t; + typename Engine::ctr_t counter{}; + typename Engine::key_t key{}; + typename Engine::res_t expected{}; + + constexpr std::size_t expected_words = + Engine::ctr_t::static_size + Engine::key_t::static_size + + std::tuple_size_v; + ASSERT_EQ(record.words.size(), expected_words) << "line " << record.line; + + std::size_t offset = 0; + for (auto& word : counter.words) { + word = static_cast(record.words[offset++]); + } + for (auto& word : key.words) { + word = static_cast(record.words[offset++]); + } + for (auto& word : expected) { + word = static_cast(record.words[offset++]); + } + + typename Engine::res_t actual; + actual.fill(std::numeric_limits::max()); + auto counter_before = counter; + auto key_before = key; + + Engine{}.generate(counter, key, actual); + + EXPECT_EQ(actual, expected) << "line " << record.line; + EXPECT_EQ(counter, counter_before) << "line " << record.line; + EXPECT_EQ(key, key_before) << "line " << record.line; +} + +template +void check_round(KatRecord const& record, std::index_sequence) { + bool matched = false; + ([&] { + if (record.rounds == Rounds) { + check_vector>(record); + matched = true; + } + }(), + ...); + EXPECT_TRUE(matched) << "unsupported round count at line " << record.line; +} + +template +consteval bool has_expected_shape() { + using word_t = typename Engine::word_t; + return std::unsigned_integral && + std::numeric_limits::digits == W && + Engine::ctr_t::static_size == N && + Engine::key_t::static_size == N / 2 && + std::tuple_size_v == N && + std::same_as; +} + +TEST(Philox, AllKnownAnswerVectors) { + auto records = read_vectors(); + ASSERT_EQ(records.size(), 4U * 17U * 3U); + + std::array counts{}; + for (auto const& record : records) { + ASSERT_LE(record.rounds, 16U) << "line " << record.line; + std::size_t family_index; + if (record.family == "philox2x32") { + family_index = 0; + check_round<2, 32>(record, std::make_index_sequence<17>{}); + } else if (record.family == "philox4x32") { + family_index = 1; + check_round<4, 32>(record, std::make_index_sequence<17>{}); + } else if (record.family == "philox2x64") { + family_index = 2; + check_round<2, 64>(record, std::make_index_sequence<17>{}); + } else if (record.family == "philox4x64") { + family_index = 3; + check_round<4, 64>(record, std::make_index_sequence<17>{}); + } else { + FAIL() << "unknown family at line " << record.line; + continue; + } + ++counts[family_index * 17 + record.rounds]; + } + + for (auto count : counts) { + EXPECT_EQ(count, 3U); + } +} + +TEST(Philox, RoundZeroCopiesCounterAndOverwritesOutput) { + using Engine = RandBLAS::rng::Philox<4, 64, 0>; + Engine::ctr_t counter{{0, 1, UINT64_C(0x8000000000000000), + UINT64_C(0x0123456789abcdef)}}; + Engine::key_t key{{UINT64_MAX, UINT64_C(0x3141592653589793)}}; + Engine::res_t output; + output.fill(UINT64_MAX); + + Engine{}.generate(counter, key, output); + + EXPECT_EQ(output, (Engine::res_t{counter[0], counter[1], counter[2], + counter[3]})); +} + +TEST(Philox, PublicTypesHaveExpectedShapes) { + static_assert(has_expected_shape, 2, 32>()); + static_assert(has_expected_shape, 4, 32>()); + static_assert(has_expected_shape, 2, 64>()); + static_assert(has_expected_shape, 4, 64>()); + SUCCEED(); +} + +TEST(Philox, MakeKeyUsesLittleEndianSeedAddition) { + constexpr std::uint64_t seed = UINT64_C(0x0123456789abcdef); + + EXPECT_EQ((RandBLAS::rng::Philox<2, 32, 10>::make_key(0)), + (RandBLAS::rng::Philox<2, 32, 10>::key_t{})); + EXPECT_EQ((RandBLAS::rng::Philox<2, 32, 10>::make_key(seed)), + (RandBLAS::rng::Philox<2, 32, 10>::key_t{{0x89abcdef}})); + EXPECT_EQ((RandBLAS::rng::Philox<4, 32, 10>::make_key(seed)), + (RandBLAS::rng::Philox<4, 32, 10>::key_t{{0x89abcdef, + 0x01234567}})); + EXPECT_EQ((RandBLAS::rng::Philox<2, 64, 10>::make_key(seed)), + (RandBLAS::rng::Philox<2, 64, 10>::key_t{{seed}})); + EXPECT_EQ((RandBLAS::rng::Philox<4, 64, 10>::make_key(seed)), + (RandBLAS::rng::Philox<4, 64, 10>::key_t{{seed, 0}})); +} + +} // namespace diff --git a/test/basic_rng/test_r123.cc b/test/basic_rng/test_r123.cc deleted file mode 100644 index 90922d1c..00000000 --- a/test/basic_rng/test_r123.cc +++ /dev/null @@ -1,810 +0,0 @@ -/* -Copyright 2010-2011, D. E. Shaw Research. -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 D. E. Shaw Research nor the names of its - contributors may be used to endorse or promote products derived from - this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR -A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT -OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - -#include - -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include - -#include -#include -#include -#include -#include -#include - - -#define LINESIZE 1024 -#ifdef _MSC_FULL_VER -#define strtoull _strtoui64 -// ^ Needed to define the strtou32 and strtou64 functions. -#pragma warning (disable : 4521) -// ^ Engines have multiple copy constructors, quite legal C++, disable MSVC complaint -#endif - -int verbose = 0; -int debug = 0; - -// MARK: I/O and conversions - -/* strdup may or may not be in string.h, depending on the value - of the pp-symbol _XOPEN_SOURCE and other arcana. Just - do it ourselves. - Mnemonic: "ntcs" = "nul-terminated character string" */ -char *ntcsdup(const char *s){ - char *p = (char *)malloc(strlen(s)+1); - strcpy(p, s); - return p; -} - -// Functions to read a (portion of) a string in a given base and convert it to -// an unsigned integer. -// -// These functions differ from std::from_chars in how they handle white space. -// Specifically, they strip leading whitespace, and then they stop reading as -// soon as they reach a non-numeric character. (Note that the "a" in 257a3673 -// counts as a numeric character if we're reading in hexadecimal format.) -uint32_t strtou32(const char *p, char **endp, int base){ - uint32_t ret; - errno = 0; - ret = strtoul(p, endp, base); - assert(errno==0); - return ret; -} -uint64_t strtou64(const char *p, char **endp, int base){ - uint64_t ret; - errno = 0; - ret = strtoull(p, endp, base); - assert(errno==0); - return ret; -} - -// A helper function to print unsigned integers in hexadecimal format, with leading zeros if necessary. -template -void prtu(std::ostream& os, T val) { - os << std::hex << std::setw(std::numeric_limits::digits / 4) << std::setfill('0') << val; - assert(!os.bad()); -} -void prtu32(std::ostream& os, uint32_t v) { prtu(os, v); } -void prtu64(std::ostream& os, uint64_t v) { prtu(os, v); } - -#define PRINTARRAY(ARR, fp) \ -do { \ - char ofmt[64]; \ - size_t xj; \ - /* use %lu and the cast (instead of z) for portability to Microsoft, sizeof(v[0]) should fit easily in an unsigned long. Avoid inttypes for the same reason. */ \ - snprintf(ofmt, sizeof(ofmt), " %%0%lullx", (unsigned long)sizeof(ARR.v[0])*2UL); \ - for (xj = 0; xj < sizeof(ARR.v)/sizeof(ARR.v[0]); xj++) { \ - fprintf(fp, ofmt, (unsigned long long) ARR.v[xj]); \ - } \ -} while(0) - -#define PRINTLINE(NAME, N, W, R, ictr, ukey, octr, fp) \ -do { \ - fprintf(fp, "%s %d ", #NAME #N "x" #W, R); \ - PRINTARRAY(ictr, fp); \ - putc(' ', fp); \ - PRINTARRAY(ukey, fp); \ - putc(' ', fp); \ - PRINTARRAY(octr, fp); \ - putc('\n', fp); \ - fflush(fp); \ -} while(0) - -// MARK: Base generator test -// -// There's a lot of code involved in this test. The code can roughly -// be broken down into three categories. -// -// Category 1: code generated from compiler directives -// -// This pattern is left over from our adaptation of Random123 tests, -// which have to compile whether interpreted as C or C++ source. It -// uses compiler directives to accomplish what something roughly -// equivalent to C++ templating and metaprogramming. -// -// The code specifically generates the following identifiers. -// -// method_e::NxW_e (enum members) -// NxW_kat (structs) -// kat_instance.NxW_data (members of type NxW_kat) -// read_NxW (functions) -// report_NxWerror (functions) -// -// Category 2: helper functions -// -// The base_rng_test_[arrange,act,assert] functions are slight adaptations -// of functions that appeared in Random123 testing infrastructure. Their -// names indicte their roles in the common "arrange, act, assert" pattern of -// writing unit tests. Their precise descriptions are complicated. -// -// Category 3: the main runner -// -// This manages all calls to the helper functions defined in Category 2. -// - -enum method_e{ -#define RNGNxW_TPL(base, N, W) base##N##x##W##_e, -#include "r123_rngNxW.mm" -#undef RNGNxW_TPL - last -}; - -#define RNGNxW_TPL(base, N, W) \ - struct base##N##x##W##_kat { \ - base##N##x##W##_ctr_t ctr; \ - base##N##x##W##_ukey_t ukey; \ - base##N##x##W##_ctr_t expected; \ - base##N##x##W##_ctr_t computed; \ - }; -#include "r123_rngNxW.mm" -#undef RNGNxW_TPL - -struct kat_instance { - enum method_e method; - unsigned nrounds; - union{ -#define RNGNxW_TPL(base, N, W) base##N##x##W##_kat base##N##x##W##_data; -#include "r123_rngNxW.mm" -#undef RNGNxW_TPL - // Sigh... For those platforms that lack uint64_t, carve - // out 128 bytes for the counter, key, expected, and computed. - char justbytes[128]; - }u; -}; - -#define RNGNxW_TPL(base, N, W) \ -int read_##base##N##x##W(const char *line, kat_instance* tinst){ \ - size_t i; \ - int nchar; \ - const char *p = line; \ - char *newp; \ - size_t nkey = sizeof(tinst->u.base##N##x##W##_data.ukey.v)/sizeof(tinst->u.base##N##x##W##_data.ukey.v[0]); \ - tinst->method = base##N##x##W##_e; \ - sscanf(p, "%u%n", &tinst->nrounds, &nchar); \ - p += nchar; \ - for(i=0; iu.base##N##x##W##_data.ctr.v[i] = strtou##W(p, &newp, 16); \ - p = newp; \ - } \ - for(i=0; iu.base##N##x##W##_data.ukey.v[i] = strtou##W(p, &newp, 16); \ - p = newp; \ - } \ - for(i=0; iu.base##N##x##W##_data.expected.v[i] = strtou##W(p, &newp, 16); \ - p = newp; \ - } \ - /* set the computed to 0xca. If the test fails to set computed, we'll see cacacaca in the FAILURE notices */ \ - memset(tinst->u.base##N##x##W##_data.computed.v, 0xca, sizeof(tinst->u.base##N##x##W##_data.computed.v)); \ - return 1; \ -} -#include "r123_rngNxW.mm" -#undef RNGNxW_TPL - -#define RNGNxW_TPL(base, N, W) \ -void report_##base##N##x##W##error(int &nfailed, const kat_instance *ti){ \ - size_t i; \ - size_t nkey = sizeof(ti->u.base##N##x##W##_data.ukey.v)/sizeof(ti->u.base##N##x##W##_data.ukey.v[0]); \ - std::stringstream ss; \ - ss << "FAIL: expected: "; \ - ss << #base #N "x" #W " " << ti->nrounds; \ - for(i=0; iu.base##N##x##W##_data.ctr.v[i]); \ - } \ - for(i=0; iu.base##N##x##W##_data.ukey.v[i]); \ - } \ - for(i=0; iu.base##N##x##W##_data.expected.v[i]); \ - } \ - ss << "\n"; \ - \ - ss << "FAIL: computed: "; \ - ss << #base #N "x" #W " " << ti->nrounds; \ - for(i=0; iu.base##N##x##W##_data.ctr.v[i]); \ - } \ - for(i=0; iu.base##N##x##W##_data.ukey.v[i]); \ - } \ - for(i=0; iu.base##N##x##W##_data.computed.v[i]); \ - } \ - ss << "\n"; \ - FAIL() << ss.str(); \ - nfailed++; \ -} -#include "r123_rngNxW.mm" -#undef RNGNxW_TPL - -struct UnknownKatTracker { - const static int MAXUNKNOWNS = 20; - int num_unknowns = 0; - const char *unknown_names[MAXUNKNOWNS]; - int unknown_counts[MAXUNKNOWNS]; -}; - -void register_unknown(UnknownKatTracker &ukt, const char *name){ - int i; - for(i=0; i< ukt.num_unknowns; ++i){ - if( strcmp(name, ukt.unknown_names[i]) == 0 ){ - ukt.unknown_counts[i]++; - return; - } - } - if( i >= ukt.MAXUNKNOWNS ){ - FAIL() << "Too many unknown rng types. Bye.\n"; - } - ukt.num_unknowns++; - ukt.unknown_names[i] = ntcsdup(name); - ukt.unknown_counts[i] = 1; -} - -void base_rng_test_arrange(const char *line, kat_instance* tinst, UnknownKatTracker &ukt, bool &flag){ - int nchar; - char name[LINESIZE]; - if( line[0] == '#') { - flag = false; - return; - } - sscanf(line, "%s%n", name, &nchar); - /* skip any tests that require AESNI */ - if(strncmp(name, "aes", 3)==0 || strncmp(name, "ars", 3)==0){ - register_unknown(ukt, name); - flag = false; - return; - } -#define RNGNxW_TPL(base, N, W) \ - if(strcmp(name, #base #N "x" #W) == 0) { \ - flag = (bool) read_##base##N##x##W(line+nchar, tinst); \ - return; \ - } -#include "r123_rngNxW.mm" -#undef RNGNxW_TPL - - register_unknown(ukt, name); - flag = false; - return; -} - -static int murng_reported; -static int engine_reported; - -template -void base_rng_test_act(kat_instance* ti){ - GEN g; - struct gdata{ - typename GEN::ctr_type ctr; - typename GEN::ukey_type ukey; - typename GEN::ctr_type expected; - typename GEN::ctr_type computed; - }; - gdata data; - // use memcpy. A reinterpret_cast would violate strict aliasing. - std::memcpy(&data, &ti->u, sizeof(data)); - data.computed = g(data.ctr, data.ukey); - - // Before we return, let's make sure that MicroURNG and - // Engine work as expeccted. This doesn't really "fit" the - // execution model of kat.c, which just expects us to fill in - // ti->u.computed, so we report the error by failing to write back - // the computed data item in the (hopefully unlikely) event that - // things don't match up as expected. - int errs = 0; - - // MicroURNG: throws if the top 32 bits of the high word of ctr - // are non-zero. - typedef typename GEN::ctr_type::value_type value_type; - - value_type hibits = data.ctr[data.ctr.size()-1]>>( std::numeric_limits::digits - 32 ); - try{ - r123::MicroURNG urng(data.ctr, data.ukey); - if(hibits) - errs++; // Should have thrown. - for (size_t i = 0; i < data.expected.size(); i++) { - size_t j = data.expected.size() - i - 1; - if (data.expected[j] != urng()) { - errs++; - } - } - }catch(std::runtime_error& /*ignored*/){ - // A runtime_error is expected from the constructor - // when hibit is set. - if(!hibits) - errs++; - } - if(errs && (murng_reported++ == 0)) - std::cerr << "Error in MicroURNG, will appear as \"computed\" value of zero in error summary\n"; - - // Engine - // N.B. exercising discard() arguably belongs in ut_Engine.cpp - typedef r123::Engine Etype; - typedef typename GEN::ctr_type::value_type value_type; - Etype e(data.ukey); - typename GEN::ctr_type c = data.ctr; - value_type c0; - if( c[0] > 0 ){ - c0 = c[0]-1; - }else{ - // N.B. Assume that if c[0] is 0, then so are all the - // others. Arrange to "roll over" to {0,..,0} on the first - // counter-increment. Alternatively, we could just - // skip the test for this case... - c.fill(std::numeric_limits::max()); - c0 = c[0]; - } - c[0] /= 3; - e.setcounter(c, 0); - if( c0 > c[0] ){ - // skip one value by calling e() - (void)e(); - if (c0 > c[0]+1) { - // skip many values by calling discard() - R123_ULONG_LONG ndiscard = (c0 - c[0] - 1); - // Take care not to overflow the long long - if( ndiscard >= std::numeric_limits::max() / c.size() ){ - for(size_t j=0; j, will appear as \"computed\" value of zero in error summary\n"; - } - } - - // Signal an error to the caller by *not* copying back - // the computed data object into the ti - if(errs == 0) - std::memcpy(&ti->u, &data, sizeof(data)); -} - -void base_rng_test_assert(int &nfailed, const kat_instance *tests, unsigned ntests){ - unsigned i; - char zeros[512] = {0}; - for(i=0; iu.base##N##x##W##_data.expected.v, N*W/8)==0){ \ - FAIL() << "kat expected all zeros? Something is wrong with the test harness!\n"; \ - nfailed++; \ - } \ - if (memcmp(ti->u.base##N##x##W##_data.computed.v, ti->u.base##N##x##W##_data.expected.v, N*W/8)) \ - report_##base##N##x##W##error(nfailed, ti); \ - break; -#include "r123_rngNxW.mm" -#undef RNGNxW_TPL - case last: ; - } - } -} - -void run_all_base_rng_kats() { - kat_instance *tests; - unsigned t, ntests = 1000; - char linebuf[LINESIZE]; - FILE *inpfile; - const char *p; - const char *inname; - int nfailed = 0; - - UnknownKatTracker ukt{}; - - inname = KAT_VECTORS_PATH; - inpfile = fopen(inname, "r"); - if (inpfile == NULL) - FAIL() << "Error opening input file " << inname << " for reading. Received error code " << errno << "\n"; - - if ((p = getenv("KATC_VERBOSE")) != NULL) - verbose = atoi(p); - - if ((p = getenv("KATC_DEBUG")) != NULL) - debug = atoi(p); - - tests = (kat_instance *) malloc(sizeof(tests[0])*ntests); - if (tests == NULL) { - FAIL() << "Could not allocate " << (unsigned long) ntests << " bytes for tests\n"; - } - t = 0; - while (fgets(linebuf, sizeof linebuf, inpfile) != NULL) { - if( t == ntests ) { - ntests *= 2; - tests = (kat_instance *)realloc(tests, sizeof(tests[0])*ntests); - if (tests == NULL) { - FAIL() << "Could not grow tests to " << (unsigned long) ntests << " bytes.\n"; - } - } - bool flag = false; - base_rng_test_arrange(linebuf, &tests[t], ukt, flag); - if( flag ) - ++t; - } - if(t==ntests){ - FAIL() << "No more space for tests? Recompile with a larger ntests\n"; - } - tests[t].method = last; - - for(int i=0; i< ukt.num_unknowns; ++i){ - printf("%d test vectors of type %s skipped\n", ukt.unknown_counts[i], ukt.unknown_names[i]); - } - printf("Perform %lu tests.\n", (unsigned long)t); - using std::map; - using std::pair; - using std::make_pair; - typedef map, void (*)(kat_instance *)> genmap_t; - genmap_t genmap; - // In C++1x, this could be staticly declared with an initializer list. - genmap[make_pair(threefry2x32_e, 13u)] = base_rng_test_act >; - genmap[make_pair(threefry2x32_e, 20u)] = base_rng_test_act >; - genmap[make_pair(threefry2x32_e, 32u)] = base_rng_test_act >; - - genmap[make_pair(threefry4x32_e, 13u)] = base_rng_test_act >; - genmap[make_pair(threefry4x32_e, 20u)] = base_rng_test_act >; - genmap[make_pair(threefry4x32_e, 72u)] = base_rng_test_act >; - - #if R123_USE_64BIT - genmap[make_pair(threefry2x64_e, 13u)] = base_rng_test_act >; - genmap[make_pair(threefry2x64_e, 20u)] = base_rng_test_act >; - genmap[make_pair(threefry2x64_e, 32u)] = base_rng_test_act >; - - genmap[make_pair(threefry4x64_e, 13u)] = base_rng_test_act >; - genmap[make_pair(threefry4x64_e, 20u)] = base_rng_test_act >; - genmap[make_pair(threefry4x64_e, 72u)] = base_rng_test_act >; - #endif - - genmap[make_pair(philox2x32_e, 7u)] = base_rng_test_act >; - genmap[make_pair(philox2x32_e, 10u)] = base_rng_test_act >; - genmap[make_pair(philox4x32_e, 7u)] = base_rng_test_act >; - genmap[make_pair(philox4x32_e, 10u)] = base_rng_test_act >; - - #if R123_USE_PHILOX_64BIT - genmap[make_pair(philox2x64_e, 7u)] = base_rng_test_act >; - genmap[make_pair(philox2x64_e, 10u)] = base_rng_test_act >; - genmap[make_pair(philox4x64_e, 7u)] = base_rng_test_act >; - genmap[make_pair(philox4x64_e, 10u)] = base_rng_test_act >; - #endif - - unsigned i; - for(i=0; imethod, ti->nrounds)); - if(p == genmap.end()) - throw std::runtime_error("pair not in map. You probably need to add more genmap entries."); - - p->second(ti); - // ^ That prepares the test data - } - - base_rng_test_assert(nfailed, tests, t); - free(tests); - if(nfailed != 0) - FAIL() << "Failed " << nfailed << " out of " << t << std::endl; - return; -} - -// MARK: histogram test - -using namespace r123; - -template -typename r123::make_unsigned::type U(T x){ return x; } - -template -typename r123::make_signed::type S(T x){ return x; } - -#define Chk(u, Rng, Ftype, _nfail_, _refhist_) do{ \ - chk(#u, #Rng, #Ftype, &u, _nfail_, _refhist_); \ - }while(0) - -template -void chk(const std::string& fname, const std::string& rngname, const std::string& ftypename, Utype f, int &nfail, std::map &refmap){ - std::string key = fname + " " + rngname + " " + ftypename; - RNG rng; - typedef typename RNG::ukey_type ukey_type; - typedef typename RNG::ctr_type ctr_type; - typedef typename RNG::key_type key_type; - - ctr_type c = {{}}; - ukey_type uk = {{}}; - key_type k = uk; - // 26 bins - 13 greater than 0 and 13 less. Why 13? Because a - // prime number seems less likely to tickle the rounding-related - // corner cases, which is aruably both good and bad. - const int NBINS=26; - - int hist[NBINS] = {}; - for(int i=0; i<1000; ++i){ - c = c.incr(); - ctr_type r = rng(c, k); - for(int j=0; j= -1.); - R123_ASSERT( u <= 1.); - int idx = (int) ((u + Ftype(1.))*Ftype(NBINS/2)); - hist[idx]++; - } - } - std::ostringstream oss; - for(int i=0; i s; - ASSERT_EQ(s.key[0], 0); - ASSERT_EQ(s.key[1], 0); - for (int i = 0; i < 4; ++i) { - ASSERT_EQ(s.counter[i], 0) << "Failed at index " << i; - } - // unsigned-int constructor - RandBLAS::RNGState t(42); - ASSERT_EQ(t.key[0], 42); - ASSERT_EQ(t.key[1], 0); - for (int i = 0; i < 4; ++i) { - ASSERT_EQ(t.counter[i], 0) << "Failed at index " << i; - } - return; - } -}; - -TEST_F(TestRNGState, uint_key_constructors) { - test_uint_key_constructors(); -} - - -class TestRandom123 : public ::testing::Test { - - protected: - - static void test_incr() { - using RNG = r123::Philox4x32; - RandBLAS::RNGState s(0); - // The "counter" array of s is a 4*32=128 bit unsigned integer. - // - // Each block is interpreted in the usual way (i.e., no need to consider differences - // between big-endian and little-endian representations). - // - // Looking across blocks, we read as as a little-endian number in base IMAX = 2^32 - 1. - // That is, if we initialize s.counter = {0,0,0,0} and then call s.counter.incr(IMAX), - // we should have s.counter = {IMAX, 0, 0, 0}, and if we make another call - // s.counter.incr(9), then we should see s.counter = {8, 1, 0, 0}. Put another way, - // if c = s.counter, then we have - // - // (128-bit integer) c == c[0] + 2^{32}*c[1] + 2^{64}*c[2] + 2^{96}*c[3] (mod 2^128 - 1) - // - // where 0 <= c[i] <= IMAX - // - uint64_t i32max = std::numeric_limits::max(); - auto c = s.counter; - ASSERT_EQ(c[0], 0); - ASSERT_EQ(c[1], 0); - ASSERT_EQ(c[2], 0); - ASSERT_EQ(c[3], 0); - - c.incr(i32max); - ASSERT_EQ(c[0], i32max); - ASSERT_EQ(c[1], 0); - ASSERT_EQ(c[2], 0); - ASSERT_EQ(c[3], 0); - - c.incr(1); - ASSERT_EQ(c[0], 0); - ASSERT_EQ(c[1], 1); - ASSERT_EQ(c[2], 0); - ASSERT_EQ(c[3], 0); - - c.incr(3); - ASSERT_EQ(c[0], 3); - ASSERT_EQ(c[1], 1); - ASSERT_EQ(c[2], 0); - ASSERT_EQ(c[3], 0); - - uint64_t two32 = ((uint64_t) 1) << 32; - - c = {0,0,0,0}; - c.incr(two32-1); - ASSERT_EQ(c[0], i32max); - ASSERT_EQ(c[1], 0); - ASSERT_EQ(c[2], 0); - ASSERT_EQ(c[3], 0); - - c = {0,0,0,0}; - c.incr(two32); - ASSERT_EQ(c[0], 0); - ASSERT_EQ(c[1], 1); - ASSERT_EQ(c[2], 0); - ASSERT_EQ(c[3], 0); - - // Let's construct 2^32 * (2^32 - 1), which is equal to (ctr_type) {0, (uint32_t) i32max, 0, 0}. - // - // Do this using the identity - // 2^32 * (2^32 - 1) == 2^64 - 2^32 - // == 2^63 + 2^63 - 2^32. - // - // Then construct 2^64, using 2^64 = (2^63) + (2^63 - 2^32) + (2^32) - uint64_t two63 = ((uint64_t) 1) << 63; - c = {0,0,0,0}; - c.incr(two63); - c.incr(two63 - two32); - ASSERT_EQ(c[0], 0); - ASSERT_EQ(c[1], i32max); - ASSERT_EQ(c[2], 0); - ASSERT_EQ(c[3], 0); - c.incr(two32); - ASSERT_EQ(c[0], 0); - ASSERT_EQ(c[1], 0); - ASSERT_EQ(c[2], 1); - ASSERT_EQ(c[3], 0); - - c = {(uint32_t) i32max, (uint32_t) i32max, (uint32_t) i32max, 0}; - c.incr(1); - ASSERT_EQ(c[0], 0); - ASSERT_EQ(c[1], 0); - ASSERT_EQ(c[2], 0); - ASSERT_EQ(c[3], 1); - return; - } -}; - -TEST_F(TestRandom123, base_generators) { - run_all_base_rng_kats(); -} - -TEST_F(TestRandom123, uniform_histograms) { - run_ut_uniform(); -} - -TEST_F(TestRandom123, big_incr) { - test_incr(); -} \ No newline at end of file diff --git a/test/basic_rng/test_repacked_output.cc b/test/basic_rng/test_repacked_output.cc new file mode 100644 index 00000000..b841cf4d --- /dev/null +++ b/test/basic_rng/test_repacked_output.cc @@ -0,0 +1,188 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace { + +struct FixedEngine { + using ctr_t = RandBLAS::rng::WordArray; + using key_t = RandBLAS::rng::WordArray; + using res_t = std::array; + + static constexpr key_t make_key(std::uint64_t seed) noexcept { + key_t key{}; + key.advance(seed); + return key; + } + + constexpr void generate(ctr_t const&, key_t const&, + res_t& output) const noexcept { + output = {0xaabbccddu, 0x01234567u}; + } +}; + +struct EngineWithoutMakeKey { + using ctr_t = FixedEngine::ctr_t; + using key_t = FixedEngine::key_t; + using res_t = FixedEngine::res_t; + + constexpr void generate(ctr_t const&, key_t const&, + res_t& output) const noexcept { + output = {0xaabbccddu, 0x01234567u}; + } +}; + +template +concept CanRepack = requires { + typename RandBLAS::rng::RepackedOutput; +}; + +template +concept HasMakeKey = requires(std::uint64_t seed) { + { Engine::make_key(seed) } -> std::same_as; +}; + +template +concept HasPublicWrappedEngine = requires(Engine engine) { + engine.engine; +}; + +using PublicRepacked = + RandBLAS::rng::RepackedOutput; +static_assert(HasPublicWrappedEngine); +static_assert(std::equality_comparable< + RandBLAS::rng::Philox<4, 32, 10>>); +static_assert(std::equality_comparable< + RandBLAS::rng::RepackedOutput< + RandBLAS::rng::Philox<4, 32, 10>, std::uint16_t>>); + +TEST(RepackedOutput, SplitsEachSourceWordLeastSignificantChunkFirst) { + FixedEngine::ctr_t counter{}; + FixedEngine::key_t key{}; + + using Engine16 = RandBLAS::rng::RepackedOutput; + Engine16::res_t out16{}; + Engine16{}.generate(counter, key, out16); + EXPECT_EQ(out16, (std::array{ + 0xccddu, 0xaabbu, 0x4567u, 0x0123u})); + + using Engine8 = RandBLAS::rng::RepackedOutput; + Engine8::res_t out8{}; + Engine8{}.generate(counter, key, out8); + EXPECT_EQ(out8, (std::array{ + 0xddu, 0xccu, 0xbbu, 0xaau, + 0x67u, 0x45u, 0x23u, 0x01u})); +} + +TEST(RepackedOutput, RepackingPhiloxHasExactDirectAndNestedResults) { + using Base = RandBLAS::rng::Philox<4, 32, 10>; + using Direct16 = RandBLAS::rng::RepackedOutput; + using Direct8 = RandBLAS::rng::RepackedOutput; + using Nested8 = RandBLAS::rng::RepackedOutput; + + Base::ctr_t counter{}; + Base::key_t key{}; + Direct16::res_t out16{}; + Direct8::res_t out8{}; + Nested8::res_t nested8{}; + Direct16{}.generate(counter, key, out16); + Direct8{}.generate(counter, key, out8); + Nested8{}.generate(counter, key, nested8); + + EXPECT_EQ(out16, (Direct16::res_t{ + 0xe8d5u, 0x6627u, 0xc58du, 0xe169u, + 0xac4cu, 0xbc57u, 0xdbd8u, 0x9b00u})); + EXPECT_EQ(out8, (Direct8::res_t{ + 0xd5u, 0xe8u, 0x27u, 0x66u, + 0x8du, 0xc5u, 0x69u, 0xe1u, + 0x4cu, 0xacu, 0x57u, 0xbcu, + 0xd8u, 0xdbu, 0x00u, 0x9bu})); + EXPECT_EQ(nested8, out8); +} + +TEST(RepackedOutput, PreservesBlockBitsCounterAndKeyTypes) { + using Base = RandBLAS::rng::Philox<4, 32, 10>; + using Repacked = RandBLAS::rng::RepackedOutput; + + static_assert(std::same_as); + static_assert(std::same_as); + static_assert(std::tuple_size_v * + std::numeric_limits::digits == + std::tuple_size_v * + std::numeric_limits::digits); + SUCCEED(); +} + +TEST(RepackedOutput, EqualWidthIsAnIdentityAdaptor) { + using Identity = + RandBLAS::rng::RepackedOutput; + FixedEngine::ctr_t counter{}; + FixedEngine::key_t key{}; + Identity::res_t output{}; + + Identity{}.generate(counter, key, output); + + EXPECT_EQ(output, (Identity::res_t{0xaabbccddu, 0x01234567u})); +} + +TEST(RepackedOutput, ForwardsMakeKeyOnlyWhenTheWrappedEngineHasIt) { + using With = RandBLAS::rng::RepackedOutput; + using Without = + RandBLAS::rng::RepackedOutput; + static_assert(HasMakeKey); + static_assert(!HasMakeKey); + EXPECT_EQ(With::make_key(UINT64_C(0x0123456789abcdef)), + (With::key_t{{0x89abcdefu}})); +} + +TEST(RepackedOutput, RejectsInvalidWordTypesAndWidths) { + static_assert(CanRepack); + static_assert(CanRepack); + static_assert(CanRepack); + static_assert(!CanRepack); + static_assert(!CanRepack); + + // Standard unsigned integer widths cannot express these two cases on the + // supported hosts, so exercise the width predicate directly. + static_assert(!RandBLAS::rng::detail::valid_repacking_widths<32, 12>); + static_assert(!RandBLAS::rng::detail::valid_repacking_widths<24, 8>); + SUCCEED(); +} + +} // namespace diff --git a/test/basic_rng/test_rng_state.cc b/test/basic_rng/test_rng_state.cc new file mode 100644 index 00000000..c5792d6c --- /dev/null +++ b/test/basic_rng/test_rng_state.cc @@ -0,0 +1,203 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace { + +class OpaqueCounter { +public: + constexpr void advance(std::uint64_t blocks) noexcept { + value_ += blocks; + } + + friend constexpr bool operator==(OpaqueCounter const&, + OpaqueCounter const&) = default; + +private: + std::uint64_t value_ = 0; + friend struct OpaqueEngine; +}; + +struct OpaqueEngine { + using ctr_t = OpaqueCounter; + using key_t = std::array; + using res_t = std::array; + + static constexpr key_t make_key(std::uint64_t seed) noexcept { + return {static_cast(seed)}; + } + + constexpr void generate(ctr_t const& counter, key_t const& key, + res_t& output) const noexcept { + output[0] = static_cast(counter.value_); + output[1] = static_cast(counter.value_ >> 32) ^ key[0]; + } + + friend constexpr bool operator==(OpaqueEngine const&, + OpaqueEngine const&) = default; +}; + +struct EngineWithoutMakeKey { + using ctr_t = OpaqueEngine::ctr_t; + using key_t = OpaqueEngine::key_t; + using res_t = OpaqueEngine::res_t; + + constexpr void generate(ctr_t const& counter, key_t const& key, + res_t& output) const noexcept { + OpaqueEngine{}.generate(counter, key, output); + } +}; + +static_assert(RandBLAS::rng::CounterBasedEngine); +using OpaqueState = RandBLAS::RNGState; + +template +concept HasPublicStateData = requires(state_t state) { + state.counter; + state.key; + state.engine; +}; + +static_assert(RandBLAS::GeneratorState); +static_assert(HasPublicStateData); +static_assert(std::equality_comparable); +static_assert(std::convertible_to); +static_assert(!std::uniform_random_bit_generator); +static_assert(!std::uniform_random_bit_generator< + RandBLAS::RNGState>); + +TEST(RNGState, SupportsAllApprovedConstructionForms) { + using State = RandBLAS::RNGState; + State default_state; + State scalar_seeded(UINT64_C(0x0123456789abcdef)); + OpaqueEngine::key_t explicit_key{UINT32_C(0x31415926)}; + State key_seeded(explicit_key); + OpaqueCounter counter; + counter.advance(UINT64_C(0x100000002)); + State explicit_state(counter, explicit_key); + + OpaqueEngine::res_t output{}; + default_state.generate(output); + EXPECT_EQ(output, (OpaqueEngine::res_t{0, 0})); + scalar_seeded.generate(output); + EXPECT_EQ(output, (OpaqueEngine::res_t{0, UINT32_C(0x89abcdef)})); + key_seeded.generate(output); + EXPECT_EQ(output, (OpaqueEngine::res_t{0, UINT32_C(0x31415926)})); + explicit_state.generate(output); + EXPECT_EQ(output, (OpaqueEngine::res_t{2, UINT32_C(0x31415927)})); + + static_assert(!std::constructible_from< + RandBLAS::RNGState, std::uint64_t>); +} + +TEST(RNGState, HasRuleOfZeroValueSemanticsAndEquality) { + using State = RandBLAS::RNGState; + static_assert(std::copyable); + static_assert(std::movable); + static_assert(std::is_copy_assignable_v); + static_assert(std::is_move_assignable_v); + + State original(UINT64_C(0x12345678)); + original.advance(19); + State copied = original; + EXPECT_EQ(copied, original); + + State assigned; + assigned = copied; + EXPECT_EQ(assigned, original); + + State moved = std::move(copied); + EXPECT_EQ(moved, original); + State move_assigned; + move_assigned = std::move(assigned); + EXPECT_EQ(move_assigned, original); +} + +TEST(RNGState, GenerateDoesNotMutateAndAdvanceDelegatesToCounter) { + using State = RandBLAS::RNGState; + State state(UINT64_C(0xa5a5a5a5)); + auto before = state; + State::res_t output{}; + + state.generate(output); + + EXPECT_EQ(state, before); + EXPECT_EQ(output, (State::res_t{0, UINT32_C(0xa5a5a5a5)})); + + OpaqueCounter expected_counter; + expected_counter.advance(UINT64_C(0x100000003)); + state.advance(UINT64_C(0x100000003)); + EXPECT_EQ(state.counter, expected_counter); + state.generate(output); + EXPECT_EQ(output, (State::res_t{3, UINT32_C(0xa5a5a5a4)})); +} + +TEST(RNGState, ExposesItsValueStateAsPublicData) { + OpaqueState state(UINT64_C(0x0123456789abcdef)); + EXPECT_EQ(state.counter, OpaqueCounter{}); + EXPECT_EQ(state.key, + (OpaqueEngine::key_t{UINT32_C(0x89abcdef)})); +} + +TEST(RNGState, RepackedStatePreservesBitsAndBlockAdvancement) { + using BaseEngine = RandBLAS::DefaultRNG; + using RepackedEngine = + RandBLAS::rng::RepackedOutput; + RandBLAS::RNGState base(UINT64_C(0x0123456789abcdef)); + RandBLAS::RNGState repacked( + UINT64_C(0x0123456789abcdef)); + BaseEngine::res_t base_output{}; + RepackedEngine::res_t repacked_output{}; + + base.generate(base_output); + repacked.generate(repacked_output); + for (std::size_t i = 0; i < base_output.size(); ++i) { + EXPECT_EQ(repacked_output[2 * i], + static_cast(base_output[i])); + EXPECT_EQ(repacked_output[2 * i + 1], + static_cast(base_output[i] >> 16)); + } + + base.advance(37); + repacked.advance(37); + EXPECT_EQ(base.counter, repacked.counter); + EXPECT_EQ(base.key, repacked.key); +} + +} // namespace diff --git a/test/basic_rng/test_sampler_regression.cc b/test/basic_rng/test_sampler_regression.cc new file mode 100644 index 00000000..c526e73e --- /dev/null +++ b/test/basic_rng/test_sampler_regression.cc @@ -0,0 +1,202 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace { + +using State = RandBLAS::RNGState<>; + +constexpr std::uint64_t seed = 0x0123456789abcdefULL; + +void expect_state(State const& actual, std::uint32_t counter_word_zero) { + EXPECT_EQ(actual.counter[0], counter_word_zero); + EXPECT_EQ(actual.counter[1], 0u); + EXPECT_EQ(actual.counter[2], 0u); + EXPECT_EQ(actual.counter[3], 0u); + EXPECT_EQ(actual.key[0], 0x89abcdefu); + EXPECT_EQ(actual.key[1], 0x01234567u); +} + +template +void expect_float_bits( + std::vector const& actual, + std::array const& expected +) { + ASSERT_EQ(actual.size(), N); + for (std::size_t i = 0; i < N; ++i) { + EXPECT_EQ(std::bit_cast(actual[i]), expected[i]) + << "mismatch at index " << i; + } +} + +template +void expect_gaussian_values( + std::vector const& actual, + std::array const& expected_bits +) { + ASSERT_EQ(actual.size(), N); + constexpr float eps_scale = 32 * std::numeric_limits::epsilon(); + for (std::size_t i = 0; i < N; ++i) { + float expected = std::bit_cast(expected_bits[i]); + float tolerance = eps_scale * std::max(1.0f, std::abs(expected)); + EXPECT_NEAR(actual[i], expected, tolerance) << "mismatch at index " << i; + } +} + +template +void expect_sparse_case( + RandBLAS::Axis axis, + std::int64_t vec_nnz, + std::array const& expected_rows, + std::array const& expected_cols, + std::array const& expected_value_bits, + std::uint32_t expected_counter +) { + State initial{seed}; + RandBLAS::SparseDist dist{5, 11, vec_nnz, axis}; + std::vector values(dist.full_nnz); + std::vector rows(dist.full_nnz); + std::vector cols(dist.full_nnz); + std::int64_t nnz = 0; + + auto next = RandBLAS::fill_sparse_unpacked( + dist, dist.n_rows, dist.n_cols, 0, 0, nnz, + values.data(), rows.data(), cols.data(), initial + ); + + ASSERT_EQ(nnz, static_cast(N)); + for (std::size_t i = 0; i < N; ++i) { + EXPECT_EQ(rows[i], expected_rows[i]) << "row mismatch at index " << i; + EXPECT_EQ(cols[i], expected_cols[i]) << "column mismatch at index " << i; + EXPECT_EQ(std::bit_cast(values[i]), expected_value_bits[i]) + << "value mismatch at index " << i; + } + expect_state(next, expected_counter); +} + +TEST(SamplerRegression, DenseUniformDefaultStream) { + State initial{seed}; + RandBLAS::DenseDist dist{3, 7, RandBLAS::ScalarDist::Uniform}; + std::vector values(21); + + auto next = RandBLAS::fill_dense(dist, values.data(), initial); + + constexpr std::array expected{ + 0xbf7854bbu, 0xbf4a7a6eu, 0x3e8f19beu, 0x3fd435c6u, + 0xbf8e649au, 0x3f8e72b0u, 0x3faf5ee4u, 0xbde318f9u, + 0xbfd9d3bcu, 0x3f8a747fu, 0x3eb80719u, 0xbf912a10u, + 0x3fdd4273u, 0x3de59866u, 0xbfa5fdceu, 0xbf1f001eu, + 0xbf877332u, 0xbf8fac99u, 0x3f9760a8u, 0xbf22c0e8u, + 0xbf47b7a7u + }; + expect_float_bits(values, expected); + expect_state(next, 6u); +} + +TEST(SamplerRegression, DenseGaussianDefaultStream) { + State initial{seed}; + RandBLAS::DenseDist dist{3, 7, RandBLAS::ScalarDist::Gaussian}; + std::vector values(21); + + auto next = RandBLAS::fill_dense(dist, values.data(), initial); + + constexpr std::array expected{ + 0xbf350b89u, 0xbe0a45ffu, 0x3f16e3a3u, 0x3f87d977u, + 0xbfadf224u, 0xbf26bf29u, 0x3f0497a1u, 0xbe6dd4e4u, + 0x3f91d018u, 0x3ffbe4e1u, 0xbf4fc3d1u, 0xbf856f0eu, + 0xbf0d0d9cu, 0x3ccb308du, 0xbee48253u, 0xbee2ab01u, + 0xbf54ee6bu, 0xbe9ac357u, 0x3f08df34u, 0xbeb11d6du, + 0xbed820b3u + }; + expect_gaussian_values(values, expected); + expect_state(next, 6u); +} + +TEST(SamplerRegression, SparseShortAxisDefaultStream) { + constexpr std::array rows{ + 0, 2, 3, 0, 2, 3, 0, 1, 4, 1, 3, 4, 2, 3, 4, 0, 2, + 4, 1, 2, 3, 0, 1, 2, 0, 1, 3, 0, 2, 4, 0, 3, 4 + }; + constexpr std::array cols{ + 0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, + 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10 + }; + constexpr std::array bits{ + 0x3f800000u, 0x3f800000u, 0xbf800000u, 0x3f800000u, + 0x3f800000u, 0x3f800000u, 0x3f800000u, 0x3f800000u, + 0xbf800000u, 0xbf800000u, 0x3f800000u, 0xbf800000u, + 0xbf800000u, 0x3f800000u, 0xbf800000u, 0xbf800000u, + 0xbf800000u, 0xbf800000u, 0x3f800000u, 0xbf800000u, + 0xbf800000u, 0xbf800000u, 0xbf800000u, 0x3f800000u, + 0xbf800000u, 0x3f800000u, 0x3f800000u, 0x3f800000u, + 0xbf800000u, 0x3f800000u, 0xbf800000u, 0x3f800000u, + 0xbf800000u + }; + expect_sparse_case(RandBLAS::Axis::Short, 3, rows, cols, bits, 33u); +} + +TEST(SamplerRegression, SparseLongAxisDefaultStream) { + constexpr std::array rows{ + 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4 + }; + constexpr std::array cols{ + 2, 10, 0, 3, 8, 2, 4, 9, 0, 5, 1, 5, 9 + }; + constexpr std::array bits{ + 0x3fb504f3u, 0x3f800000u, 0x3f800000u, 0x3f800000u, + 0x3f800000u, 0xbf800000u, 0x3f800000u, 0x3f800000u, + 0xbfb504f3u, 0x3f800000u, 0xbf800000u, 0x3f800000u, + 0xbf800000u + }; + expect_sparse_case(RandBLAS::Axis::Long, 3, rows, cols, bits, 15u); +} + +TEST(SamplerRegression, SparseOneNonzeroDefaultStream) { + constexpr std::array rows{2, 2, 0, 0, 1, 2, 4, 1, 2, 4, 1}; + constexpr std::array cols{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; + constexpr std::array bits{ + 0x3f800000u, 0xbf800000u, 0x3f800000u, 0x3f800000u, + 0x3f800000u, 0x3f800000u, 0xbf800000u, 0x3f800000u, + 0x3f800000u, 0xbf800000u, 0xbf800000u + }; + expect_sparse_case(RandBLAS::Axis::Short, 1, rows, cols, bits, 11u); +} + +} // namespace diff --git a/test/basic_rng/test_word_array.cc b/test/basic_rng/test_word_array.cc new file mode 100644 index 00000000..13cb830c --- /dev/null +++ b/test/basic_rng/test_word_array.cc @@ -0,0 +1,128 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#include + +#include + +#include +#include +#include + +namespace { + +using Array32x4 = RandBLAS::rng::WordArray; +using Array64x2 = RandBLAS::rng::WordArray; + +static_assert(std::is_trivially_copyable_v); +static_assert(Array32x4::static_size == 4); +static_assert(Array64x2::static_size == 2); + +TEST(WordArray, ValueInitializesAndSupportsIndexedObservation) { + Array32x4 value{}; + + EXPECT_EQ(value.size(), 4u); + for (std::size_t i = 0; i < value.size(); ++i) { + EXPECT_EQ(value[i], 0u); + } + + value[2] = 17u; + Array32x4 copy = value; + EXPECT_EQ(copy, value); + EXPECT_EQ(copy[2], 17u); +} + +TEST(WordArray, ZeroAdvanceDoesNotChangeValue) { + Array32x4 value{{1u, 2u, 3u, 4u}}; + auto expected = value; + + value.advance(0); + + EXPECT_EQ(value, expected); +} + +TEST(WordArray, AdvancesWithCarryFromLeastSignificantWord) { + Array32x4 value{{0xffffffffu, 7u, 9u, 11u}}; + + value.advance(2); + + EXPECT_EQ(value, (Array32x4{{1u, 8u, 9u, 11u}})); +} + +TEST(WordArray, PropagatesCarryThroughMultipleWords) { + Array32x4 value{{0xffffffffu, 0xffffffffu, 7u, 9u}}; + + value.advance(1); + + EXPECT_EQ(value, (Array32x4{{0u, 0u, 8u, 9u}})); +} + +TEST(WordArray, AddsAllBitsOfUint64To32BitWords) { + Array32x4 value{}; + + value.advance(0x0000000100000001ULL); + + EXPECT_EQ(value, (Array32x4{{1u, 1u, 0u, 0u}})); +} + +TEST(WordArray, CombinesAmountWordsWithExistingCarry) { + Array32x4 value{{0xffffffffu, 0u, 0u, 0u}}; + + value.advance(0x0000000100000001ULL); + + EXPECT_EQ(value, (Array32x4{{0u, 2u, 0u, 0u}})); +} + +TEST(WordArray, WrapsAtFullWidth) { + constexpr auto max = std::numeric_limits::max(); + Array32x4 value{{max, max, max, max}}; + + value.advance(1); + + EXPECT_EQ(value, Array32x4{}); +} + +TEST(WordArray, Advances64BitWordsWithCarry) { + constexpr auto max = std::numeric_limits::max(); + Array64x2 value{{max, 5u}}; + + value.advance(2); + + EXPECT_EQ(value, (Array64x2{{1u, 6u}})); +} + +TEST(WordArray, Wraps64BitWordsAtFullWidth) { + constexpr auto max = std::numeric_limits::max(); + Array64x2 value{{max, max}}; + + value.advance(1); + + EXPECT_EQ(value, Array64x2{}); +} + +} // namespace diff --git a/test/datastructures/test_coo_matrix.cc b/test/datastructures/test_coo_matrix.cc index 715fed32..9215ea5c 100644 --- a/test/datastructures/test_coo_matrix.cc +++ b/test/datastructures/test_coo_matrix.cc @@ -47,9 +47,9 @@ using RandBLAS::SignedInteger; #endif -template +template void sparseskop_to_dense( - RandBLAS::SparseSkOp &S0, + RandBLAS::SparseSkOp &S0, T *mat, Layout layout ) { @@ -154,8 +154,10 @@ class TestCOO : public ::testing::Test { // arrange RandBLAS::DenseDist D(n, n); std::vector buff(n*n); - fill_dense(D, buff.data(), {0}); - iid_sparsify_random_dense(n, n, Layout::ColMajor, buff.data(), prob_zero, {94}); + fill_dense(D, buff.data(), RandBLAS::DefaultRNGState{0}); + iid_sparsify_random_dense( + n, n, Layout::ColMajor, buff.data(), prob_zero, + RandBLAS::DefaultRNGState{94}); std::vector perm(n); for (i = 0; i < n; ++i) { buff[i + i*n] = 2*(i+1); // diagonal is 2, 4, ..., 2*n @@ -312,7 +314,8 @@ class Test_SkOp_to_COO : public ::testing::Test { template void sparse_skop_to_coo(int64_t d, int64_t m, int64_t key_index, int64_t nnz_index, Axis major_axis) { RandBLAS::SparseDist D(d, m, vec_nnzs[nnz_index], major_axis); - RandBLAS::SparseSkOp S(D, keys[key_index]); + RandBLAS::SparseSkOp S( + D, RandBLAS::DefaultRNGState{keys[key_index]}); fill_sparse(S); auto A = RandBLAS::sparse::coo_view_of_skop(S); diff --git a/test/datastructures/test_denseskop.cc b/test/datastructures/test_denseskop.cc index 6cdd913c..f89ed0bf 100644 --- a/test/datastructures/test_denseskop.cc +++ b/test/datastructures/test_denseskop.cc @@ -39,42 +39,41 @@ #include #include #include +#include // Fill a random matrix and truncate at the end of each row so that each row starts with a fresh counter. -template +template static void fill_dense_rmat_trunc( T* mat, int64_t n_rows, int64_t n_cols, - const RandBLAS::RNGState & seed + const State &seed ) { - - RNG rng; - typename RNG::ctr_type c = seed.counter; - typename RNG::key_type k = seed.key; + State work = seed; + constexpr int block_size = std::tuple_size_v; int ind = 0; - int cts = n_cols / RNG::ctr_type::static_size; + int cts = n_cols / block_size; // ^ number of counters per row, where all the random numbers are to be filled in the array. - int res = n_cols % RNG::ctr_type::static_size; + int res = n_cols % block_size; // ^ Number of random numbers to be filled at the end of each row the the last counter of the row for (int i = 0; i < n_rows; i++) { for (int ctr = 0; ctr < cts; ctr++){ - auto rv = OP::generate(rng, c, k); - for (int j = 0; j < RNG::ctr_type::static_size; j++) { + auto rv = OP::generate(work); + for (int j = 0; j < block_size; j++) { mat[ind] = rv[j]; ind++; } - c.incr(); + work.advance(1); } if (res != 0) { for (int j = 0; j < res; j++) { - auto rv = OP::generate(rng, c, k); + auto rv = OP::generate(work); mat[ind] = rv[j]; ind++; } - c.incr(); + work.advance(1); } } } @@ -168,23 +167,23 @@ class TestSubmatGeneration : public ::testing::Test virtual void TearDown(){}; - template + template static void test_colwise_smat_gen( int64_t n_cols, int64_t n_rows, int64_t n_scols, int64_t n_srows, int64_t ptr, - const RandBLAS::RNGState &seed + const State &seed ) { int stride = n_cols / 50; T* mat = new T[n_rows * n_cols]; T* smat = new T[n_srows * n_scols]; - fill_dense_rmat_trunc(mat, n_rows, n_cols, seed); + fill_dense_rmat_trunc(mat, n_rows, n_cols, seed); int ind = 0; // used for indexing smat when comparing to rmat for (int nptr = ptr; nptr < n_cols*(n_rows-n_srows-1); nptr += stride*n_cols) { // ^ Loop through various pointer locations.- goes down the random matrix by amount stride. - RandBLAS::dense::fill_dense_submat_impl(n_cols, smat, n_srows, n_scols, nptr, seed); + RandBLAS::dense::fill_dense_submat_impl(n_cols, smat, n_srows, n_scols, nptr, seed); ind = 0; for (int i = 0; i + template static void test_rowwise_smat_gen( int64_t n_cols, int64_t n_rows, int64_t n_scols, int64_t n_srows, int64_t ptr, - const RandBLAS::RNGState &seed + const State &seed ) { int stride = n_cols / 50; T* mat = new T[n_rows * n_cols]; T* smat = new T[n_srows * n_scols]; - fill_dense_rmat_trunc(mat, n_rows, n_cols, seed); + fill_dense_rmat_trunc(mat, n_rows, n_cols, seed); int ind = 0; // variable used for indexing smat when comparing to rmat for (int nptr = ptr; nptr < (n_cols - n_scols - 1); nptr += stride) { // ^ Loop through various pointer locations.- goes across the random matrix by amount stride. - RandBLAS::dense::fill_dense_submat_impl(n_cols, smat, n_srows, n_scols, nptr, seed); + RandBLAS::dense::fill_dense_submat_impl(n_cols, smat, n_srows, n_scols, nptr, seed); ind = 0; for (int i = 0; i + template static void test_diag_smat_gen( int64_t n_cols, int64_t n_rows, - const RandBLAS::RNGState &seed + const State &seed ) { T* mat = new T[n_rows * n_cols]; T* smat = new T[n_rows * n_cols]{}; - fill_dense_rmat_trunc(mat, n_rows, n_cols, seed); + fill_dense_rmat_trunc(mat, n_rows, n_cols, seed); int ind = 0; int64_t n_scols = 1; int64_t n_srows = 1; for (int ptr = 0; ptr + n_scols + n_cols*n_srows < n_cols*n_rows; ptr += n_rows+1) { // Loop through the diagonal of the matrix RandBLAS::util::safe_scal(n_srows * n_scols, (T) 0.0, smat); - RandBLAS::dense::fill_dense_submat_impl(n_cols, smat, n_srows, n_scols, ptr, seed); + RandBLAS::dense::fill_dense_submat_impl(n_cols, smat, n_srows, n_scols, ptr, seed); ind = 0; for (int i = 0; i seed(k); - test_colwise_smat_gen(n_cols, n_rows, n_scols, n_srows, ptr, seed); + RandBLAS::DefaultRNGState seed(k); + test_colwise_smat_gen(n_cols, n_rows, n_scols, n_srows, ptr, seed); } } @@ -280,8 +279,8 @@ TEST_F(TestSubmatGeneration, row_wise) int64_t n_scols = 100; int64_t ptr = n_rows + 2; for (int k = 0; k < 3; k++) { - RandBLAS::RNGState seed(k); - test_rowwise_smat_gen(n_cols, n_rows, n_scols, n_srows, ptr, seed); + RandBLAS::DefaultRNGState seed(k); + test_rowwise_smat_gen(n_cols, n_rows, n_scols, n_srows, ptr, seed); } } @@ -290,14 +289,14 @@ TEST_F(TestSubmatGeneration, diag) int64_t n_rows = 100; int64_t n_cols = 2000; for (int k = 0; k < 3; k++) { - RandBLAS::RNGState seed(k); - test_diag_smat_gen(n_cols, n_rows, seed); + RandBLAS::DefaultRNGState seed(k); + test_diag_smat_gen(n_cols, n_rows, seed); } } #if defined(RandBLAS_HAS_OpenMP) -template +template void DenseThreadTest(int64_t m, int64_t n) { int64_t d = m*n; @@ -306,8 +305,8 @@ void DenseThreadTest(int64_t m, int64_t n) { // generate the base state with 1 thread. omp_set_num_threads(1); - RandBLAS::RNGState state(0); - RandBLAS::dense::fill_dense_submat_impl(n, base.data(), m, n, 0, state); + State state(0); + RandBLAS::dense::fill_dense_submat_impl(n, base.data(), m, n, 0, state); std::cerr << "with 1 thread: " << base << std::endl; // run with different numbers of threads, and check that the result is the same @@ -315,7 +314,7 @@ void DenseThreadTest(int64_t m, int64_t n) { for (int i = 2; i <= n_threads; ++i) { std::fill(test.begin(), test.end(), (T) 0.0); omp_set_num_threads(i); - RandBLAS::dense::fill_dense_submat_impl(n, test.data(), m, n, 0, state); + RandBLAS::dense::fill_dense_submat_impl(n, test.data(), m, n, 0, state); std::cerr << "with " << i << " threads: " << test << std::endl; for (int64_t i = 0; i < d; ++i) { EXPECT_FLOAT_EQ( base[i], test[i] ); @@ -325,17 +324,17 @@ void DenseThreadTest(int64_t m, int64_t n) { TEST(TestDenseThreading, UniformPhilox) { for (int i = 0; i < 10; ++i) { - DenseThreadTest(32, 8); - DenseThreadTest(1, 5); - DenseThreadTest(5, 1); + DenseThreadTest(32, 8); + DenseThreadTest(1, 5); + DenseThreadTest(5, 1); } } TEST(TestDenseThreading, GaussianPhilox) { for (int i = 0; i < 10; ++i) { - DenseThreadTest(32, 8); - DenseThreadTest(1, 5); - DenseThreadTest(5, 1); + DenseThreadTest(32, 8); + DenseThreadTest(1, 5); + DenseThreadTest(5, 1); } } #endif @@ -352,12 +351,14 @@ class TestFillAxis : public::testing::Test // make the wide sketching operator RandBLAS::DenseDist D_wide(short_dim, long_dim, distname, major_axis); - RandBLAS::DenseSkOp S_wide(D_wide, seed); + RandBLAS::DenseSkOp S_wide( + D_wide, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_dense(S_wide); // make the tall sketching operator RandBLAS::DenseDist D_tall(long_dim, short_dim, distname, major_axis); - RandBLAS::DenseSkOp S_tall(D_tall, seed); + RandBLAS::DenseSkOp S_tall( + D_tall, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_dense(S_tall); // Sanity check: layouts are opposite. @@ -442,7 +443,7 @@ class TestDenseSkOpStates : public ::testing::Test } } - template + template static void test_compute_next_state( uint32_t key, int64_t n_rows, @@ -460,7 +461,7 @@ class TestDenseSkOpStates : public ::testing::Test auto expect_final_state = RandBLAS::dense::compute_next_state(D, state); auto expect_c = expect_final_state.counter; - for (int i = 0; i < RNG::ctr_type::static_size; i++) { + for (std::size_t i = 0; i < std::tuple_size_v; i++) { ASSERT_EQ(actual_c[i], expect_c[i]); } @@ -483,10 +484,10 @@ TEST_F(TestDenseSkOpStates, concat_tall_with_long_major_axis) { TEST_F(TestDenseSkOpStates, compare_skopless_fill_dense_to_compute_next_state) { for (uint32_t key : {0, 1, 2}) { auto sd = RandBLAS::ScalarDist::Gaussian; - test_compute_next_state(key, 13, 7, sd); - test_compute_next_state(key, 11, 5, sd); - test_compute_next_state(key, 131, 71, sd); - test_compute_next_state(key, 80, 40, sd); - test_compute_next_state(key, 91, 43, sd); + test_compute_next_state(key, 13, 7, sd); + test_compute_next_state(key, 11, 5, sd); + test_compute_next_state(key, 131, 71, sd); + test_compute_next_state(key, 80, 40, sd); + test_compute_next_state(key, 91, 43, sd); } } diff --git a/test/datastructures/test_sparseskop.cc b/test/datastructures/test_sparseskop.cc index faee4004..b8c9c2b1 100644 --- a/test/datastructures/test_sparseskop.cc +++ b/test/datastructures/test_sparseskop.cc @@ -94,9 +94,10 @@ class TestSparseSkOpConstruction : public ::testing::Test template void proper_saso_construction(int64_t d, int64_t m, int64_t key_index, int64_t nnz_index) { - using RNG = SparseSkOp::state_t::generator; + using State = SparseSkOp::state_t; SparseDist D0(d, m, vec_nnzs[nnz_index], Axis::Short); - SparseSkOp S0(D0, keys[key_index]); + SparseSkOp S0( + D0, State{keys[key_index]}); fill_sparse(S0); if (d < m) { check_fixed_nnz_per_col(S0); @@ -107,10 +108,11 @@ class TestSparseSkOpConstruction : public ::testing::Test template void proper_laso_construction(int64_t d, int64_t m, int64_t key_index) { - using RNG = SparseSkOp::state_t::generator; + using State = SparseSkOp::state_t; int64_t vec_nnz = 1; SparseDist D0(d, m, vec_nnz, Axis::Long); - SparseSkOp S0(D0, keys[key_index]); + SparseSkOp S0( + D0, State{keys[key_index]}); fill_sparse(S0); if (d < m) { check_fixed_nnz_per_row(S0); @@ -175,7 +177,7 @@ class TestSparseSkOpConstruction : public ::testing::Test } void unpacked_nosub(const SparseDist &D) { - RNGState s(1); + RandBLAS::DefaultRNGState s(1); SparseSkOp S(D, s); auto expect_next = S.next_state; fill_sparse(S); @@ -221,7 +223,7 @@ class TestSparseSkOpConstruction : public ::testing::Test int64_t ro_s, int64_t co_s, int64_t n_rows_sub, int64_t n_cols_sub, uint32_t key ) { SparseDist D {d, m, vec_nnz, major_axis}; - RNGState seed(key); + RandBLAS::DefaultRNGState seed(key); // Full operator via the no-submatrix path. int64_t full_nnz_out = -1; diff --git a/test/downstream/main.cc b/test/downstream/main.cc index 142fdbbe..4bbf4c64 100644 --- a/test/downstream/main.cc +++ b/test/downstream/main.cc @@ -14,8 +14,11 @@ int main() { double A[m * n]; RandBLAS::DenseDist Dist_A(m, n); - RandBLAS::RNGState state(0); - RandBLAS::fill_dense(Dist_A, A, state); + RandBLAS::DefaultRNGState state(0); + RandBLAS::DefaultRNGState::res_t block{}; + state.generate(block); + state.advance(1); + auto next_state = RandBLAS::fill_dense(Dist_A, A, state); - return 0; + return next_state == state; } diff --git a/test/linops/test_lskge3.cc b/test/linops/test_lskge3.cc index 28874b67..cc4ed6b5 100644 --- a/test/linops/test_lskge3.cc +++ b/test/linops/test_lskge3.cc @@ -51,7 +51,7 @@ class TestLSKGE3 : public ::testing::Test blas::Layout layout ) { DenseDist D(d, m); - DenseSkOp S0(D, seed); + DenseSkOp S0(D, RandBLAS::DefaultRNGState{seed}); if (preallocate) RandBLAS::fill_dense(S0); test_left_apply_submatrix_to_eye(1.0, S0, d, m, 0, 0, layout, 0.0); @@ -65,7 +65,7 @@ class TestLSKGE3 : public ::testing::Test blas::Layout layout ) { DenseDist Dt(m, d); - DenseSkOp S0(Dt, seed); + DenseSkOp S0(Dt, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_dense(S0); test_left_apply_transpose_to_eye(S0, layout); } @@ -84,7 +84,7 @@ class TestLSKGE3 : public ::testing::Test randblas_require(d0 > d); randblas_require(m0 > m); DenseDist D(d0, m0); - DenseSkOp S0(D, seed); + DenseSkOp S0(D, RandBLAS::DefaultRNGState{seed}); test_left_apply_submatrix_to_eye(1.0, S0, d, m, S_ro, S_co, layout, 0.0); } @@ -103,7 +103,7 @@ class TestLSKGE3 : public ::testing::Test randblas_require(m0 > m); randblas_require(n0 > n); DenseDist D(d, m); - DenseSkOp S0(D, seed_S0); + DenseSkOp S0(D, RandBLAS::DefaultRNGState{seed_S0}); test_left_apply_to_submatrix(S0, n, m0, n0, A_ro, A_co, layout); } @@ -304,4 +304,3 @@ TEST_F(TestLSKGE3, submatrix_a_single) blas::Layout::ColMajor ); } - diff --git a/test/linops/test_lskges.cc b/test/linops/test_lskges.cc index ff92b0cd..1de2af65 100644 --- a/test/linops/test_lskges.cc +++ b/test/linops/test_lskges.cc @@ -58,7 +58,8 @@ class TestLSKGES : public ::testing::Test int threads ) { SparseDist D0(d, m, vec_nnzs[nnz_index], major_axis); - SparseSkOp S0(D0, keys[key_index]); + SparseSkOp S0( + D0, RandBLAS::DefaultRNGState{keys[key_index]}); test_left_apply_to_random(1.0, S0, n, 0.0, layout, threads); } @@ -75,7 +76,7 @@ class TestLSKGES : public ::testing::Test ) { int64_t vec_nnz = d0 / 3; // this is actually quite dense. SparseDist D0(d0, m0, vec_nnz, Axis::Short); - SparseSkOp S0(D0, seed); + SparseSkOp S0(D0, RandBLAS::DefaultRNGState{seed}); test_left_apply_submatrix_to_eye(1.0, S0, d1, m1, S_ro, S_co, layout, 0.0); } @@ -90,7 +91,7 @@ class TestLSKGES : public ::testing::Test ) { int64_t vec_nnz = d / 2; SparseDist DS(d, m, vec_nnz, Axis::Short); - SparseSkOp S(DS, key); + SparseSkOp S(DS, RandBLAS::DefaultRNGState{key}); test_left_apply_submatrix_to_eye(alpha, S, d, m, 0, 0, layout, beta); } @@ -106,7 +107,7 @@ class TestLSKGES : public ::testing::Test bool is_saso = (major_axis == Axis::Short); int64_t vec_nnz = (is_saso) ? d/2 : m/2; SparseDist Dt(m, d, vec_nnz, major_axis); - SparseSkOp S0(Dt, key); + SparseSkOp S0(Dt, RandBLAS::DefaultRNGState{key}); test_left_apply_transpose_to_eye(S0, layout); } @@ -127,7 +128,7 @@ class TestLSKGES : public ::testing::Test bool is_saso = (major_axis == Axis::Short); int64_t vec_nnz = (is_saso) ? d/2 : m/2; SparseDist D(d, m, vec_nnz, major_axis); - SparseSkOp S0(D, seed_S0); + SparseSkOp S0(D, RandBLAS::DefaultRNGState{seed_S0}); test_left_apply_to_submatrix(S0, n, m0, n0, A_ro, A_co, layout); } @@ -144,7 +145,7 @@ class TestLSKGES : public ::testing::Test bool is_saso = (major_axis == Axis::Short); int64_t vec_nnz = (is_saso) ? d/2 : m/2; SparseDist D(d, m, vec_nnz, major_axis); - SparseSkOp S0(D, seed_S0); + SparseSkOp S0(D, RandBLAS::DefaultRNGState{seed_S0}); test_left_apply_to_transposed(S0, n, layout); } }; @@ -599,12 +600,12 @@ class TestLSKGES_SubmatrixPath : public ::testing::Test { if (opS == blas::Op::NoTrans) { big_rows = d1 + S_ro + 1; big_cols = m1 + S_co + 2; } else { big_rows = m1 + S_ro + 1; big_cols = d1 + S_co + 2; } SparseDist D0 {big_rows, big_cols, vec_nnz, major_axis}; - RandBLAS::RNGState seed((uint32_t) 7); + RandBLAS::DefaultRNGState seed((uint32_t) 7); // Dense input A: op(A) is m1-by-n with opA = NoTrans, so A is m1-by-n. int64_t lda = (layout == blas::Layout::ColMajor) ? m1 : n; std::vector A(m1 * n); - RandBLAS::RNGState a_state((uint32_t) 99); + RandBLAS::DefaultRNGState a_state((uint32_t) 99); RandBLAS::DenseDist DA {m1, n}; RandBLAS::fill_dense(DA, A.data(), a_state); diff --git a/test/linops/test_rskge3.cc b/test/linops/test_rskge3.cc index d47d0ee3..59cffef7 100644 --- a/test/linops/test_rskge3.cc +++ b/test/linops/test_rskge3.cc @@ -52,7 +52,7 @@ class TestRSKGE3 : public ::testing::Test Layout layout ) { DenseDist D(m, d); - DenseSkOp S0(D, seed); + DenseSkOp S0(D, RandBLAS::DefaultRNGState{seed}); if (preallocate) RandBLAS::fill_dense(S0); test_right_apply_submatrix_to_eye(1.0, S0, m, d, 0, 0, layout, 0.0, 0); @@ -66,7 +66,7 @@ class TestRSKGE3 : public ::testing::Test Layout layout ) { DenseDist Dt(d, m); - DenseSkOp S0(Dt, seed); + DenseSkOp S0(Dt, RandBLAS::DefaultRNGState{seed}); test_right_apply_transpose_to_eye(S0, layout); } @@ -82,7 +82,7 @@ class TestRSKGE3 : public ::testing::Test Layout layout ) { DenseDist D(m0, d0); - DenseSkOp S0(D, seed); + DenseSkOp S0(D, RandBLAS::DefaultRNGState{seed}); test_right_apply_submatrix_to_eye(1.0, S0, m, d, S_ro, S_co, layout, 0.0, 0); } @@ -99,7 +99,7 @@ class TestRSKGE3 : public ::testing::Test Layout layout ) { DenseDist D(n, d); - DenseSkOp S0(D, seed_S0); + DenseSkOp S0(D, RandBLAS::DefaultRNGState{seed_S0}); test_right_apply_to_submatrix(S0, m, m0, n0, A_ro, A_co, layout); } diff --git a/test/linops/test_rskges.cc b/test/linops/test_rskges.cc index da1ab06a..1bcebcb6 100644 --- a/test/linops/test_rskges.cc +++ b/test/linops/test_rskges.cc @@ -56,7 +56,7 @@ class TestRSKGES : public ::testing::Test Layout layout ) { SparseDist D(m, d, vec_nnz, major_axis); - SparseSkOp S0(D, seed); + SparseSkOp S0(D, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_sparse(S0); test_right_apply_submatrix_to_eye(1.0, S0, m, d, 0, 0, layout, 0.0, 0); } @@ -73,7 +73,8 @@ class TestRSKGES : public ::testing::Test int threads ) { SparseDist D(n, d, vec_nnzs[nnz_index], major_axis); - SparseSkOp S0(D, keys[key_index]); + SparseSkOp S0( + D, RandBLAS::DefaultRNGState{keys[key_index]}); RandBLAS::fill_sparse(S0); test_right_apply_to_random(1.0, S0, m, layout, 0.0, threads); } @@ -93,7 +94,7 @@ class TestRSKGES : public ::testing::Test randblas_require(n0 >= n1); int64_t vec_nnz = d0 / 3; // this is actually quite dense. SparseDist D0(n0, d0, vec_nnz, RandBLAS::Axis::Short); - SparseSkOp S0(D0, seed); + SparseSkOp S0(D0, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_sparse(S0); test_right_apply_submatrix_to_eye(1.0, S0, n1, d1, S_ro, S_co, layout, 0.0, 0); } @@ -410,12 +411,12 @@ class TestRSKGES_SubmatrixPath : public ::testing::Test { if (opS == blas::Op::NoTrans) { big_rows = n + S_ro + 1; big_cols = d1 + S_co + 2; } else { big_rows = d1 + S_ro + 1; big_cols = n + S_co + 2; } SparseDist D0 {big_rows, big_cols, vec_nnz, major_axis}; - RandBLAS::RNGState seed((uint32_t) 7); + RandBLAS::DefaultRNGState seed((uint32_t) 7); // Dense input A: op(A) is m-by-n with opA = NoTrans, so A is m-by-n. int64_t lda = (layout == blas::Layout::ColMajor) ? m : n; std::vector A(m * n); - RandBLAS::RNGState a_state((uint32_t) 99); + RandBLAS::DefaultRNGState a_state((uint32_t) 99); RandBLAS::DenseDist DA {m, n}; RandBLAS::fill_dense(DA, A.data(), a_state); diff --git a/test/linops/test_sketch_sparse.cc b/test/linops/test_sketch_sparse.cc index 2afa4d7f..5eb4b14f 100644 --- a/test/linops/test_sketch_sparse.cc +++ b/test/linops/test_sketch_sparse.cc @@ -212,7 +212,7 @@ class TestLSKSP3 : public ::testing::Test template static void sketch_eye(uint32_t seed, int64_t m, int64_t d, bool preallocate, Layout layout) { DenseDist D(d, m); - DenseSkOp S0(D, seed); + DenseSkOp S0(D, RandBLAS::DefaultRNGState{seed}); if (preallocate) RandBLAS::fill_dense(S0); test_left_submat_sketch_of_eye(1.0, S0, d, m, 0, 0, layout, 0.0); @@ -221,7 +221,7 @@ class TestLSKSP3 : public ::testing::Test template static void transpose_S(uint32_t seed, int64_t m, int64_t d, Layout layout) { DenseDist Dt(m, d); - DenseSkOp S0(Dt, seed); + DenseSkOp S0(Dt, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_dense(S0); test_left_transposed_sketch_of_eye(S0, layout); } @@ -240,7 +240,7 @@ class TestLSKSP3 : public ::testing::Test randblas_require(d0 > d); randblas_require(m0 > m); DenseDist D(d0, m0); - DenseSkOp S0(D, seed); + DenseSkOp S0(D, RandBLAS::DefaultRNGState{seed}); test_left_submat_sketch_of_eye(1.0, S0, d, m, S_ro, S_co, layout, 0.0); } @@ -388,7 +388,7 @@ class TestRSKSP3 : public ::testing::Test template static void sketch_eye(uint32_t seed, int64_t m, int64_t d, bool preallocate, Layout layout) { DenseDist D(m, d); - DenseSkOp S0(D, seed); + DenseSkOp S0(D, RandBLAS::DefaultRNGState{seed}); if (preallocate) RandBLAS::fill_dense(S0); test_right_submat_sketch_of_eye(1.0, S0, m, d, 0, 0, layout, 0.0); @@ -397,7 +397,7 @@ class TestRSKSP3 : public ::testing::Test template static void transpose_S(uint32_t seed, int64_t m, int64_t d, Layout layout) { DenseDist Dt(d, m); - DenseSkOp S0(Dt, seed); + DenseSkOp S0(Dt, RandBLAS::DefaultRNGState{seed}); test_right_transposed_sketch_of_eye(S0, layout); } @@ -413,7 +413,7 @@ class TestRSKSP3 : public ::testing::Test Layout layout ) { DenseDist D(m0, d0); - DenseSkOp S0(D, seed); + DenseSkOp S0(D, RandBLAS::DefaultRNGState{seed}); test_right_submat_sketch_of_eye(1.0, S0, m, d, S_ro, S_co, layout, 0.0); } diff --git a/test/linops/test_sketch_symmetric.cc b/test/linops/test_sketch_symmetric.cc index 4156c3a5..ecd77a61 100644 --- a/test/linops/test_sketch_symmetric.cc +++ b/test/linops/test_sketch_symmetric.cc @@ -94,7 +94,7 @@ class TestSketchSymmetric : public ::testing::Test { std::vector A(lda*lda, 0.0); random_symmetric_mat(n, A.data(), lda, RNGState(seed_a)); DenseDist D(rows_out, cols_out, ScalarDist::Uniform, major_axis); - DenseSkOp S(D, seed_skop); + DenseSkOp S(D, RandBLAS::DefaultRNGState{seed_skop}); RandBLAS::fill_dense(S); int64_t lds = (S.layout == Layout::RowMajor) ? cols_out : rows_out; int64_t ldb = lds; @@ -126,7 +126,7 @@ class TestSketchSymmetric : public ::testing::Test { std::vector A(lda*lda, 0.0); random_symmetric_mat(n, A.data(), lda, RNGState(seed_a)); DenseDist D(rows_out, cols_out, ScalarDist::Uniform, major_axis); - DenseSkOp S(D, seed_skop); + DenseSkOp S(D, RandBLAS::DefaultRNGState{seed_skop}); RandBLAS::fill_dense(S); int64_t lds_init, ldb; Layout layout_B; @@ -177,7 +177,7 @@ class TestSketchSymmetric : public ::testing::Test { 0.0, 0.0, 3.0 }; DenseDist D(d, n, ScalarDist::Uniform, Axis::Short); - DenseSkOp S(D, 42); + DenseSkOp S(D, RandBLAS::DefaultRNGState{42}); RandBLAS::fill_dense(S); std::vector B(d * n, 0.0); try { diff --git a/test/linops/test_sketch_vector.cc b/test/linops/test_sketch_vector.cc index 4dc5040f..565437dc 100644 --- a/test/linops/test_sketch_vector.cc +++ b/test/linops/test_sketch_vector.cc @@ -66,7 +66,8 @@ class TestSketchVector : public ::testing::Test x[incx*i] = 1.0; RandBLAS::DenseDist D(d, m); - RandBLAS::DenseSkOp S(D, seed); + RandBLAS::DenseSkOp S( + D, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_dense(S); int64_t lds = (S.layout == blas::Layout::RowMajor) ? m : d; @@ -100,7 +101,8 @@ class TestSketchVector : public ::testing::Test x[incx*i] = 1.0; RandBLAS::DenseDist D(m, d); - RandBLAS::DenseSkOp S(D, seed); + RandBLAS::DenseSkOp S( + D, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_dense(S); int64_t lds = (S.layout == blas::Layout::RowMajor) ? d : m; @@ -137,9 +139,11 @@ class TestSketchVector : public ::testing::Test // Generate wide and tall sketching operator using same seed RandBLAS::DenseDist D_wide(d, m); RandBLAS::DenseDist D_tall(m, d); - RandBLAS::DenseSkOp S_wide(D_wide, seed); + RandBLAS::DenseSkOp S_wide( + D_wide, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_dense(S_wide); - RandBLAS::DenseSkOp S_tall(D_tall, seed); + RandBLAS::DenseSkOp S_tall( + D_tall, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_dense(S_tall); // Perform wide sketch with Op::NoTrans and tall sketch with Op::Trans. Should be the same operation @@ -173,7 +177,8 @@ class TestSketchVector : public ::testing::Test // Generate tall sketching operator RandBLAS::DenseDist D(d, m); - RandBLAS::DenseSkOp S(D, seed); + RandBLAS::DenseSkOp S( + D, RandBLAS::DefaultRNGState{seed}); RandBLAS::fill_dense(S); int64_t lds = (S.layout == blas::Layout::RowMajor) ? m : d; diff --git a/test/linops/test_spgemm.cc b/test/linops/test_spgemm.cc index 2db69fa4..9468f654 100644 --- a/test/linops/test_spgemm.cc +++ b/test/linops/test_spgemm.cc @@ -36,6 +36,7 @@ #include "RandBLAS/testing/comparison.hh" #include #include +#include #include using namespace RandBLAS::sparse_data; @@ -122,15 +123,15 @@ class TestSpGEMM : public ::testing::Test { // Error model: |C_actual - C_ref| <= |alpha| * k * 2*eps * |A_dense| * |B_dense| + |beta| * eps * |C_orig| // We compute the error bound via gemm on absolute values. T eps = std::numeric_limits::epsilon(); - T err_alpha = abs(alpha) * k * 2 * eps; - T err_beta = abs(beta) * eps; + T err_alpha = std::abs(alpha) * k * 2 * eps; + T err_beta = std::abs(beta) * eps; std::vector A_abs(rows_A * cols_A); std::vector B_abs(k * n); for (int64_t i = 0; i < (int64_t)A_abs.size(); ++i) - A_abs[i] = abs(A_dense[i]); + A_abs[i] = std::abs(A_dense[i]); for (int64_t i = 0; i < (int64_t)B_abs.size(); ++i) - B_abs[i] = abs(B_dense[i]); + B_abs[i] = std::abs(B_dense[i]); // Start error bound with |beta| * eps * |C_orig| std::vector E(m * n, 0.0); @@ -140,7 +141,7 @@ class TestSpGEMM : public ::testing::Test { // we need |C_orig|. Re-generate it. auto C_orig = std::get<0>(random_matrix(m, n, RandBLAS::RNGState(42))); for (int64_t i = 0; i < m * n; ++i) - E[i] = abs(C_orig[i]); + E[i] = std::abs(C_orig[i]); } blas::gemm(layout, opA, Op::NoTrans, m, n, k, diff --git a/test/meta/test_rng_stream.cc b/test/meta/test_rng_stream.cc new file mode 100644 index 00000000..6b4089d3 --- /dev/null +++ b/test/meta/test_rng_stream.cc @@ -0,0 +1,97 @@ +// Copyright, 2026. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) 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. +// +// (3) 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 OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +#include + +#include + +#include +#include +#include + +namespace { + +struct SequenceState { + using res_t = std::array; + + res_t first{}; + std::uint64_t block{}; + + void generate(res_t& output) const { + output = { + static_cast(first[0] + 2 * block), + static_cast(first[1] + 2 * block) + }; + } + + void advance(std::uint64_t blocks) { block += blocks; } +}; + +static_assert(RandBLAS::GeneratorState); + +using RNGStream = RandBLAS::testing::detail::RNGStream; + +TEST(RNGStream, NextWordBuffersOneBlockAndAdvancesOnRefill) { + RNGStream stream(SequenceState{{10, 20}}); + + EXPECT_EQ(stream.get_state().block, 0u); + EXPECT_EQ(stream.next_word(), 10u); + EXPECT_EQ(stream.get_state().block, 1u); + EXPECT_EQ(stream.next_word(), 20u); + EXPECT_EQ(stream.get_state().block, 1u); + EXPECT_EQ(stream.next_word(), 12u); + EXPECT_EQ(stream.get_state().block, 2u); +} + +TEST(RNGStream, GaussianCachesTheSecondValue) { + constexpr std::uint32_t angle_word = UINT32_C(0x243f6a88); + constexpr std::uint32_t radius_word = UINT32_C(0x85a308d3); + RNGStream stream(SequenceState{{angle_word, radius_word}}); + auto expected = RandBLAS::rng::boxmuller(angle_word, radius_word); + + EXPECT_EQ(stream.gaussian(), expected[0]); + EXPECT_EQ(stream.get_state().block, 1u); + EXPECT_EQ(stream.gaussian(), expected[1]); + EXPECT_EQ(stream.get_state().block, 1u); +} + +TEST(RNGStream, UniformAndGeometricUseScalarConversions) { + constexpr std::uint32_t word = UINT32_C(0x243f6a88); + RNGStream uniform_stream(SequenceState{{word, 0}}); + EXPECT_EQ(uniform_stream.uniform_01(), + RandBLAS::rng::u01(word)); + + RNGStream geometric_stream(SequenceState{{word, 0}}); + double log_1_minus_p = std::log(0.75); + double u = RandBLAS::rng::u01(word); + auto expected = static_cast( + std::floor(std::log(1.0 - u) / log_1_minus_p)); + EXPECT_EQ(geometric_stream.geometric(log_1_minus_p), expected); +} + +} // namespace diff --git a/test/test_io.cc b/test/test_io.cc index 2a12ff5d..9e025270 100644 --- a/test/test_io.cc +++ b/test/test_io.cc @@ -69,6 +69,6 @@ TEST_F(TestIO, test_rngstate_insertion_operator) { int64_t two_pow_33 = static_cast(1) << 33; int64_t two_pow_61 = static_cast(1) << 61; RandBLAS::RNGState state(two_pow_33); // expect {0, 2} - state.counter.incr(two_pow_61); // expect {0, 536870912, 0, 0} + state.advance(two_pow_61); // expect {0, 536870912, 0, 0} std::cout << state << std::endl; }