0
# Random Number Generation
1

2
High-quality parallel random number generation using cryptographically secure algorithms (Philox, Threefry) suitable for Monte Carlo simulations, stochastic computations, and statistical sampling with excellent statistical properties and performance.
3

4
## Capabilities
5

6
### Random Number Generators
7

8
Cryptographically secure parallel random number generators based on counter-based algorithms.
9

10
```python { .api }
11
class Random123GeneratorBase:
12
    """
13
    Base class for Random123 family generators (Philox, Threefry).
14
    
15
    Attributes:
16
    - context (Context): OpenCL context
17
    - dtype: Output data type
18
    """
19

20
class PhiloxGenerator(Random123GeneratorBase):
21
    """
22
    Philox random number generator with excellent statistical properties.
23
    
24
    The Philox generator uses a Feistel-like cipher structure and provides
25
    high-quality random numbers suitable for parallel Monte Carlo simulations.
26
    """
27
    
28
    def __init__(self, context, dtype=np.float32):
29
        """
30
        Create Philox random number generator.
31
        
32
        Parameters:
33
        - context (Context): OpenCL context
34
        - dtype: Output data type (float32, float64, int32, int64)
35
        """
36
    
37
    def fill_uniform(self, queue, ary, luxury=None):
38
        """
39
        Fill array with uniform random numbers in [0, 1).
40
        
41
        Parameters:
42
        - queue (CommandQueue): Command queue
43
        - ary (Array): Target array to fill
44
        - luxury (int, optional): Quality level (higher = better quality)
45
        """
46
    
47
    def uniform(self, queue, shape, dtype=None, luxury=None):
48
        """
49
        Generate uniform random array.
50
        
51
        Parameters:
52
        - queue (CommandQueue): Command queue
53
        - shape (tuple[int, ...]): Array shape
54
        - dtype: Data type (uses generator default if None)
55
        - luxury (int, optional): Quality level
56
        
57
        Returns:
58
        Array: Array of uniform random numbers in [0, 1)
59
        """
60

61
class ThreefryGenerator(Random123GeneratorBase):
62
    """
63
    Threefry random number generator based on the Threefish cipher.
64
    
65
    The Threefry generator provides high-quality random numbers with
66
    strong cryptographic properties and excellent parallel performance.
67
    """
68
    
69
    def __init__(self, context, dtype=np.float32):
70
        """
71
        Create Threefry random number generator.
72
        
73
        Parameters:
74
        - context (Context): OpenCL context
75
        - dtype: Output data type
76
        """
77
    
78
    def fill_uniform(self, queue, ary, luxury=None):
79
        """
80
        Fill array with uniform random numbers in [0, 1).
81
        
82
        Parameters:
83
        - queue (CommandQueue): Command queue
84
        - ary (Array): Target array to fill
85
        - luxury (int, optional): Quality level
86
        """
87
    
88
    def uniform(self, queue, shape, dtype=None, luxury=None):
89
        """
90
        Generate uniform random array.
91
        
92
        Parameters:
93
        - queue (CommandQueue): Command queue  
94
        - shape (tuple[int, ...]): Array shape
95
        - dtype: Data type
96
        - luxury (int, optional): Quality level
97
        
98
        Returns:
99
        Array: Array of uniform random numbers in [0, 1)
100
        """
101
```
102

103
### High-level Random Functions
104

105
Convenient functions for generating random arrays with various distributions.
106

107
```python { .api }
108
def rand(queue, shape, dtype=float, luxury=None, generator=None):
109
    """
110
    Generate uniform random array in [0, 1).
111
    
112
    Parameters:
113
    - queue (CommandQueue): Command queue
114
    - shape (int | tuple[int, ...]): Array shape
115
    - dtype: Data type (float32, float64)
116
    - luxury (int, optional): Quality level for generator
117
    - generator (Random123GeneratorBase, optional): Specific generator to use
118
    
119
    Returns:
120
    Array: Array of uniform random numbers
121
    """
122

123
def fill_rand(result, queue=None, luxury=None, generator=None):
124
    """
125
    Fill existing array with uniform random numbers.
126
    
127
    Parameters:
128
    - result (Array): Target array to fill
129
    - queue (CommandQueue, optional): Command queue
130
    - luxury (int, optional): Quality level
131
    - generator (Random123GeneratorBase, optional): Specific generator to use
132
    """
133
```
134

