Skip to content

vPhilox theory.md

Latest commit

 

History

History
executable file
·
107 lines (73 loc) · 6.97 KB

vPhilox theory.md

File metadata and controls

executable file
·
107 lines (73 loc) · 6.97 KB
Raw

Here is the theoretical foundation of Counter-Based PRNGs (CBRNGs) and the Philox algorithm—covering the exact mathematics, hardware-level mechanics, vectorization theory, and float-conversion arithmetic required to build vphilox.

1. Counter-Based PRNG Theory (CBRNG)

Traditional stateful PRNGs (e.g., std::mt19937 or LCGs) evaluate a sequence via recurrence:

$$S_{i+1} = T(S_i) \quad \implies \quad R_i = g(S_i)$$
Because state $S_i$ depends strictly on $S_{i-1}$, generating random number $N$ requires running the state transition $N$ times.
In contrast, a Counter-Based PRNG treats random number generation as a bijection (permutation) over an $n$-bit state domain:

$$R = f(K, C)$$

  • $K$ (Key): A user-supplied seed or sequence modifier (e.g., thread ID or tree index).
  • $C$ (Counter): A stateless integer index ($0, 1, 2, 3, \dots$).
  • $f$: A cryptographic-like mixing function that maps $(K, C) \to R$ deterministically.

Theoretical Advantages

  1. $O(1)$ Random Seeking: You can jump to any element $C_{100,000,000}$ in constant time $O(1)$ without computing the preceding steps.
  2. Embarrassingly Parallel: Threads $T_0, T_1, \dots, T_k$ compute outputs for counters $C_0, C_1, \dots, C_k$ concurrently with zero memory synchronization or locks.
  3. Hardware Agnostic: Running counter $C_i$ on a CPU or GPU guarantees bit-for-bit identical outputs.

2. Philox Algorithm Mechanics (Philox4x32-10)

Philox is a Feistel-like Substitution-Permutation Network. It breaks a 128-bit counter state into four 32-bit words:

$$C = [R_1, L_1, R_0, L_0]$$

The Round Function Mechanics

For each round $r \in [0, 9]$ (10 rounds total):

  1. Wide Multiplication ($32 \times 32 \to 64$):
    Multiply the two lower word channels ($R_0$ and $R_1$) by fixed Philox multiplicative constants $M_0$ and $M_1$:
    $$\text{prod}_0 = R_0 \times M_0 \quad \implies \quad [\text{hi}_0, \text{lo}_0] = \text{split64}(\text{prod}_0)$$
    $$\text{prod}_1 = R_1 \times M_1 \quad \implies \quad [\text{hi}_1, \text{lo}_1] = \text{split64}(\text{prod}_1)$$
  2. Key XOR & Non-Linear Mixing:
    XOR the high 32-bits ($\text{hi}_i$) with the current round key $K_i$:
    $$X_0 = \text{hi}_0 \oplus K_0$$
    $$X_1 = \text{hi}_1 \oplus K_1$$
  3. Permutation (Word Swapping):
    Permute the lower product halves ($\text{lo}_i$), high-XOR results ($X_i$), and remaining unmultiplied words ($L_0, L_1$) to form the state for the next round:
    $$R_0' = L_0, \quad L_0' = X_0, \quad R_1' = L_1, \quad L_1' = X_1 \quad \dots \text{(Word Permutation)}$$
  4. Weyl Sequence Key Update:
    To prevent structural symmetries across rounds, update the key components using odd Weyl constants:
    $$K_0^{(r+1)} = K_0^{(r)} + W_0 \pmod{2^{32}}$$
    $$K_1^{(r+1)} = K_1^{(r)} + W_1 \pmod{2^{32}}$$

The Philox Constants

  • Multipliers:
    • $M_0 = \texttt{0xCD9E8D57}$
    • $M_1 = \texttt{0xD2511F53}$
  • Weyl Addition Constants:
    • $W_0 = \texttt{0x9E3779B9}$ ($\lfloor (\sqrt{5}-1)/2 \cdot 2^{32} \rfloor$, Golden Ratio constant)
    • $W_1 = \texttt{0xBB67AE85}$ ($\lfloor (\sqrt{3}-1)/2 \cdot 2^{32} \rfloor$)

