> On Sat, Apr 11, 2026 at 06:35:18PM +1200, Thomas Munro wrote: > On Wed, Jul 30, 2025 at 10:15 PM Dmitry Dolgov <[email protected]> wrote: > > As a side note, I was trying to experiment with this patch using > > dm-mapper's delay feature to introduce an arbitrary large io latency and > > see how the io queue is growing. > > FWIW, here's what I came up with while experimenting with that sort of thing: > > shared_preload_libraries=io_limit > io_limit.ios_per_second=6000 > > That differs from eg dm-mapper delays by making everything seem like > slow direct I/O, which seemed more interesting for this project. For > example if you run some continuous workload while you SET > io_limit.ios_per_second to various numbers, with > io_workers_idle_timeout set fairly low, you can monitor the pool > adjustments.
Yeah, sounds like a good idea. Do you plan to introduce such an extension long term for testing, or is it just one off? As to me it looks worth keeping, maybe even use injections points to allow for more flexibility. And I know I sound like a broken record, but if I understand correctly the delays introduced via ios_per_seconds and others are constant in time -- I've experimented a bit and found some reference implementations in numpy for geometric distribution sampling, which allow to make the delay a random variable. Since the geometric distribution is a discrete analog of the exponential one, and the latter represents delays between events in Poisson distribution, such random variable would give an approximation for more real load.
>From 2ffba516edba39e7119e27c50e1c03296e860da9 Mon Sep 17 00:00:00 2001 From: Thomas Munro <[email protected]> Date: Sat, 11 Apr 2026 17:31:13 +1200 Subject: [PATCH v1 1/2] contrib/io_limit: Simulation of slow storage. Only affects IOs submitted to io_method=worker. Configured as: shared_preload_libraries=io_limit io_limit.ios_per_second=1000 io_limit.read_per_second=200MB io_limit_write_per_second=100MB Zero means no limit. XXX Experimental hack --- contrib/Makefile | 1 + contrib/io_limit/Makefile | 20 ++ contrib/io_limit/io_limit.c | 275 ++++++++++++++++++++++++ contrib/io_limit/io_limit.control | 5 + contrib/io_limit/meson.build | 28 +++ contrib/meson.build | 1 + src/backend/storage/aio/method_worker.c | 13 ++ src/include/storage/io_worker.h | 5 + 8 files changed, 348 insertions(+) create mode 100644 contrib/io_limit/Makefile create mode 100644 contrib/io_limit/io_limit.c create mode 100644 contrib/io_limit/io_limit.control create mode 100644 contrib/io_limit/meson.build diff --git a/contrib/Makefile b/contrib/Makefile index 7d91fe77db3..48e82c53333 100644 --- a/contrib/Makefile +++ b/contrib/Makefile @@ -24,6 +24,7 @@ SUBDIRS = \ hstore \ intagg \ intarray \ + io_limit \ isn \ lo \ ltree \ diff --git a/contrib/io_limit/Makefile b/contrib/io_limit/Makefile new file mode 100644 index 00000000000..da176698a17 --- /dev/null +++ b/contrib/io_limit/Makefile @@ -0,0 +1,20 @@ +# contrib/io_limit/Makefile + +MODULE_big = io_limit +OBJS = \ + $(WIN32RES) \ + io_limit.o + +EXTENSION = io_limit +PGFILEDESC = "io_limit - io_limit - artificially limit asynchronous I/O for tesing" + +ifdef USE_PGXS +PG_CONFIG = pg_config +PGXS := $(shell $(PG_CONFIG) --pgxs) +include $(PGXS) +else +subdir = contrib/pg_prewarm +top_builddir = ../.. +include $(top_builddir)/src/Makefile.global +include $(top_srcdir)/contrib/contrib-global.mk +endif diff --git a/contrib/io_limit/io_limit.c b/contrib/io_limit/io_limit.c new file mode 100644 index 00000000000..fa2ec6f1ff2 --- /dev/null +++ b/contrib/io_limit/io_limit.c @@ -0,0 +1,275 @@ +#include "postgres.h" + +#include "miscadmin.h" +#include "port/atomics.h" +#include "portability/instr_time.h" +#include "storage/aio_internal.h" +#include "storage/io_worker.h" +#include "storage/lwlock.h" +#include "storage/shmem.h" +#include "utils/guc.h" + +/* GUCs. */ +static int io_limit_ios_per_second = 0; +static int io_limit_read_per_second = 0; +static int io_limit_write_per_second = 0; + +typedef struct io_limit_control_data +{ + /* Whether any GUC is set to a non-zero value. */ + bool enabled; + + /* Absolute time to wait until. */ + pg_atomic_uint64 op_next_ns; + pg_atomic_uint64 read_next_ns; + pg_atomic_uint64 write_next_ns; + + /* Limits expressed as delay intervals. */ + LWLock lock; + int op_ns; + int read_block_ns; + int write_block_ns; +} io_limit_control_data; + +static io_limit_control_data * io_limit_control; + +static void io_limit_shmem_request(void *arg); +static void io_limit_shmem_init(void *arg); + +static void assign_io_limit_ios_per_second(int newval, void *extra); +static void assign_io_limit_read_per_second(int newval, void *extra); +static void assign_io_limit_write_per_second(int newval, void *extra); +static const char *show_io_limit_ios_per_second(void); +static const char *show_io_limit_read_per_second(void); +static const char *show_io_limit_write_per_second(void); + +static void io_limit_on_perform(PgAioHandle *ioh); + +static const ShmemCallbacks io_limit_shmem_callbacks = { + .request_fn = io_limit_shmem_request, + .init_fn = io_limit_shmem_init, +}; + +PG_MODULE_MAGIC_EXT( + .name = "io_limit", + .version = PG_VERSION +); + +void +_PG_init(void) +{ + /* Bail out if not configured in shared_preload_libraries. */ + if (!process_shared_preload_libraries_in_progress) + return; + + DefineCustomIntVariable("io_limit.ios_per_second", + "Limits IOs per second.", + "If set to zero, there is no limit.", + &io_limit_ios_per_second, + 0, + 0, INT_MAX, + PGC_USERSET, + 0, + NULL, + assign_io_limit_ios_per_second, + show_io_limit_ios_per_second); + DefineCustomIntVariable("io_limit.read_per_second", + "Limits read bandwidth.", + "If set to zero, there is no limit.", + &io_limit_read_per_second, + 0, + 0, INT_MAX, + PGC_USERSET, + GUC_UNIT_BLOCKS, + NULL, + assign_io_limit_read_per_second, + show_io_limit_read_per_second); + DefineCustomIntVariable("io_limit.write_per_second", + "Limits write bandwidth.", + "If set to zero, there is no limit.", + &io_limit_write_per_second, + 0, + 0, INT_MAX, + PGC_USERSET, + GUC_UNIT_BLOCKS, + NULL, + assign_io_limit_write_per_second, + show_io_limit_write_per_second); + + MarkGUCPrefixReserved("io_limit"); + RegisterShmemCallbacks(&io_limit_shmem_callbacks); + pgaio_worker_set_on_perform_hook(io_limit_on_perform); +} + +static void +io_limit_shmem_request(void *arg) +{ + ShmemRequestStruct(.name = "io_limit", + .size = sizeof(io_limit_control_data), + .ptr = (void **) &io_limit_control); +} + +static void +io_limit_shmem_init(void *arg) +{ + memset(io_limit_control, 0, sizeof(*io_limit_control)); + pg_atomic_init_u64(&io_limit_control->op_next_ns, 0); + pg_atomic_init_u64(&io_limit_control->read_next_ns, 0); + pg_atomic_init_u64(&io_limit_control->write_next_ns, 0); + LWLockInitialize(&io_limit_control->lock, LWLockNewTrancheId("io_limit")); + + /* Assign initial values. */ + assign_io_limit_ios_per_second(io_limit_ios_per_second, NULL); + assign_io_limit_read_per_second(io_limit_read_per_second, NULL); + assign_io_limit_write_per_second(io_limit_write_per_second, NULL); +} + +static void +assign_io_limit(int *wait_ns, int per_second) +{ + /* Ignore call from _PG_init() before ready. */ + if (!io_limit_control) + return; + + LWLockAcquire(&io_limit_control->lock, LW_EXCLUSIVE); + *wait_ns = per_second == 0 ? 0 : NS_PER_S / per_second; + io_limit_control->enabled = + io_limit_control->op_ns > 0 || + io_limit_control->read_block_ns > 0 || + io_limit_control->write_block_ns > 0; + LWLockRelease(&io_limit_control->lock); +} + +static void +assign_io_limit_ios_per_second(int newval, void *extra) +{ + assign_io_limit(&io_limit_control->op_ns, newval); +} + +static void +assign_io_limit_read_per_second(int newval, void *extra) +{ + assign_io_limit(&io_limit_control->read_block_ns, newval); +} + +static void +assign_io_limit_write_per_second(int newval, void *extra) +{ + assign_io_limit(&io_limit_control->write_block_ns, newval); +} + +static const char * +show_io_limit(const int *wait_ns) +{ + int per_second; + + LWLockAcquire(&io_limit_control->lock, LW_SHARED); + per_second = *wait_ns == 0 ? 0 : NS_PER_S / *wait_ns; + LWLockRelease(&io_limit_control->lock); + + return psprintf("%d", per_second); +} + +static const char * +show_io_limit_ios_per_second(void) +{ + return show_io_limit(&io_limit_control->op_ns); +} + +static const char * +show_io_limit_read_per_second(void) +{ + return show_io_limit(&io_limit_control->read_block_ns); +} + +static const char * +show_io_limit_write_per_second(void) +{ + return show_io_limit(&io_limit_control->write_block_ns); +} + +static BlockNumber +io_limit_get_block_count(PgAioHandle *ioh) +{ + if (ioh->op == PGAIO_OP_READV || + ioh->op == PGAIO_OP_WRITEV) + { + struct iovec *iov; + size_t size; + int iovcnt; + + size = 0; + iovcnt = pgaio_io_get_iovec_length(ioh, &iov); + for (int i = 0; i < iovcnt; ++i) + size += iov[i].iov_len; + + return size / BLCKSZ; + } + + return 0; +} + +/* + * Wait until *next_ns_p and advance *next_ns_p by delay_ns. + */ +static void +io_limit_wait(pg_atomic_uint64 *next_ns_p, int delay_ns) +{ + instr_time now; + uint64 now_ns; + uint64 next_ns; + + INSTR_TIME_SET_CURRENT(now); + now_ns = INSTR_TIME_GET_NANOSEC(now); + next_ns = pg_atomic_read_u64(next_ns_p); + + for (;;) + { + if (next_ns > now_ns) + { + /* Need to wait. Delay the next op further. */ + next_ns = pg_atomic_fetch_add_u64(next_ns_p, delay_ns); + + /* Average rate maintained even with low-res sleep or EINTR. */ + pg_usleep(((next_ns - now_ns) + 999) / 1000); + break; + } + else + { + /* Don't need to wait. New next_ns is relative to now. */ + if (pg_atomic_compare_exchange_u64(next_ns_p, + &next_ns, + now_ns + delay_ns)) + break; + } + } +} + +static void +io_limit_on_perform(PgAioHandle *ioh) +{ + int op_ns; + int read_block_ns; + int write_block_ns; + + if (!io_limit_control->enabled) + return; + + op_ns = io_limit_control->op_ns; + if (op_ns) + io_limit_wait(&io_limit_control->op_next_ns, op_ns); + + if (ioh->op == PGAIO_OP_READV) + { + read_block_ns = io_limit_control->read_block_ns; + if (read_block_ns) + io_limit_wait(&io_limit_control->read_next_ns, + io_limit_get_block_count(ioh) * read_block_ns); + } + else if (ioh->op == PGAIO_OP_WRITEV) + { + write_block_ns = io_limit_control->write_block_ns; + io_limit_wait(&io_limit_control->write_next_ns, + io_limit_get_block_count(ioh) * write_block_ns); + } +} diff --git a/contrib/io_limit/io_limit.control b/contrib/io_limit/io_limit.control new file mode 100644 index 00000000000..2f8f06c9e87 --- /dev/null +++ b/contrib/io_limit/io_limit.control @@ -0,0 +1,5 @@ +# io_limit extension +comment = 'io_limit' +default_version = '1.0' +module_pathname = '$libdir/io_limit' +relocatable = true diff --git a/contrib/io_limit/meson.build b/contrib/io_limit/meson.build new file mode 100644 index 00000000000..1d26a08de83 --- /dev/null +++ b/contrib/io_limit/meson.build @@ -0,0 +1,28 @@ +# Copyright (c) 2022-2026, PostgreSQL Global Development Group + +io_limit_sources = files( + 'io_limit.c', +) + +if host_system == 'windows' + io_limit_sources += rc_lib_gen.process(win32ver_rc, extra_args: [ + '--NAME', 'io_limit', + '--FILEDESC', 'io_limit - artificially limit asynchronous I/O for tesing',]) +endif + +io_limit = shared_module('io_limit', + io_limit_sources, + kwargs: contrib_mod_args, +) +contrib_targets += io_limit + +install_data( + 'io_limit.control', + kwargs: contrib_data_args, +) + +tests += { + 'name': 'io_limit', + 'sd': meson.current_source_dir(), + 'bd': meson.current_build_dir(), +} diff --git a/contrib/meson.build b/contrib/meson.build index ebb7f83d8c5..398b0d704b5 100644 --- a/contrib/meson.build +++ b/contrib/meson.build @@ -34,6 +34,7 @@ subdir('hstore_plperl') subdir('hstore_plpython') subdir('intagg') subdir('intarray') +subdir('io_limit') subdir('isn') subdir('jsonb_plperl') subdir('jsonb_plpython') diff --git a/src/backend/storage/aio/method_worker.c b/src/backend/storage/aio/method_worker.c index a5ccd506d8c..87afcf856e1 100644 --- a/src/backend/storage/aio/method_worker.c +++ b/src/backend/storage/aio/method_worker.c @@ -139,6 +139,7 @@ static int MyIoWorkerId = -1; static PgAioWorkerSubmissionQueue *io_worker_submission_queue; static PgAioWorkerControl *io_worker_control; +static io_worker_on_perform_fn io_worker_on_perform_hook; static void pgaio_workerset_initialize(PgAioWorkerSet *set) @@ -529,6 +530,9 @@ pgaio_worker_submit(uint16 num_staged_ios, PgAioHandle **staged_ios) for (int i = 0; i < nsync; ++i) { pgaio_io_perform_synchronously(synchronous_ios[i]); + + if (io_worker_on_perform_hook) + io_worker_on_perform_hook(synchronous_ios[i]); } } @@ -929,6 +933,9 @@ IoWorkerMain(const void *startup_data, size_t startup_data_len) */ pgaio_io_perform_synchronously(ioh); + if (io_worker_on_perform_hook) + io_worker_on_perform_hook(ioh); + RESUME_INTERRUPTS(); errcallback.arg = NULL; } @@ -1024,6 +1031,12 @@ IoWorkerMain(const void *startup_data, size_t startup_data_len) proc_exit(0); } +void +pgaio_worker_set_on_perform_hook(io_worker_on_perform_fn fn) +{ + io_worker_on_perform_hook = fn; +} + bool pgaio_workers_enabled(void) { diff --git a/src/include/storage/io_worker.h b/src/include/storage/io_worker.h index c852c9f3741..c9ef49a585d 100644 --- a/src/include/storage/io_worker.h +++ b/src/include/storage/io_worker.h @@ -28,4 +28,9 @@ extern bool pgaio_worker_pm_test_grow_signal_sent(void); extern void pgaio_worker_pm_clear_grow_signal_sent(void); extern bool pgaio_worker_pm_test_grow(void); +/* Hook to support contrib/io_limit. */ +typedef void (*io_worker_on_perform_fn) (PgAioHandle *handle); +extern void pgaio_worker_set_on_perform_hook(io_worker_on_perform_fn fn); + + #endif /* IO_WORKER_H */ base-commit: 80156cee06b9d257251d72379ac43f9b88bd13e1 -- 2.52.0
>From d2bbf06936c7a519d31555d845d964b00ee63503 Mon Sep 17 00:00:00 2001 From: Dmitrii Dolgov <[email protected]> Date: Mon, 13 Apr 2026 20:49:11 +0200 Subject: [PATCH v1 2/2] Add delays from geometric distribution Introduce delays between IOs following a geometric distribution (a discrete analog of exponential distribution) using ziggurat algorithm. --- contrib/io_limit/io_limit.c | 47 +++- contrib/io_limit/ziggurat.h | 457 ++++++++++++++++++++++++++++++++++++ 2 files changed, 503 insertions(+), 1 deletion(-) create mode 100644 contrib/io_limit/ziggurat.h diff --git a/contrib/io_limit/io_limit.c b/contrib/io_limit/io_limit.c index fa2ec6f1ff2..7e36af03138 100644 --- a/contrib/io_limit/io_limit.c +++ b/contrib/io_limit/io_limit.c @@ -8,11 +8,13 @@ #include "storage/lwlock.h" #include "storage/shmem.h" #include "utils/guc.h" +#include "ziggurat.h" /* GUCs. */ static int io_limit_ios_per_second = 0; static int io_limit_read_per_second = 0; static int io_limit_write_per_second = 0; +static int io_limit_mean_delay = 0; typedef struct io_limit_control_data { @@ -29,6 +31,8 @@ typedef struct io_limit_control_data int op_ns; int read_block_ns; int write_block_ns; + + double ios_probability; } io_limit_control_data; static io_limit_control_data * io_limit_control; @@ -39,6 +43,7 @@ static void io_limit_shmem_init(void *arg); static void assign_io_limit_ios_per_second(int newval, void *extra); static void assign_io_limit_read_per_second(int newval, void *extra); static void assign_io_limit_write_per_second(int newval, void *extra); +static void assign_io_mean(int newval, void *extra); static const char *show_io_limit_ios_per_second(void); static const char *show_io_limit_read_per_second(void); static const char *show_io_limit_write_per_second(void); @@ -95,6 +100,18 @@ _PG_init(void) NULL, assign_io_limit_write_per_second, show_io_limit_write_per_second); + DefineCustomIntVariable("io_limit.mean_delay", + "For sampling delays from the geometric distribution, " + "mean delay in ns (1/p).", + "If set to zero, no random sampling is performed.", + &io_limit_mean_delay, + 0, + 0, INT_MAX, + PGC_USERSET, + 0, + NULL, + assign_io_mean, + show_io_limit_ios_per_second); MarkGUCPrefixReserved("io_limit"); RegisterShmemCallbacks(&io_limit_shmem_callbacks); @@ -122,6 +139,8 @@ io_limit_shmem_init(void *arg) assign_io_limit_ios_per_second(io_limit_ios_per_second, NULL); assign_io_limit_read_per_second(io_limit_read_per_second, NULL); assign_io_limit_write_per_second(io_limit_write_per_second, NULL); + + assign_io_mean(io_limit_mean_delay, NULL); } static void @@ -136,7 +155,25 @@ assign_io_limit(int *wait_ns, int per_second) io_limit_control->enabled = io_limit_control->op_ns > 0 || io_limit_control->read_block_ns > 0 || - io_limit_control->write_block_ns > 0; + io_limit_control->write_block_ns > 0 || + io_limit_control->ios_probability > 0; + LWLockRelease(&io_limit_control->lock); +} + +static void +assign_io_mean(int newval, void *extra) +{ + /* Ignore call from _PG_init() before ready. */ + if (!io_limit_control) + return; + + LWLockAcquire(&io_limit_control->lock, LW_EXCLUSIVE); + io_limit_control->ios_probability = 1.0 / (double) newval; + io_limit_control->enabled = + io_limit_control->op_ns > 0 || + io_limit_control->read_block_ns > 0 || + io_limit_control->write_block_ns > 0 || + io_limit_control->ios_probability > 0; LWLockRelease(&io_limit_control->lock); } @@ -251,10 +288,18 @@ io_limit_on_perform(PgAioHandle *ioh) int op_ns; int read_block_ns; int write_block_ns; + double ios_probability; if (!io_limit_control->enabled) return; + ios_probability = io_limit_control->ios_probability; + if (ios_probability) + { + int delay = random_geometric_inversion(io_limit_control->ios_probability); + io_limit_wait(&io_limit_control->op_next_ns, delay); + } + op_ns = io_limit_control->op_ns; if (op_ns) io_limit_wait(&io_limit_control->op_next_ns, op_ns); diff --git a/contrib/io_limit/ziggurat.h b/contrib/io_limit/ziggurat.h new file mode 100644 index 00000000000..61331d11676 --- /dev/null +++ b/contrib/io_limit/ziggurat.h @@ -0,0 +1,457 @@ +/*------------------------------------------------------------------------- + * + * ziggurat.h + * Constants and functionality for ziggurat method of generating random + * variables from: + * + * Marsaglia, George, and Wai Wan Tsang. "The ziggurat method for + * generating random variables." Journal of statistical software 5 + * (2000):1-7. + * + * Reference implementation is taken from numpy: + * numpy/random/src/distributions/ziggurat_constants.h + * + * Copyright (c) 2016-2026, PostgreSQL Global Development Group + * + * IDENTIFICATION + * contrib/io_limit/ziggurat.h + * + *------------------------------------------------------------------------- + */ +#ifndef _ZIGGURAT_H_ +#define _ZIGGURAT_H_ + +#include <math.h> + +double std_uniform(); +double random_standard_exponential(); + +static const uint64_t ke_double[] = { + 0x001C5214272497C6, 0x0000000000000000, 0x00137D5BD79C317E, + 0x00186EF58E3F3C10, 0x001A9BB7320EB0AE, 0x001BD127F719447C, + 0x001C951D0F88651A, 0x001D1BFE2D5C3972, 0x001D7E5BD56B18B2, + 0x001DC934DD172C70, 0x001E0409DFAC9DC8, 0x001E337B71D47836, + 0x001E5A8B177CB7A2, 0x001E7B42096F046C, 0x001E970DAF08AE3E, + 0x001EAEF5B14EF09E, 0x001EC3BD07B46556, 0x001ED5F6F08799CE, + 0x001EE614AE6E5688, 0x001EF46ECA361CD0, 0x001F014B76DDD4A4, + 0x001F0CE313A796B6, 0x001F176369F1F77A, 0x001F20F20C452570, + 0x001F29AE1951A874, 0x001F31B18FB95532, 0x001F39125157C106, + 0x001F3FE2EB6E694C, 0x001F463332D788FA, 0x001F4C10BF1D3A0E, + 0x001F51874C5C3322, 0x001F56A109C3ECC0, 0x001F5B66D9099996, + 0x001F5FE08210D08C, 0x001F6414DD445772, 0x001F6809F6859678, + 0x001F6BC52A2B02E6, 0x001F6F4B3D32E4F4, 0x001F72A07190F13A, + 0x001F75C8974D09D6, 0x001F78C71B045CC0, 0x001F7B9F12413FF4, + 0x001F7E5346079F8A, 0x001F80E63BE21138, 0x001F835A3DAD9162, + 0x001F85B16056B912, 0x001F87ED89B24262, 0x001F8A10759374FA, + 0x001F8C1BBA3D39AC, 0x001F8E10CC45D04A, 0x001F8FF102013E16, + 0x001F91BD968358E0, 0x001F9377AC47AFD8, 0x001F95204F8B64DA, + 0x001F96B878633892, 0x001F98410C968892, 0x001F99BAE146BA80, + 0x001F9B26BC697F00, 0x001F9C85561B717A, 0x001F9DD759CFD802, + 0x001F9F1D6761A1CE, 0x001FA058140936C0, 0x001FA187EB3A3338, + 0x001FA2AD6F6BC4FC, 0x001FA3C91ACE0682, 0x001FA4DB5FEE6AA2, + 0x001FA5E4AA4D097C, 0x001FA6E55EE46782, 0x001FA7DDDCA51EC4, + 0x001FA8CE7CE6A874, 0x001FA9B793CE5FEE, 0x001FAA9970ADB858, + 0x001FAB745E588232, 0x001FAC48A3740584, 0x001FAD1682BF9FE8, + 0x001FADDE3B5782C0, 0x001FAEA008F21D6C, 0x001FAF5C2418B07E, + 0x001FB012C25B7A12, 0x001FB0C41681DFF4, 0x001FB17050B6F1FA, + 0x001FB2179EB2963A, 0x001FB2BA2BDFA84A, 0x001FB358217F4E18, + 0x001FB3F1A6C9BE0C, 0x001FB486E10CACD6, 0x001FB517F3C793FC, + 0x001FB5A500C5FDAA, 0x001FB62E2837FE58, 0x001FB6B388C9010A, + 0x001FB7353FB50798, 0x001FB7B368DC7DA8, 0x001FB82E1ED6BA08, + 0x001FB8A57B0347F6, 0x001FB919959A0F74, 0x001FB98A85BA7204, + 0x001FB9F861796F26, 0x001FBA633DEEE286, 0x001FBACB2F41EC16, + 0x001FBB3048B49144, 0x001FBB929CAEA4E2, 0x001FBBF23CC8029E, + 0x001FBC4F39D22994, 0x001FBCA9A3E140D4, 0x001FBD018A548F9E, + 0x001FBD56FBDE729C, 0x001FBDAA068BD66A, 0x001FBDFAB7CB3F40, + 0x001FBE491C7364DE, 0x001FBE9540C9695E, 0x001FBEDF3086B128, + 0x001FBF26F6DE6174, 0x001FBF6C9E828AE2, 0x001FBFB031A904C4, + 0x001FBFF1BA0FFDB0, 0x001FC03141024588, 0x001FC06ECF5B54B2, + 0x001FC0AA6D8B1426, 0x001FC0E42399698A, 0x001FC11BF9298A64, + 0x001FC151F57D1942, 0x001FC1861F770F4A, 0x001FC1B87D9E74B4, + 0x001FC1E91620EA42, 0x001FC217EED505DE, 0x001FC2450D3C83FE, + 0x001FC27076864FC2, 0x001FC29A2F90630E, 0x001FC2C23CE98046, + 0x001FC2E8A2D2C6B4, 0x001FC30D654122EC, 0x001FC33087DE9C0E, + 0x001FC3520E0B7EC6, 0x001FC371FADF66F8, 0x001FC390512A2886, + 0x001FC3AD137497FA, 0x001FC3C844013348, 0x001FC3E1E4CCAB40, + 0x001FC3F9F78E4DA8, 0x001FC4107DB85060, 0x001FC4257877FD68, + 0x001FC438E8B5BFC6, 0x001FC44ACF15112A, 0x001FC45B2BF447E8, + 0x001FC469FF6C4504, 0x001FC477495001B2, 0x001FC483092BFBB8, + 0x001FC48D3E457FF6, 0x001FC495E799D21A, 0x001FC49D03DD30B0, + 0x001FC4A29179B432, 0x001FC4A68E8E07FC, 0x001FC4A8F8EBFB8C, + 0x001FC4A9CE16EA9E, 0x001FC4A90B41FA34, 0x001FC4A6AD4E28A0, + 0x001FC4A2B0C82E74, 0x001FC49D11E62DE2, 0x001FC495CC852DF4, + 0x001FC48CDC265EC0, 0x001FC4823BEC237A, 0x001FC475E696DEE6, + 0x001FC467D6817E82, 0x001FC458059DC036, 0x001FC4466D702E20, + 0x001FC433070BCB98, 0x001FC41DCB0D6E0E, 0x001FC406B196BBF6, + 0x001FC3EDB248CB62, 0x001FC3D2C43E593C, 0x001FC3B5DE0591B4, + 0x001FC396F599614C, 0x001FC376005A4592, 0x001FC352F3069370, + 0x001FC32DC1B22818, 0x001FC3065FBD7888, 0x001FC2DCBFCBF262, + 0x001FC2B0D3B99F9E, 0x001FC2828C8FFCF0, 0x001FC251DA79F164, + 0x001FC21EACB6D39E, 0x001FC1E8F18C6756, 0x001FC1B09637BB3C, + 0x001FC17586DCCD10, 0x001FC137AE74D6B6, 0x001FC0F6F6BB2414, + 0x001FC0B348184DA4, 0x001FC06C898BAFF0, 0x001FC022A092F364, + 0x001FBFD5710F72B8, 0x001FBF84DD29488E, 0x001FBF30C52FC60A, + 0x001FBED907770CC6, 0x001FBE7D80327DDA, 0x001FBE1E094BA614, + 0x001FBDBA7A354408, 0x001FBD52A7B9F826, 0x001FBCE663C6201A, + 0x001FBC757D2C4DE4, 0x001FBBFFBF63B7AA, 0x001FBB84F23FE6A2, + 0x001FBB04D9A0D18C, 0x001FBA7F351A70AC, 0x001FB9F3BF92B618, + 0x001FB9622ED4ABFC, 0x001FB8CA33174A16, 0x001FB82B76765B54, + 0x001FB7859C5B895C, 0x001FB6D840D55594, 0x001FB622F7D96942, + 0x001FB5654C6F37E0, 0x001FB49EBFBF69D2, 0x001FB3CEC803E746, + 0x001FB2F4CF539C3E, 0x001FB21032442852, 0x001FB1203E5A9604, + 0x001FB0243042E1C2, 0x001FAF1B31C479A6, 0x001FAE045767E104, + 0x001FACDE9DBF2D72, 0x001FABA8E640060A, 0x001FAA61F399FF28, + 0x001FA908656F66A2, 0x001FA79AB3508D3C, 0x001FA61726D1F214, + 0x001FA47BD48BEA00, 0x001FA2C693C5C094, 0x001FA0F4F47DF314, + 0x001F9F04336BBE0A, 0x001F9CF12B79F9BC, 0x001F9AB84415ABC4, + 0x001F98555B782FB8, 0x001F95C3ABD03F78, 0x001F92FDA9CEF1F2, + 0x001F8FFCDA9AE41C, 0x001F8CB99E7385F8, 0x001F892AEC479606, + 0x001F8545F904DB8E, 0x001F80FDC336039A, 0x001F7C427839E926, + 0x001F7700A3582ACC, 0x001F71200F1A241C, 0x001F6A8234B7352A, + 0x001F630000A8E266, 0x001F5A66904FE3C4, 0x001F50724ECE1172, + 0x001F44C7665C6FDA, 0x001F36E5A38A59A2, 0x001F26143450340A, + 0x001F113E047B0414, 0x001EF6AEFA57CBE6, 0x001ED38CA188151E, + 0x001EA2A61E122DB0, 0x001E5961C78B267C, 0x001DDDF62BAC0BB0, +0x001CDB4DD9E4E8C0}; + +static const double we_double[] = { + 9.655740063209182975e-16, 7.089014243955414331e-18, + 1.163941249669122378e-17, 1.524391512353216015e-17, + 1.833284885723743916e-17, 2.108965109464486630e-17, + 2.361128077843138196e-17, 2.595595772310893952e-17, + 2.816173554197752338e-17, 3.025504130321382330e-17, + 3.225508254836375280e-17, 3.417632340185027033e-17, + 3.602996978734452488e-17, 3.782490776869649048e-17, + 3.956832198097553231e-17, 4.126611778175946428e-17, + 4.292321808442525631e-17, 4.454377743282371417e-17, + 4.613133981483185932e-17, 4.768895725264635940e-17, + 4.921928043727962847e-17, 5.072462904503147014e-17, + 5.220704702792671737e-17, 5.366834661718192181e-17, + 5.511014372835094717e-17, 5.653388673239667134e-17, + 5.794088004852766616e-17, 5.933230365208943081e-17, + 6.070922932847179572e-17, 6.207263431163193485e-17, + 6.342341280303076511e-17, 6.476238575956142121e-17, + 6.609030925769405241e-17, 6.740788167872722244e-17, + 6.871574991183812442e-17, 7.001451473403929616e-17, + 7.130473549660643409e-17, 7.258693422414648352e-17, + 7.386159921381791997e-17, 7.512918820723728089e-17, + 7.639013119550825792e-17, 7.764483290797848102e-17, + 7.889367502729790548e-17, 8.013701816675454434e-17, + 8.137520364041762206e-17, 8.260855505210038174e-17, + 8.383737972539139383e-17, 8.506196999385323132e-17, + 8.628260436784112996e-17, 8.749954859216182511e-17, + 8.871305660690252281e-17, 8.992337142215357066e-17, + 9.113072591597909173e-17, 9.233534356381788123e-17, + 9.353743910649128938e-17, 9.473721916312949566e-17, + 9.593488279457997317e-17, 9.713062202221521206e-17, + 9.832462230649511362e-17, 9.951706298915071878e-17, + 1.007081177024294931e-16, 1.018979547484694078e-16, + 1.030867374515421954e-16, 1.042746244856188556e-16, + 1.054617701794576406e-16, 1.066483248011914702e-16, + 1.078344348241948498e-16, 1.090202431758350473e-16, + 1.102058894705578110e-16, 1.113915102286197502e-16, + 1.125772390816567488e-16, 1.137632069661684705e-16, + 1.149495423059009298e-16, 1.161363711840218308e-16, + 1.173238175059045788e-16, 1.185120031532669434e-16, + 1.197010481303465158e-16, 1.208910707027385520e-16, + 1.220821875294706151e-16, 1.232745137888415193e-16, + 1.244681632985112523e-16, 1.256632486302898513e-16, + 1.268598812200397542e-16, 1.280581714730749379e-16, + 1.292582288654119552e-16, 1.304601620412028847e-16, + 1.316640789066572582e-16, 1.328700867207380889e-16, + 1.340782921828999433e-16, 1.352888015181175458e-16, + 1.365017205594397770e-16, 1.377171548282880964e-16, + 1.389352096127063919e-16, 1.401559900437571538e-16, + 1.413796011702485188e-16, 1.426061480319665444e-16, + 1.438357357315790180e-16, 1.450684695053687684e-16, + 1.463044547929475721e-16, 1.475437973060951633e-16, + 1.487866030968626066e-16, 1.500329786250736949e-16, + 1.512830308253539427e-16, 1.525368671738125550e-16, + 1.537945957544996933e-16, 1.550563253257577148e-16, + 1.563221653865837505e-16, 1.575922262431176140e-16, + 1.588666190753684151e-16, 1.601454560042916733e-16, + 1.614288501593278662e-16, 1.627169157465130500e-16, + 1.640097681172717950e-16, 1.653075238380036909e-16, + 1.666103007605742067e-16, 1.679182180938228863e-16, + 1.692313964762022267e-16, 1.705499580496629830e-16, + 1.718740265349031656e-16, 1.732037273081008369e-16, + 1.745391874792533975e-16, 1.758805359722491379e-16, + 1.772279036068006489e-16, 1.785814231823732619e-16, + 1.799412295642463721e-16, 1.813074597718501559e-16, + 1.826802530695252266e-16, 1.840597510598587828e-16, + 1.854460977797569461e-16, 1.868394397994192684e-16, + 1.882399263243892051e-16, 1.896477093008616722e-16, + 1.910629435244376536e-16, 1.924857867525243818e-16, + 1.939163998205899420e-16, 1.953549467624909132e-16, + 1.968015949351037382e-16, 1.982565151475019047e-16, + 1.997198817949342081e-16, 2.011918729978734671e-16, + 2.026726707464198289e-16, 2.041624610503588774e-16, + 2.056614340951917875e-16, 2.071697844044737034e-16, + 2.086877110088159721e-16, 2.102154176219292789e-16, + 2.117531128241075913e-16, 2.133010102535779087e-16, + 2.148593288061663316e-16, 2.164282928437604723e-16, + 2.180081324120784027e-16, 2.195990834682870728e-16, + 2.212013881190495942e-16, 2.228152948696180545e-16, + 2.244410588846308588e-16, 2.260789422613173739e-16, + 2.277292143158621037e-16, 2.293921518837311354e-16, + 2.310680396348213318e-16, 2.327571704043534613e-16, + 2.344598455404957859e-16, 2.361763752697773994e-16, + 2.379070790814276700e-16, 2.396522861318623520e-16, + 2.414123356706293277e-16, 2.431875774892255956e-16, + 2.449783723943070217e-16, 2.467850927069288738e-16, + 2.486081227895851719e-16, 2.504478596029557040e-16, + 2.523047132944217013e-16, 2.541791078205812227e-16, + 2.560714816061770759e-16, 2.579822882420530896e-16, + 2.599119972249746917e-16, 2.618610947423924219e-16, + 2.638300845054942823e-16, 2.658194886341845120e-16, + 2.678298485979525166e-16, 2.698617262169488933e-16, + 2.719157047279818500e-16, 2.739923899205814823e-16, + 2.760924113487617126e-16, 2.782164236246436081e-16, + 2.803651078006983464e-16, 2.825391728480253184e-16, + 2.847393572388174091e-16, 2.869664306419817679e-16, + 2.892211957417995598e-16, 2.915044901905293183e-16, + 2.938171887070028633e-16, 2.961602053345465687e-16, + 2.985344958730045276e-16, 3.009410605012618141e-16, + 3.033809466085003416e-16, 3.058552518544860874e-16, + 3.083651274815310004e-16, 3.109117819034266344e-16, + 3.134964845996663118e-16, 3.161205703467105734e-16, + 3.187854438219713117e-16, 3.214925846206797361e-16, + 3.242435527309451638e-16, 3.270399945182240440e-16, + 3.298836492772283149e-16, 3.327763564171671408e-16, + 3.357200633553244075e-16, 3.387168342045505162e-16, + 3.417688593525636996e-16, 3.448784660453423890e-16, + 3.480481301037442286e-16, 3.512804889222979418e-16, + 3.545783559224791863e-16, 3.579447366604276541e-16, + 3.613828468219060593e-16, 3.648961323764542545e-16, + 3.684882922095621322e-16, 3.721633036080207290e-16, + 3.759254510416256532e-16, 3.797793587668874387e-16, + 3.837300278789213687e-16, 3.877828785607895292e-16, + 3.919437984311428867e-16, 3.962191980786774996e-16, + 4.006160751056541688e-16, 4.051420882956573177e-16, + 4.098056438903062509e-16, 4.146159964290904582e-16, + 4.195833672073398926e-16, 4.247190841824385048e-16, + 4.300357481667470702e-16, 4.355474314693952008e-16, + 4.412699169036069903e-16, 4.472209874259932285e-16, + 4.534207798565834480e-16, 4.598922204905932469e-16, + 4.666615664711475780e-16, 4.737590853262492027e-16, + 4.812199172829237933e-16, 4.890851827392209900e-16, + 4.974034236191939753e-16, 5.062325072144159699e-16, + 5.156421828878082953e-16, 5.257175802022274839e-16, + 5.365640977112021618e-16, 5.483144034258703912e-16, + 5.611387454675159622e-16, 5.752606481503331688e-16, + 5.909817641652102998e-16, 6.087231416180907671e-16, + 6.290979034877557049e-16, 6.530492053564040799e-16, + 6.821393079028928626e-16, 7.192444966089361564e-16, +7.706095350032096755e-16, 8.545517038584027421e-16}; + +static const double fe_double[] = { + 1.000000000000000000e+00, 9.381436808621747003e-01, + 9.004699299257464817e-01, 8.717043323812035949e-01, + 8.477855006239896074e-01, 8.269932966430503241e-01, + 8.084216515230083777e-01, 7.915276369724956185e-01, + 7.759568520401155522e-01, 7.614633888498962833e-01, + 7.478686219851951034e-01, 7.350380924314234843e-01, + 7.228676595935720206e-01, 7.112747608050760117e-01, + 7.001926550827881623e-01, 6.895664961170779872e-01, + 6.793505722647653622e-01, 6.695063167319247333e-01, + 6.600008410789997004e-01, 6.508058334145710999e-01, + 6.418967164272660897e-01, 6.332519942143660652e-01, + 6.248527387036659775e-01, 6.166821809152076561e-01, + 6.087253820796220127e-01, 6.009689663652322267e-01, + 5.934009016917334289e-01, 5.860103184772680329e-01, + 5.787873586028450257e-01, 5.717230486648258170e-01, + 5.648091929124001709e-01, 5.580382822625874484e-01, + 5.514034165406412891e-01, 5.448982376724396115e-01, + 5.385168720028619127e-01, 5.322538802630433219e-01, + 5.261042139836197284e-01, 5.200631773682335979e-01, + 5.141263938147485613e-01, 5.082897764106428795e-01, + 5.025495018413477233e-01, 4.969019872415495476e-01, + 4.913438695940325340e-01, 4.858719873418849144e-01, + 4.804833639304542103e-01, 4.751751930373773747e-01, + 4.699448252839599771e-01, 4.647897562504261781e-01, + 4.597076156421376902e-01, 4.546961574746155033e-01, + 4.497532511627549967e-01, 4.448768734145485126e-01, + 4.400651008423538957e-01, 4.353161032156365740e-01, + 4.306281372884588343e-01, 4.259995411430343437e-01, + 4.214287289976165751e-01, 4.169141864330028757e-01, + 4.124544659971611793e-01, 4.080481831520323954e-01, + 4.036940125305302773e-01, 3.993906844752310725e-01, + 3.951369818332901573e-01, 3.909317369847971069e-01, + 3.867738290841376547e-01, 3.826621814960098344e-01, + 3.785957594095807899e-01, 3.745735676159021588e-01, + 3.705946484351460013e-01, 3.666580797815141568e-01, + 3.627629733548177748e-01, 3.589084729487497794e-01, + 3.550937528667874599e-01, 3.513180164374833381e-01, + 3.475804946216369817e-01, 3.438804447045024082e-01, + 3.402171490667800224e-01, 3.365899140286776059e-01, + 3.329980687618089852e-01, 3.294409642641363267e-01, + 3.259179723935561879e-01, 3.224284849560891675e-01, + 3.189719128449572394e-01, 3.155476852271289490e-01, + 3.121552487741795501e-01, 3.087940669345601852e-01, + 3.054636192445902565e-01, 3.021634006756935276e-01, + 2.988929210155817917e-01, 2.956517042812611962e-01, + 2.924392881618925744e-01, 2.892552234896777485e-01, + 2.860990737370768255e-01, 2.829704145387807457e-01, + 2.798688332369729248e-01, 2.767939284485173568e-01, + 2.737453096528029706e-01, 2.707225967990600224e-01, + 2.677254199320447947e-01, 2.647534188350622042e-01, + 2.618062426893629779e-01, 2.588835497490162285e-01, + 2.559850070304153791e-01, 2.531102900156294577e-01, + 2.502590823688622956e-01, 2.474310756653276266e-01, + 2.446259691318921070e-01, 2.418434693988772144e-01, + 2.390832902624491774e-01, 2.363451524570596429e-01, + 2.336287834374333461e-01, 2.309339171696274118e-01, + 2.282602939307167011e-01, 2.256076601166840667e-01, + 2.229757680581201940e-01, 2.203643758433594946e-01, + 2.177732471487005272e-01, 2.152021510753786837e-01, + 2.126508619929782795e-01, 2.101191593889882581e-01, + 2.076068277242220372e-01, 2.051136562938377095e-01, + 2.026394390937090173e-01, 2.001839746919112650e-01, + 1.977470661050988732e-01, 1.953285206795632167e-01, + 1.929281499767713515e-01, 1.905457696631953912e-01, + 1.881811994042543179e-01, 1.858342627621971110e-01, + 1.835047870977674633e-01, 1.811926034754962889e-01, + 1.788975465724783054e-01, 1.766194545904948843e-01, + 1.743581691713534942e-01, 1.721135353153200598e-01, + 1.698854013025276610e-01, 1.676736186172501919e-01, + 1.654780418749360049e-01, 1.632985287519018169e-01, + 1.611349399175920349e-01, 1.589871389693142123e-01, + 1.568549923693652315e-01, 1.547383693844680830e-01, + 1.526371420274428570e-01, 1.505511850010398944e-01, + 1.484803756438667910e-01, 1.464245938783449441e-01, + 1.443837221606347754e-01, 1.423576454324722018e-01, + 1.403462510748624548e-01, 1.383494288635802039e-01, + 1.363670709264288572e-01, 1.343990717022136294e-01, + 1.324453279013875218e-01, 1.305057384683307731e-01, + 1.285802045452281717e-01, 1.266686294375106714e-01, + 1.247709185808309612e-01, 1.228869795095451356e-01, + 1.210167218266748335e-01, 1.191600571753276827e-01, + 1.173168992115555670e-01, 1.154871635786335338e-01, + 1.136707678827443141e-01, 1.118676316700562973e-01, + 1.100776764051853845e-01, 1.083008254510337970e-01, + 1.065370040500016602e-01, 1.047861393065701724e-01, + 1.030481601712577161e-01, 1.013229974259536315e-01, + 9.961058367063713170e-02, 9.791085331149219917e-02, + 9.622374255043279756e-02, 9.454918937605585882e-02, + 9.288713355604354127e-02, 9.123751663104015530e-02, + 8.960028191003285847e-02, 8.797537446727021759e-02, + 8.636274114075691288e-02, 8.476233053236811865e-02, + 8.317409300963238272e-02, 8.159798070923741931e-02, + 8.003394754231990538e-02, 7.848194920160642130e-02, + 7.694194317048050347e-02, 7.541388873405840965e-02, + 7.389774699236474620e-02, 7.239348087570873780e-02, + 7.090105516237182881e-02, 6.942043649872875477e-02, + 6.795159342193660135e-02, 6.649449638533977414e-02, + 6.504911778675374900e-02, 6.361543199980733421e-02, + 6.219341540854099459e-02, 6.078304644547963265e-02, + 5.938430563342026597e-02, 5.799717563120065922e-02, + 5.662164128374287675e-02, 5.525768967669703741e-02, + 5.390531019604608703e-02, 5.256449459307169225e-02, + 5.123523705512628146e-02, 4.991753428270637172e-02, + 4.861138557337949667e-02, 4.731679291318154762e-02, + 4.603376107617516977e-02, 4.476229773294328196e-02, + 4.350241356888818328e-02, 4.225412241331623353e-02, + 4.101744138041481941e-02, 3.979239102337412542e-02, + 3.857899550307485742e-02, 3.737728277295936097e-02, + 3.618728478193142251e-02, 3.500903769739741045e-02, + 3.384258215087432992e-02, 3.268796350895953468e-02, + 3.154523217289360859e-02, 3.041444391046660423e-02, + 2.929566022463739317e-02, 2.818894876397863569e-02, + 2.709438378095579969e-02, 2.601204664513421735e-02, + 2.494202641973178314e-02, 2.388442051155817078e-02, + 2.283933540638524023e-02, 2.180688750428358066e-02, + 2.078720407257811723e-02, 1.978042433800974303e-02, + 1.878670074469603046e-02, 1.780620041091136169e-02, + 1.683910682603994777e-02, 1.588562183997316302e-02, + 1.494596801169114850e-02, 1.402039140318193759e-02, + 1.310916493125499106e-02, 1.221259242625538123e-02, + 1.133101359783459695e-02, 1.046481018102997894e-02, + 9.614413642502209895e-03, 8.780314985808975251e-03, + 7.963077438017040002e-03, 7.163353183634983863e-03, + 6.381905937319179087e-03, 5.619642207205483020e-03, + 4.877655983542392333e-03, 4.157295120833795314e-03, + 3.460264777836904049e-03, 2.788798793574076128e-03, + 2.145967743718906265e-03, 1.536299780301572356e-03, +9.672692823271745359e-04, 4.541343538414967652e-04}; + +static const double ziggurat_exp_r = 7.6971174701310497140446280481; + +/* Returns a standard uniform random variable */ +double +std_uniform() +{ + return (double) rand() / (double) ((unsigned) RAND_MAX); +} + +static uint64_t +rand_uint64(void) +{ + uint64_t r = 0; + + for (int i = 0; i < 64; i += 15) + { + r = r * ((uint64_t) RAND_MAX + 1) + rand(); + } + + return r; +} + +static double +standard_exponential_unlikely(uint8_t idx, double x) +{ + if (idx == 0) + { + /* Switch to 1.0 - U to avoid log(0.0) */ + return ziggurat_exp_r - log1p(-std_uniform()); + } + else if ((fe_double[idx - 1] - fe_double[idx]) * std_uniform() + fe_double[idx] < exp(-x)) + { + return x; + } + else + { + return random_standard_exponential(); + } +} + +double +random_standard_exponential() +{ + uint64_t ri; + uint8_t idx; + double x; + + ri = rand_uint64(); + ri >>= 3; + idx = ri & 0xFF; + ri >>= 8; + x = ri * we_double[idx]; + + if (ri < ke_double[idx]) + { + /* 98.9% of the time we return here 1st try */ + return x; + } + + return standard_exponential_unlikely(idx, x); +} + +/* Only for probabilities less than 0.33(3) */ +static long int +random_geometric_inversion(double p) +{ + double z = ceil(-random_standard_exponential() / log1p(-p)); + + /* + * The constant 9.223372036854776e+18 is the smallest double that is + * larger than INT64_MAX. + */ + if (z >= 9.223372036854776e+18) + { + return INT64_MAX; + } + + return (int64_t) z; +} + +#endif -- 2.52.0