135
## Usage Examples
136

137
### Basic Random Number Generation
138

139
```python
140
import pyopencl as cl
141
import pyopencl.array as cl_array
142
from pyopencl.clrandom import PhiloxGenerator, ThreefryGenerator
143
import pyopencl.clrandom as clrand
144
import numpy as np
145

146
# Setup
147
ctx = cl.create_some_context()
148
queue = cl.CommandQueue(ctx)
149

150
# Generate uniform random numbers using convenience function
151
uniform_data = clrand.rand(queue, (10000,), dtype=np.float32)
152
print(f"Random data shape: {uniform_data.shape}")
153
print(f"Random samples: {uniform_data.get()[:5]}")
154
print(f"Range: [{uniform_data.get().min():.3f}, {uniform_data.get().max():.3f}]")
155

156
# Fill existing array with random data
157
existing_array = cl_array.empty(queue, (5000,), dtype=np.float32)
158
clrand.fill_rand(existing_array, queue)
159
print(f"Filled array samples: {existing_array.get()[:5]}")
160
```
161

162
### Using Specific Generators
163

164
```python
165
import pyopencl as cl
166
import pyopencl.array as cl_array
167
from pyopencl.clrandom import PhiloxGenerator, ThreefryGenerator
168
import numpy as np
169

170
# Setup
171
ctx = cl.create_some_context()
172
queue = cl.CommandQueue(ctx)
173

174
# Create Philox generator
175
philox_gen = PhiloxGenerator(ctx, dtype=np.float32)
176

177
# Generate random array
178
philox_data = philox_gen.uniform(queue, (10000,))
179
print(f"Philox generator data: {philox_data.get()[:5]}")
180

181
# Create Threefry generator
182
threefry_gen = ThreefryGenerator(ctx, dtype=np.float32)
183

184
# Generate random array
185
threefry_data = threefry_gen.uniform(queue, (10000,))
186
print(f"Threefry generator data: {threefry_data.get()[:5]}")
187

188
# Fill existing array
189
target_array = cl_array.empty(queue, (8000,), dtype=np.float32)
190
philox_gen.fill_uniform(queue, target_array)
191
print(f"Filled with Philox: {target_array.get()[:5]}")
192
```
193

194
### Monte Carlo Simulation Example
195

196
```python
197
import pyopencl as cl
198
import pyopencl.array as cl_array
199
from pyopencl.clrandom import PhiloxGenerator
200
from pyopencl.elementwise import ElementwiseKernel
201
import numpy as np
202

203
# Setup
204
ctx = cl.create_some_context()
205
queue = cl.CommandQueue(ctx)
206

207
# Monte Carlo estimation of π using random points in unit square
208
n_samples = 1000000
209

210
# Create generator
211
generator = PhiloxGenerator(ctx, dtype=np.float32)
212

213
# Generate random points
214
x_coords = generator.uniform(queue, (n_samples,))
215
y_coords = generator.uniform(queue, (n_samples,))
216

217
# Create kernel to check if points are inside unit circle
218
inside_circle_kernel = ElementwiseKernel(ctx,
219
    "__global float *x, __global float *y, __global int *inside",
220
    "inside[i] = (x[i]*x[i] + y[i]*y[i] <= 1.0f) ? 1 : 0",
221
    "check_inside_circle")
222

223
# Count points inside circle
224
inside_flags = cl_array.empty(queue, (n_samples,), dtype=np.int32)
225
inside_circle_kernel(x_coords, y_coords, inside_flags)
226

227
# Estimate π
228
points_inside = cl_array.sum(inside_flags).get()
229
pi_estimate = 4.0 * points_inside / n_samples
230

231
print(f"Samples: {n_samples}")
232
print(f"Points inside circle: {points_inside}")
233
print(f"π estimate: {pi_estimate}")
234
print(f"Error: {abs(pi_estimate - np.pi):.6f}")
235
```
236