3. Why Scalar Philox Fails on CPUs & How Vectorization Solves It

The Scalar Bottleneck (Why XGBoost saw $10\times$ slowdown)

Modern CPU execution units (ALUs) process 32-bit scalar instructions sequentially.

  1. Instruction Latency: A scalar $32 \times 32 \to 64$-bit wide multiply (mul / mulhi) requires 3–4 clock cycles of execution latency on x86 execution pipelines.
  2. Dependency Chains: Because $R^{(r+1)}$ depends directly on the result of $\text{prod}^{(r)}$, CPU instruction-level parallelism (ILP) stalls waiting for wide multiplications to complete.
  3. Auto-Vectorization Failure: Compilers cannot automatically turn a single scalar Philox sequence into vector code because the algorithm operates on scalar integer states sequentially.

SIMD Interleaving Theory (The Solution)

Instead of vectorizing the internal operations of a single Philox counter state, interleave multiple independent Philox state evaluation streams across SIMD lanes.

  • AVX2 (256-bit Vector Registers):
    • Holds 8 32-bit words (or 4 64-bit integer lanes).
    • Load 4 separate Philox counters into vector registers:
      $$\vec{C}_{\text{AVX2}} = [C_0, C_1, C_2, C_3]$$
    • Execute _mm256_mul_epu32: This single instruction multiplies the lower 32 bits of 4 separate 64-bit vector slots simultaneously, yielding four 64-bit wide-multiplies in 1 clock cycle.

By evaluating 4 counters concurrently, you amortize multiply latency across 4 streams, eliminating the sequential pipeline stalls and matching/exceeding std::mt19937 throughput.

4. IEEE-754 Fast Floating-Point Generation Mechanics

Converting an unsigned 32-bit integer $U \in [0, 2^{32}-1]$ to a normalized float $F \in [0.0, 1.0)$ traditionally uses integer-to-float conversion and floating-point division:

$$F = \frac{(\text{float})U}{2^{32}}$$
Floating-point division (vdivps) is computationally expensive (~10–14 cycles latency).

Bitwise Mantissa Injection Theory

An IEEE-754 single-precision float consists of:

$$\text{Sign (1 bit)} \mid \text{Exponent (8 bits)} \mid \text{Mantissa (23 bits)}$$

  1. Fixed Exponent: Set the sign bit to 0 and exponent bits to 127 (0x3F800000 in hex), representing $2^{127-127} = 2^0 = 1.0$.
  2. Bit Injection: Take the top 23 bits of random integer $U$ and bitwise OR them directly into the mantissa field:
    $$F_{\text{raw}} = \text{bit\_cast<float>}(\texttt{0x3F800000} \mid (U \gg 9))$$
    This constructs a uniform floating-point value in the half-open interval $[1.0, 2.0)$.
  3. Subtraction Shift: Subtract $1.0f$:
    $$F = F_{\text{raw}} - 1.0f \quad \in [0.0, 1.0)$$

This bitwise transformation executes in 1 cycle via vector OR (vpor) and vector floating-point subtract (vsubps), bypassing integer-to-float divisions entirely.

5. C++20 Concept Integration Theory

To seamlessly integrate with modern C++ standard library distribution templates (std::uniform_real_distribution, std::normal_distribution), a generator must satisfy the std::uniform_random_bit_generator concept:

$$\text{requires } G \implies \begin{cases} \text{unsigned integral type } \text{G::result\_type} \\ \text{static constexpr } \text{G::min()} \\ \text{static constexpr } \text{G::max()} \\ \text{instance call } \text{g()} \to \text{G::result\_type} \end{cases}$$
By implementing an internal ring-buffer that caches the 16 bytes generated per vectorized SIMD iteration, the C++20 wrapper provides a $O(1)$ scalar operator() interface while keeping the underlying CPU execution pipeline fed with SIMD vector blocks.