237
### Statistical Quality Testing
238

239
```python
240
import pyopencl as cl
241
import pyopencl.array as cl_array
242
from pyopencl.clrandom import PhiloxGenerator, ThreefryGenerator
243
import numpy as np
244
import matplotlib.pyplot as plt
245

246
# Setup
247
ctx = cl.create_some_context()
248
queue = cl.CommandQueue(ctx)
249

250
# Generate large sample for statistical testing
251
n_samples = 100000
252

253
# Test both generators
254
philox_gen = PhiloxGenerator(ctx, dtype=np.float32)
255
threefry_gen = ThreefryGenerator(ctx, dtype=np.float32)
256

257
philox_data = philox_gen.uniform(queue, (n_samples,)).get()
258
threefry_data = threefry_gen.uniform(queue, (n_samples,)).get()
259

260
# Statistical tests
261
print("Statistical Analysis:")
262
print(f"Philox - Mean: {np.mean(philox_data):.6f} (expected: 0.5)")
263
print(f"Philox - Std:  {np.std(philox_data):.6f} (expected: {1/np.sqrt(12):.6f})")
264
print(f"Threefry - Mean: {np.mean(threefry_data):.6f}")
265
print(f"Threefry - Std:  {np.std(threefry_data):.6f}")
266

267
# Histogram comparison
268
bins = np.linspace(0, 1, 50)
269
philox_hist, _ = np.histogram(philox_data, bins)
270
threefry_hist, _ = np.histogram(threefry_data, bins)
271

272
print(f"Philox histogram uniformity (chi-square): {np.var(philox_hist):.2f}")
273
print(f"Threefry histogram uniformity (chi-square): {np.var(threefry_hist):.2f}")
274
```
275

276
### Generating Different Random Distributions
277

278
```python
279
import pyopencl as cl
280
import pyopencl.array as cl_array
281
from pyopencl.clrandom import PhiloxGenerator
282
from pyopencl.elementwise import ElementwiseKernel
283
import pyopencl.clmath as clmath
284
import numpy as np
285

286
# Setup
287
ctx = cl.create_some_context()
288
queue = cl.CommandQueue(ctx)
289
generator = PhiloxGenerator(ctx, dtype=np.float32)
290

291
# Generate uniform random numbers as base
292
uniform1 = generator.uniform(queue, (10000,))
293
uniform2 = generator.uniform(queue, (10000,))
294

295
# Box-Muller transform for normal distribution
296
normal_kernel = ElementwiseKernel(ctx,
297
    "__global float *u1, __global float *u2, __global float *n1, __global float *n2",
298
    """
299
    float r = sqrt(-2.0f * log(u1[i]));
300
    float theta = 2.0f * M_PI * u2[i];
301
    n1[i] = r * cos(theta);
302
    n2[i] = r * sin(theta);
303
    """,
304
    "box_muller_transform")
305

306
normal1 = cl_array.empty_like(uniform1)
307
normal2 = cl_array.empty_like(uniform2)
308
normal_kernel(uniform1, uniform2, normal1, normal2)
309

310
print(f"Normal distribution samples: {normal1.get()[:5]}")
311
print(f"Normal mean: {cl_array.sum(normal1).get() / normal1.size:.6f}")
312

313
# Exponential distribution: -log(1-u)
314
exponential = -clmath.log(1.0 - uniform1)
315
print(f"Exponential samples: {exponential.get()[:5]}")
316

317
# Integer random numbers in range [0, max_val)
318
max_val = 100
319
uniform_int_base = generator.uniform(queue, (10000,))
320
int_kernel = ElementwiseKernel(ctx,
321
    "__global float *u, __global int *result, int max_val",
322
    "result[i] = (int)(u[i] * max_val)",
323
    "uniform_int")
324

325
random_ints = cl_array.empty(queue, (10000,), dtype=np.int32)
326
int_kernel(uniform_int_base, random_ints, np.int32(max_val))
327
print(f"Random integers [0, {max_val}): {random_ints.get()[:10]}")
328
```
329

330
### Performance Comparison
331

332
```python
333
import pyopencl as cl
334
import pyopencl.array as cl_array
335
from pyopencl.clrandom import PhiloxGenerator, ThreefryGenerator
336
import numpy as np
337
import time
338

339
# Setup
340
ctx = cl.create_some_context()
341
queue = cl.CommandQueue(ctx)
342

343
# Large array for performance testing
344
n_samples = 10000000
345
n_iterations = 5
346

347
# Test Philox generator
348
philox_gen = PhiloxGenerator(ctx, dtype=np.float32)
349
philox_times = []
350

351
for i in range(n_iterations):
352
    start_time = time.time()
353
    data = philox_gen.uniform(queue, (n_samples,))
354
    queue.finish()
355
    end_time = time.time()
356
    philox_times.append(end_time - start_time)
357

358
# Test Threefry generator
359
threefry_gen = ThreefryGenerator(ctx, dtype=np.float32)
360
threefry_times = []
361

362
for i in range(n_iterations):
363
    start_time = time.time()
364
    data = threefry_gen.uniform(queue, (n_samples,))
365
    queue.finish()
366
    end_time = time.time()
367
    threefry_times.append(end_time - start_time)
368

369
# Compare with NumPy
370
numpy_times = []
371
for i in range(n_iterations):
372
    start_time = time.time()
373
    data = np.random.random(n_samples).astype(np.float32)
374
    end_time = time.time()
375
    numpy_times.append(end_time - start_time)
376

377
print(f"Performance Comparison ({n_samples:,} samples):")
378
print(f"Philox (GPU):   {np.mean(philox_times):.4f}s ± {np.std(philox_times):.4f}s")
379
print(f"Threefry (GPU): {np.mean(threefry_times):.4f}s ± {np.std(threefry_times):.4f}s")
380
print(f"NumPy (CPU):    {np.mean(numpy_times):.4f}s ± {np.std(numpy_times):.4f}s")
381

382
# Throughput
383
philox_throughput = n_samples / np.mean(philox_times) / 1e6
384
threefry_throughput = n_samples / np.mean(threefry_times) / 1e6
385
numpy_throughput = n_samples / np.mean(numpy_times) / 1e6
386

387
print(f"Throughput (Million samples/sec):")
388
print(f"Philox:   {philox_throughput:.1f}")
389
print(f"Threefry: {threefry_throughput:.1f}")
390
print(f"NumPy:    {numpy_throughput:.1f}")
391
```
392

393
## Generator Properties
394

395
### Quality Levels (Luxury)
396

397
The `luxury` parameter controls the quality vs. speed tradeoff:
398

399
- **luxury=None** (default): Standard quality, good for most applications
400
- **luxury=1**: Higher quality, slightly slower generation
401
- **luxury=2**: Highest quality, best statistical properties
402

403
Higher luxury levels provide better statistical properties but may reduce performance.
404

405
### Statistical Properties
406

407
Both Philox and Threefry generators provide:
408
- Period: 2^128 or larger
409
- Excellent statistical properties (pass BigCrush tests)
410
- Strong cryptographic properties
411
- Parallel generation without correlation
412
- Reproducible sequences with same seeds
413

414
### Memory Usage
415

416
Random number generation is memory-efficient:
417
- Generators store minimal state
418
- Memory usage scales with output array size
419
- No significant memory overhead compared to computation arrays