From 900de8dcbb829251de6001e8bf07af45eeb2cf18 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Tue, 11 Aug 2026 07:54:42 +0000
Subject: [PATCH 1/3] Add hv_approx_fpras R function with C wrapper, tests,
 vignette section, and NEWS update

Co-authored-by: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com>
---
 c/NEWS.md                          |   4 +
 c/common.h                         |  25 ++++
 c/config.h                         |   7 +
 c/gcc_attribs.h                    |  11 ++
 c/hvapprox.c                       |  92 +++++++++++++
 c/hvapprox.h                       |  12 ++
 c/rng.h                            |  84 ++++++++++++
 c/rng_alias.h                      | 210 +++++++++++++++++++++++++++++
 r/NAMESPACE                        |   1 +
 r/NEWS.md                          |   1 +
 r/R/hv_approx.R                    | 103 ++++++++++++++
 r/inst/REFERENCES.bib              |  14 ++
 r/src/Rmoocore.c                   |  17 +++
 r/src/init.h                       |   1 +
 r/tests/testthat/test-hv_approx.R  |  17 +++
 r/vignettes/articles/hv_approx.Rmd |  57 ++++++++
 16 files changed, 656 insertions(+)
 create mode 100644 c/rng_alias.h

diff --git a/c/NEWS.md b/c/NEWS.md
index 389c63de6..0a66bcb2c 100644
--- a/c/NEWS.md
+++ b/c/NEWS.md
@@ -3,6 +3,10 @@
 ## 0.20
 
  * radixsort.h, insort.h: New.
+ * rng_alias.h: New. Walker-Vose alias method.
+ * rng.h: New functions for sampling uniformly distributed 32-bits integer in a
+   bounded interval and for sampling with a given CDF (roulette-wheel method).
+ * hvapprox.c (hv_approx_fpras): New function.
 
 
 ## 0.19.2
diff --git a/c/common.h b/c/common.h
index d1e7deff4..dff5de98c 100644
--- a/c/common.h
+++ b/c/common.h
@@ -129,4 +129,29 @@ printf_point(const char * prefix, const double * p, dimension_t dim,
     fprintf(stderr, "%s", suffix);
 }
 
+static inline double
+cumsum_of_vector_double(double * restrict v, size_t n)
+{
+    ASSUME(n > 0);
+    for (size_t i = 1; i < n; i++)
+        v[i] += v[i-1];
+    return v[n - 1];
+}
+
+/** Uses Kahan summation. */
+static inline double
+kahan_sum_of_vector_double(const double * restrict v, size_t n)
+{
+    double sum = 0.0, compensation = 0.0;
+    for (size_t i = 0; i < n; i++) {
+        double y = v[i] - compensation;
+        double z = sum + y;
+        compensation = (z - sum) - y;
+        sum = z;
+    }
+    return sum;
+}
+
+
+
 #endif 	    /* !MOOCORE_COMMON_H_ */
diff --git a/c/config.h b/c/config.h
index e904dd378..3e4046cf0 100644
--- a/c/config.h
+++ b/c/config.h
@@ -28,6 +28,13 @@ typedef uint_fast8_t dimension_t;
 #define MOOCORE_HVAPPROX_DIMENSION_MAX 31
 #define MOOCORE_HV_DIMENSION_MAX 31
 
+#ifndef UINT32_MAX
+# define UINT32_MAX ((uint32_t)-1)
+#endif
+#ifndef UINT64_MAX
+# define UINT64_MAX ((uint64_t)-1)
+#endif
+
 // Use boolvec for boolean arrays that will be passed from/to R/Python. It also
 // helps with autovectorization.
 #ifdef R_PACKAGE
diff --git a/c/gcc_attribs.h b/c/gcc_attribs.h
index 9b22402e0..4f553fa14 100644
--- a/c/gcc_attribs.h
+++ b/c/gcc_attribs.h
@@ -236,4 +236,15 @@
 # define ASAN_POISON_MEMORY_REGION(addr, size) ((void) (addr), (void) (size))
 #endif
 
+/* GCC/Clang: use builtins to avoid MinGW isfinite/isnormal macro oddities and
+             -Werror=float-conversion false positives.
+*/
+#if defined(__clang__) || defined(__GNUC__)
+#  define is_finite(x) __builtin_isfinite(x)
+#  define is_normal(x) __builtin_isnormal(x)
+#else
+#  define is_finite(x) isfinite(x) // C99
+#  define is_normal(x) isnormal(x) // C99
+#endif
+
 #endif /* GCC_ATTRIBUTES */
diff --git a/c/hvapprox.c b/c/hvapprox.c
index 6a8636087..256200d4c 100644
--- a/c/hvapprox.c
+++ b/c/hvapprox.c
@@ -5,6 +5,8 @@
 #include "hvapprox.h"
 #include "pow_int.h"
 #include "rng.h"
+#include "rng_alias.h"
+#include "sort.h"
 
 #define ALMOST_ZERO_WEIGHT 1e-20
 
@@ -863,3 +865,93 @@ hv_approx_rphi_fang_wang_plus(const double * restrict data, size_t npoints, dime
     const long double c_m = sphere_area_div_2_pow_d_times_d[dim];
     return STATIC_CAST(double, c_m * (expected / STATIC_CAST(long double, nsamples)));
 }
+
+
+// FIXME: Make sure that this is vectorized.
+static double *
+compute_vols(const double * restrict points, size_t npoints, dimension_t dim)
+{
+    double * vol = malloc(npoints * sizeof(*vol));
+    for (size_t i = 0; i < npoints; i++) {
+        const double * restrict pi = points + i * dim;
+        double tmp = pi[0];
+        for (dimension_t d = 1; d < dim; d++)
+            tmp *= pi[d];
+        vol[i] = tmp;
+    }
+    return vol;
+}
+
+/**
+   FPRAS (fully polynomial-time randomized approximation scheme)
+
+   This function returns, with probability :math:`(1 - \delta)`, an
+   :math:`\epsilon`-approximation of the hypervolume metric with respect to the
+   given reference point.
+
+   delta : Probability of failure.
+
+   K. Bringmann, T. Friedrich. Approximating the volume of unions and
+   intersections of high-dimensional geometric objects. Computational Geometry:
+   Theory and Applications, Vol. 43, pages 601-610,. 2010.
+
+   This implementation uses Walker-Vose's Alias method to sample from a
+   discrete distribution in O(1). The naive roulette-wheel approach requires
+   O(log n).
+*/
+double
+hv_approx_fpras(const double * restrict data, size_t npoints, dimension_t dim,
+                const double * restrict ref, const boolvec * restrict maximise,
+                uint32_t random_seed, double epsilon, double delta)
+{
+    ASSUME(2 <= dim && dim <= MOOCORE_HVAPPROX_DIMENSION_MAX);
+    const double * points = transform_and_filter(data, &npoints, dim, ref, maximise);
+    if (points == NULL)
+        return 0;
+
+    ASSUME(0 < npoints && npoints < UINT32_MAX);
+    ASSUME(0 < epsilon && epsilon < 1);
+    ASSUME(0 < delta && delta < 1);
+
+    const double T_factor = 8 * (log(2) - log(delta)) * (1. + epsilon) / (epsilon * epsilon);
+    const double T_double = T_factor * (double)npoints;
+    if (T_double >= (double)UINT64_MAX)
+        return -1; // This will run for too long!
+
+    const uint64_t T = (uint64_t) ceil(T_double);
+    double * vols = compute_vols(points, npoints, dim); // VolumeQuery(B_i)
+    double total_vol = kahan_sum_of_vector_double(vols, npoints);
+
+    rng_state * rng = rng_new(random_seed);
+    rng_alias_sampler_t * sampler = rng_alias_sampler_new(vols, (uint32_t) npoints);
+
+    uint64_t t_sum = 0, m  = 0;
+    while (true) {
+        // Sampling should be done in O(1).
+        uint32_t i = rng_alias_sampler_choose(rng, sampler);
+
+        // SampleQuery(B_i)
+        double x[MOOCORE_DIMENSION_MAX + 1];
+        for (dimension_t d = 0; d < dim; d++)
+            x[d] = rng_random(rng);
+        const double * restrict p_i = points + i * dim;
+        for (dimension_t d = 0; d < dim; d++)
+            x[d] *= p_i[d];
+
+        while (true) {
+            if (unlikely(t_sum >= T)) {
+                free(vols);
+                rng_free(rng);
+                rng_alias_sampler_free(sampler);
+                free((void *) points);
+                return total_vol * (T_factor / (double) m);
+            }
+            t_sum++;
+            uint32_t j = rng_uniform_u32_ubound(rng, (uint32_t) npoints);
+            const double * restrict p_j = points + j * dim;
+            if (weakly_dominates(x, p_j, dim)) // PointQuery(x, B_j)
+                break;
+        }
+        m++;
+    }
+}
diff --git a/c/hvapprox.h b/c/hvapprox.h
index fd7af7d45..c8cba61f2 100644
--- a/c/hvapprox.h
+++ b/c/hvapprox.h
@@ -37,5 +37,17 @@ MOOCORE_API double hv_approx_rphi_fang_wang_plus(
     const double * restrict ref, const boolvec * restrict maximise,
     uint_fast32_t nsamples);
 
+/**
+   FPRAS (fully polynomial-time randomized approximation scheme)
+
+   K. Bringmann, T. Friedrich. Approximating the volume of unions and
+   intersections of high-dimensional geometric objects. Computational Geometry:
+   Theory and Applications, Vol. 43, pages 601-610,. 2010.
+*/
+MOOCORE_API double hv_approx_fpras(
+    const double * restrict data, size_t npoints, dimension_t dim,
+    const double * restrict ref, const boolvec * restrict maximise,
+    uint32_t random_seed, double epsilon, double delta);
+
 END_C_DECLS
 #endif // HV_APPROX_H_
diff --git a/c/rng.h b/c/rng.h
index 00876ec1b..c1fc0940c 100644
--- a/c/rng.h
+++ b/c/rng.h
@@ -1,3 +1,6 @@
+#ifndef RNG_H
+#define RNG_H
+
 #include "mt19937/mt19937.h"
 
 typedef mt19937_state rng_state;
@@ -34,9 +37,90 @@ rng_uniform(rng_state * rng, double low, double high)
     return low + (high - low) * rng_random(rng);
 }
 
+/**
+   Returns a uniformly distributed 32-bits integer in [0, n).
+
+   Lemire, Daniel. "Fast Random Integer Generation in an Interval", ACM
+   Transactions on Modeling and Computer Simulation (TOMACS), 29(1):1-12,
+   2019. https://doi.org/10.1145/323063
+
+*/
+static inline uint32_t
+rng_uniform_u32_ubound(rng_state * rng, uint32_t n)
+{
+    uint64_t m = ((uint64_t)mt19937_next32(rng)) * n;
+    uint32_t leftover = (uint32_t)m;
+    if (leftover < n) {
+        // t = 2^32 mod n, expressed without requiring a 64-bit 2^32 value.
+        const uint32_t t = (uint32_t)(-n) % n;
+        while (leftover < t) {
+            m = ((uint64_t)mt19937_next32(rng)) * n;
+            leftover = (uint32_t)m;
+        }
+    }
+    return (uint32_t)(m >> 32);
+}
+
+
+static inline uint32_t
+rng_uniform_u32_bounded(rng_state * rng, uint32_t low, uint32_t high)
+{
+    assert(rng != NULL);
+    if (low >= high)
+        return low;
+
+    return low + rng_uniform_u32_ubound(rng, high - low);
+}
+
+static inline void
+rng_validate_cdf(const double * cdf _attr_maybe_unused, uint32_t n)
+{
+    assert(n > 0);
+    assert(cdf[0] >= 0.0);
+    for (uint32_t i = 1; i < n; ++i) {
+        assert(cdf[i - 1] <= cdf[i]);
+        assert(0.0 <= cdf[i] && cdf[i] <= 1.0);
+    }
+    assert(cdf[n - 1] == 1.0);
+}
+
+/**
+   Returns a random value within [0, n - 1] based on the given CDF using the
+   roulette-wheel method.
+
+   The CDF values should be non-decreasing and within [0,1].
+
+   Requires O(log n).
+*/
+static inline uint32_t
+rng_random_wheel_uint32(rng_state * rng, const double * cdf, uint32_t n)
+{
+    DEBUG1(rng_validate_cdf(cdf, n)); // Check the CDF is correct.
+
+    double r = rng_random(rng);
+    if (r < cdf[0])
+        return 0;
+    if (n == 2 || r >= cdf[n - 2])
+        return n - 1;
+
+    // Binary search.
+    uint32_t low = 1, high = n - 2;
+    while (low < high) {
+        uint32_t mid = low  + (high - low) / 2;
+        if (r < cdf[mid])
+            high = mid;
+        else
+            low = mid + 1;
+    }
+
+    DEBUG1(for (uint32_t j = 0; j < low; j++) assert(cdf[j] < r));
+    DEBUG1(for (uint32_t j = low; j < n; j++) assert(r <= cdf[j]));
+    return low;
+}
 
 double rng_standard_normal(rng_state *rng);
 void rng_bivariate_normal_fill(rng_state * rng,
                                double mu1, double mu2,
                                double sigma1, double sigma2, double rho,
                                double *out, int n);
+#endif /* RNG_H */
diff --git a/c/rng_alias.h b/c/rng_alias.h
new file mode 100644
index 000000000..cbd31c1a3
--- /dev/null
+++ b/c/rng_alias.h
@@ -0,0 +1,210 @@
+/******************************************************************************
+
+  O(1)-per-sample discrete distribution sampler built on Walker-Vose's
+  algorithm for the alias method.
+
+  Given n outcomes with (possibly unnormalized) weights p[0..n-1], the
+  Walker-Vose algorithm builds a table in O(n) time that then supports drawing
+  a weighted-random outcome in O(1) time.
+
+  The main interface is:
+
+   rng_alias_sampler_t *
+   rng_alias_sampler_new(const double * probabilities, uint32_t n);
+
+   void rng_alias_sampler_free(rng_alias_sampler_t * sampler);
+
+   uint32_t rng_alias_sampler_choose(rng_state * rng,
+                                     const rng_alias_sampler_t * sampler);
+
+  Based on Keith Schwarz's Java reference implementation:
+
+     https://www.keithschwarz.com/interesting/code/?dir=alias-method
+
+  but with many significant changes. In particular, this implementation uses
+  two 32-bit random values.  It could be even faster if we had a 64-bits RNG.
+
+  See also:
+
+   https://en.wikipedia.org/wiki/Alias_method
+   https://jugit.fz-juelich.de/mlz/ransampl
+   https://www.keithschwarz.com/darts-dice-coins/
+
+*****************************************************************************/
+#ifndef RNG_ALIAS_SAMPLER_H
+#define RNG_ALIAS_SAMPLER_H
+
+#include <math.h>
+#include <float.h>
+#include "rng.h"
+
+/**
+   One row of the alias table. Kept as a single struct (rather than two
+   parallel arrays) so a sample only ever touches one cache line.
+*/
+typedef struct rng_alias_entry {
+    uint32_t cutoff;
+    uint32_t alias; // Outcome to return when the coin flip above fails.
+} rng_alias_entry_t;
+
+typedef struct rng_alias_sampler {
+    uint32_t n;
+    rng_alias_entry_t table[]; // Flexible array.
+} rng_alias_sampler_t;
+
+#define RNG_ALIAS_RESIDUAL_TOL (64.0 * DBL_EPSILON)
+
+static inline double
+rng_alias_clean_scaled(double q)
+{
+    if (q < 0.0 && q > -RNG_ALIAS_RESIDUAL_TOL)
+        return 0.0;
+
+    if (q < 1.0 && q > 1.0 - RNG_ALIAS_RESIDUAL_TOL)
+        return 1.0;
+
+    if (q > 1.0 && q < 1.0 + RNG_ALIAS_RESIDUAL_TOL)
+        return 1.0;
+
+    return q;
+}
+
+static inline uint32_t
+rng_alias_cutoff(double q)
+{
+    if (q <= 0.0)
+        return 0;
+
+    double x = q * 0x1p32; // 2^32
+    /* Prevent conversion of an out-of-range floating-point value to uint32_t
+      if rounding produces a value >= 2^32 - 1.  */
+    return (x >= (double)UINT32_MAX) ?  UINT32_MAX : (uint32_t)x;
+}
+
+static inline bool
+size_max_mul_overflows(size_t a, size_t b)
+{
+    return b != 0 && a > SIZE_MAX / b;
+}
+
+/**
+   Construct an alias sampler using Vose's algorithm.
+
+   The caller retains ownership of probabilities.
+
+   probabilities do not need to sum up to 1, but must be >= 0.
+*/
+static inline rng_alias_sampler_t *
+rng_alias_sampler_new(const double * probabilities, uint32_t n)
+{
+    if (probabilities == NULL || n == 0)
+        return NULL;
+
+    double sum = kahan_sum_of_vector_double(probabilities, n);
+    if (sum <= 0 || !is_finite(sum))
+        return NULL;
+
+    if (size_max_mul_overflows(n, sizeof(rng_alias_entry_t))
+        || size_max_mul_overflows(n, sizeof(double))) {
+        return NULL;
+    }
+
+    uint32_t small_end = 0, large_begin = n;
+    /* Temporary Vose double stack:
+         [0, small_end)    small entries
+         [large_begin, n)  large entries
+    */
+    uint32_t * stack = (uint32_t *) malloc(n * sizeof(*stack));
+    double * scaled = (double *) malloc(n * sizeof(*scaled));
+    rng_alias_sampler_t * sampler = malloc(sizeof(*sampler) + n * sizeof(rng_alias_entry_t));
+    if (stack == NULL || sampler == NULL || scaled == NULL) {
+        free(stack);
+        free(scaled);
+        free(sampler);
+        return NULL;
+    }
+    sampler->n = n;
+
+    for (uint32_t i = 0; i < n; i++) {
+        double p = probabilities[i];
+        if (p < 0 || !is_finite(p)) {
+            free(stack);
+            free(scaled);
+            free(sampler);
+            return NULL;
+        }
+        double q = (p / sum) * (double)n;
+        q = rng_alias_clean_scaled(q);
+        scaled[i] = q;
+
+        if (q < 1.0)
+            stack[small_end++] = i;
+        else
+            stack[--large_begin] = i;
+    }
+
+    while (small_end != 0 && large_begin != n) {
+        uint32_t small_index = stack[--small_end];
+        uint32_t large_index = stack[large_begin++];
+        double q = scaled[small_index];
+
+        sampler->table[small_index].cutoff = rng_alias_cutoff(q);
+        sampler->table[small_index].alias = large_index;
+        scaled[large_index] += q - 1.0;
+        scaled[large_index] = rng_alias_clean_scaled(scaled[large_index]);
+
+        if (scaled[large_index] < 1.0)
+            stack[small_end++] = large_index;
+        else
+            stack[--large_begin] = large_index;
+    }
+
+    while (large_begin != n) {
+        uint32_t i = stack[large_begin++];
+        sampler->table[i].cutoff = UINT32_MAX;
+        sampler->table[i].alias = i;
+    }
+
+    while (small_end != 0) {
+        uint32_t i = stack[--small_end];
+        sampler->table[i].cutoff = UINT32_MAX;
+        sampler->table[i].alias = i;
+    }
+
+    free(stack);
+    free(scaled);
+    return sampler;
+}
+
+
+static inline void
+rng_alias_sampler_free(rng_alias_sampler_t * sampler)
+{
+    free(sampler);
+}
+
+
+/**
+   Sample an integer between [0, sampler->n).
+
+   It consumes two 32-bits random values:
+
+    1. One for unbiased column selection.
+    2. One for the fixed-point alias decision.
+
+    Column selection may very rarely consume another 32-bits value when
+    rejection is required.
+
+    No floating-point operations are performed.
+*/
+static inline uint32_t
+rng_alias_sampler_choose(rng_state * rng, const rng_alias_sampler_t * sampler)
+{
+    assert(rng != NULL && sampler != NULL);
+    uint32_t column = rng_uniform_u32_ubound(rng, sampler->n);
+    uint32_t u = mt19937_next32(rng);
+    const rng_alias_entry_t * entry = sampler->table + column;
+    return u < entry->cutoff ? column : entry->alias;
+}
+
+#endif /* RNG_ALIAS_SAMPLER_H */
diff --git a/r/NAMESPACE b/r/NAMESPACE
index 8a331f2fa..798f64bd9 100644
--- a/r/NAMESPACE
+++ b/r/NAMESPACE
@@ -15,6 +15,7 @@ export(epsilon_mult)
 export(filter_dominated)
 export(generate_ndset)
 export(hv_approx)
+export(hv_approx_fpras)
 export(hv_contributions)
 export(hypervolume)
 export(igd)
diff --git a/r/NEWS.md b/r/NEWS.md
index 90caebaad..ace4405b0 100644
--- a/r/NEWS.md
+++ b/r/NEWS.md
@@ -1,6 +1,7 @@
 # moocore (development)
 
  * `is_nondominated` is up to 10x faster in some inputs thanks to a customized radixsort implementation.
+ * `hv_approx_fpras()`: New function. Approximates the hypervolume indicator via a fully polynomial-time randomized approximation scheme (FPRAS) [@BriFri2010approx].
 
 
 # moocore 0.3.2
diff --git a/r/R/hv_approx.R b/r/R/hv_approx.R
index 820cb4755..b903c16cf 100644
--- a/r/R/hv_approx.R
+++ b/r/R/hv_approx.R
@@ -120,3 +120,106 @@ hv_approx <- function(x, reference, maximise = FALSE, nsamples = 262144L, seed =
       as.integer(nsamples)))
   }
 }
+
+
+#' Approximate the hypervolume indicator via a fully polynomial-time randomized approximation scheme (FPRAS).
+#'
+#' This function implements the approximation algorithm by
+#' \citet{BriFri2010approx}.  This algorithm returns, with probability
+#' \eqn{(1 - \delta)}, an \eqn{\epsilon}-approximation of the hypervolume
+#' metric with respect to the given reference point, assuming minimization of
+#' all objectives by default.
+#'
+#' @details
+#'
+#' This function computes an approximation \eqn{\hat{v}} of the true
+#' hypervolume \eqn{v = \text{hyp}_r(A)} of the input points in \eqn{A
+#' \subset \mathbb{R}^m} with respect to the reference point \eqn{r \in
+#' \mathbb{R}^m}, such that
+#'
+#' \deqn{\text{Pr}[(1-\epsilon)v \leq \hat{v} \leq (1+\epsilon)v] \geq (1 - \delta)}
+#'
+#' where \eqn{\epsilon > 0} and \eqn{0 < \delta < 1}.
+#'
+#' The algorithm requires \eqn{O(\frac{nm}{\epsilon^2}\log\frac{1}{\delta})}.
+#' That is, it is linear on the number of points and dimensions, but quadratic
+#' in the approximation error.  In other words, more accurate approximations
+#' require significantly more time.
+#'
+#' In contrast to the (quasi)-Monte-Carlo methods provided by [hv_approx()],
+#' the presence of weakly-dominated points not only increases the runtime, but
+#' also changes the returned approximation for a fixed random seed.
+#'
+#' The implementation uses Walker-Vose's alias method for sampling from a
+#' discrete probability distribution \citep{Vose1991alias}, which requires
+#' \eqn{O(1)} per sample.  Using the naive roulette-wheel method would add, at
+#' least, a factor of \eqn{O(\log n)} to the above runtime.
+#'
+#' @inheritParams hv_approx
+#'
+#' @return A single numerical value.
+#'
+#' @param epsilon `double(1)`\cr Desired relative error of the approximation,
+#'   \eqn{\epsilon > 0}.
+#'
+#' @param delta `double(1)`\cr Desired failure probability
+#'   \eqn{0 < \delta < 1}; \eqn{(1 - \delta)} gives the confidence level.
+#'
+#' @warning Lower values of `epsilon` (\eqn{\epsilon}) or `delta`
+#'   (\eqn{\delta}) require significantly longer computation time.
+#'
+#' @seealso [hypervolume()], [whv_hype()], [hv_approx()]
+#'
+#' @author Manuel \enc{López-Ibáñez}{Lopez-Ibanez}
+#'
+#' @references
+#'
+#' \insertAllCited{}
+#'
+#' @doctest
+#'
+#' x <- matrix(c(5, 5, 4, 6, 2, 7, 7, 4), ncol=2, byrow=TRUE)
+#' @expect equal(38.0)
+#' hypervolume(x, ref=10)
+#' @expect equal(37.999979)
+#' hv_approx(x, ref=10, method="Rphi-FWE+")
+#' @expect equal(38.1446, tolerance=1e-4)
+#' hv_approx_fpras(x, ref=10, epsilon=0.1, delta=0.2, seed=42)
+#' @expect equal(37.9541, tolerance=1e-4)
+#' hv_approx_fpras(x, ref=10, epsilon=0.01, delta=0.2, seed=42)
+#'
+#' @export
+#' @concept metrics
+hv_approx_fpras <- function(x, reference, maximise = FALSE, seed = NULL,
+                            epsilon = 0.01, delta = 0.1)
+{
+  x <- as_double_matrix(x)
+  nobjs <- ncol(x)
+
+  if (!is.numeric(reference))
+    stop("a numerical reference vector must be provided")
+  if (length(reference) == 1L) reference <- rep_len(reference, nobjs)
+  stopifnot(length(reference) == nobjs)
+
+  if (length(maximise) == 1L) maximise <- rep_len(maximise, nobjs)
+  stopifnot(length(maximise) == nobjs)
+  check_dimension_max(nobjs, .libmoocore_constants[["MOOCORE_HVAPPROX_DIMENSION_MAX"]])
+
+  if (!is.numeric(epsilon) || length(epsilon) != 1L || epsilon <= 0)
+    stop("epsilon must be a positive numeric value")
+  if (!is.numeric(delta) || length(delta) != 1L || delta <= 0 || delta >= 1)
+    stop("delta must be strictly within (0, 1)")
+
+  seed <- if (is.null(seed)) get_seed() else as_integer(seed)
+  hv <- .Call(hv_approx_fpras_C,
+    t(x),
+    as.double(reference),
+    as.logical(maximise),
+    seed,
+    as.double(epsilon),
+    as.double(delta))
+  if (hv < 0)
+    stop("The requested approximation (epsilon=", epsilon, ", delta=", delta,
+         ") would require a very long time")
+  hv
+}
diff --git a/r/inst/REFERENCES.bib b/r/inst/REFERENCES.bib
index de557e2dc..ffe9e72dd 100644
--- a/r/inst/REFERENCES.bib
+++ b/r/inst/REFERENCES.bib
@@ -409,6 +409,20 @@ @article{SchKer2025r2v2
   doi = {10.1162/evco.a.366}
 }
 
+@article{Vose1991alias,
+  author = {Michael D. Vose},
+  title = {A linear algorithm for generating random numbers with a given
+                  distribution},
+  journal = {IEEE Transactions on Software Engineering},
+  year = 1991,
+  volume = 17,
+  number = 9,
+  pages = {972--975},
+  annote = {Proposed Walker-Vose alias method for sampling from a
+                  discrete distribution},
+  doi = {10.1109/32.92917}
+}
+
 @article{WuAza2001metrics,
   author = {J. Wu and S. Azam},
   title = {Metrics for Quality Assessment of a Multiobjective Design
diff --git a/r/src/Rmoocore.c b/r/src/Rmoocore.c
index 78fec90a6..3ed8ae6fb 100644
--- a/r/src/Rmoocore.c
+++ b/r/src/Rmoocore.c
@@ -512,6 +512,23 @@ hv_approx_rphi_fang_wang_plus_C(SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP N
     return Rf_ScalarReal(hv);
 }
 
+SEXP
+hv_approx_fpras_C(SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP SEED, SEXP EPSILON, SEXP DELTA)
+{
+    SEXP_2_DOUBLE_MATRIX(DATA, data, nobj, npoints);
+    SEXP_2_DOUBLE_VECTOR(REFERENCE, ref, reference_len);
+    SEXP_2_LOGICAL_INT_VECTOR(MAXIMISE, maximise, maximise_len);
+    SEXP_2_UINT32(SEED, seed);
+
+    assert(nobj == reference_len);
+    assert(nobj == maximise_len);
+
+    double epsilon = Rf_asReal(EPSILON);
+    double delta = Rf_asReal(DELTA);
+    double hv = hv_approx_fpras(data, npoints, nobj, ref, maximise, seed, epsilon, delta);
+    return Rf_ScalarReal(hv);
+}
+
 #include "r2_exact.h"
 
 SEXP
diff --git a/r/src/init.h b/r/src/init.h
index ec7a6eb8d..f9808d096 100644
--- a/r/src/init.h
+++ b/r/src/init.h
@@ -21,4 +21,5 @@ DECLARE_CALL(whv_hype_C, SEXP DATA, SEXP IDEAL, SEXP REFERENCE, SEXP NSAMPLES, S
 DECLARE_CALL(hv_approx_dz2019_mc_C, SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP NSAMPLES, SEXP SEED)
 DECLARE_CALL(hv_approx_dz2019_hw_C, SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP NSAMPLES)
 DECLARE_CALL(hv_approx_rphi_fang_wang_plus_C, SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP NSAMPLES)
+DECLARE_CALL(hv_approx_fpras_C, SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP SEED, SEXP EPSILON, SEXP DELTA)
 DECLARE_CALL_VOID(libmoocore_constants)
diff --git a/r/tests/testthat/test-hv_approx.R b/r/tests/testthat/test-hv_approx.R
index ed21e4bc3..9fc744409 100644
--- a/r/tests/testthat/test-hv_approx.R
+++ b/r/tests/testthat/test-hv_approx.R
@@ -5,6 +5,17 @@ test_that("hv_approx errors", {
   expect_equal(hv_approx(x, ref = 1), 0)
 })
 
+test_that("hv_approx_fpras errors", {
+  x <- matrix(c(0, 0), ncol = 2)
+  expect_error(hv_approx_fpras(x, reference = c(1, 1), epsilon = 0),
+    "epsilon must be a positive numeric value")
+  expect_error(hv_approx_fpras(x, reference = c(1, 1), delta = 0),
+    "delta must be strictly within")
+  expect_error(hv_approx_fpras(x, reference = c(1, 1), epsilon = 1e-9, delta = 0.001),
+    "would require a very long time")
+  expect_equal(hv_approx_fpras(x, reference = c(1, 1)), 0)
+})
+
 for (dim in seq(3L, 10L)) {
   test_that(paste0("hv_approx dim=", dim), {
     x <- matrix(replicate(dim, 0.5), ncol=dim)
@@ -27,5 +38,11 @@ for (dim in seq(3L, 10L)) {
     expect_equal(hv_approx(x, ref=ref), true_hv,
       tolerance = 10**-signif, info = paste0("dim=", dim, ", signif=", signif,
         " error=", log10(abs(true_hv - appr_hv) / true_hv)))
+
+    appr_hv <- hv_approx_fpras(x, ref=ref)
+    signif <- 4
+    expect_equal(appr_hv, true_hv,
+      tolerance = 10**-signif, info = paste0("dim=", dim, ", signif=", signif,
+        " error=", log10(abs(true_hv - appr_hv) / true_hv)))
   })
 }
diff --git a/r/vignettes/articles/hv_approx.Rmd b/r/vignettes/articles/hv_approx.Rmd
index a46c8b73c..6ca180f9b 100644
--- a/r/vignettes/articles/hv_approx.Rmd
+++ b/r/vignettes/articles/hv_approx.Rmd
@@ -158,4 +158,61 @@ ggplot(df, aes(x = samples, y = hverror, color = Method)) +
   theme_bw()
 ```
 
+# FPRAS: Fully Polynomial-Time Randomized Approximation Scheme
+
+`hv_approx_fpras()` allows obtaining an approximation with relative error
+smaller than $\epsilon$ (`epsilon`) with a given probability $1-\delta$
+(`delta`).  As the plot below shows, the approximation error is often better
+than the requested value, but the computation time increases very quickly for
+smaller `epsilon`.
+
+```{r fpras}
+#| fig.alt = "Boxplot of relative error and line plot of CPU time for FPRAS with varying epsilon and delta."
+library(ggplot2)
+library(scales)
+shape <- "convex-sphere"
+ref <- 1.1
+npoints <- 50
+dim <- 6
+nreps <- 10
+set.seed(42)
+res <- NULL
+for (r in seq_len(nreps)) {
+  z <- generate_ndset(npoints, dim, method = shape, seed = 42 + r)
+  exact <- hypervolume(z, reference = ref)
+  for (epsilon in c(0.1, 0.025, 0.01, 0.005)) {
+    for (delta in c(0.25, 0.1, 0.05)) {
+      t_start <- proc.time()[["elapsed"]]
+      hv <- hv_approx_fpras(z, reference = ref, epsilon = epsilon, delta = delta)
+      t_end <- proc.time()[["elapsed"]] - t_start
+      res <- rbind(res, data.frame(
+        r = r,
+        epsilon = epsilon,
+        delta = delta,
+        hverror = abs(1 - hv / exact),
+        time = t_end
+      ))
+    }
+  }
+}
+res[["delta"]] <- as.factor(res[["delta"]])
+res[["epsilon"]] <- as.factor(res[["epsilon"]])
+
+ggplot(res, aes(x = epsilon, y = hverror, fill = delta)) +
+  geom_boxplot() +
+  scale_y_log10(labels = label_log()) +
+  labs(y = "Relative error") +
+  ggtitle(paste0(shape, "-", npoints, "-", dim, "d")) +
+  theme_bw()
+
+ggplot(res, aes(x = epsilon, y = time, color = delta, group = delta)) +
+  stat_summary(fun = mean, geom = "line") +
+  stat_summary(fun = mean, geom = "point") +
+  scale_y_log10(labels = label_log()) +
+  scale_x_discrete() +
+  labs(y = "CPU time (s)") +
+  ggtitle(paste0(shape, "-", npoints, "-", dim, "d")) +
+  theme_bw()
+```
+
 # References

From 34367ffbf8208c557ceff015b76b38c04f45a680 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Tue, 11 Aug 2026 07:59:19 +0000
Subject: [PATCH 2/3] Fix division-by-zero guard, DOI, and spurious comma in C
 files

Co-authored-by: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com>
---
 c/hvapprox.c | 1 +
 c/hvapprox.h | 2 +-
 c/rng.h      | 2 +-
 3 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/c/hvapprox.c b/c/hvapprox.c
index 256200d4c..c388fa504 100644
--- a/c/hvapprox.c
+++ b/c/hvapprox.c
@@ -944,6 +944,7 @@ hv_approx_fpras(const double * restrict data, size_t npoints, dimension_t dim,
                 rng_free(rng);
                 rng_alias_sampler_free(sampler);
                 free((void *) points);
+                if (unlikely(m == 0)) return 0.0;
                 return total_vol * (T_factor / (double) m);
             }
             t_sum++;
diff --git a/c/hvapprox.h b/c/hvapprox.h
index c8cba61f2..cf2d73244 100644
--- a/c/hvapprox.h
+++ b/c/hvapprox.h
@@ -42,7 +42,7 @@ MOOCORE_API double hv_approx_rphi_fang_wang_plus(
 
    K. Bringmann, T. Friedrich. Approximating the volume of unions and
    intersections of high-dimensional geometric objects. Computational Geometry:
-   Theory and Applications, Vol. 43, pages 601-610,. 2010.
+   Theory and Applications, Vol. 43, pages 601-610. 2010.
 */
 MOOCORE_API double hv_approx_fpras(
     const double * restrict data, size_t npoints, dimension_t dim,
diff --git a/c/rng.h b/c/rng.h
index c1fc0940c..43629e2c3 100644
--- a/c/rng.h
+++ b/c/rng.h
@@ -42,7 +42,7 @@ rng_uniform(rng_state * rng, double low, double high)
 
    Lemire, Daniel. "Fast Random Integer Generation in an Interval", ACM
    Transactions on Modeling and Computer Simulation (TOMACS), 29(1):1-12,
-   2019. https://doi.org/10.1145/323063
+   2019. https://doi.org/10.1145/32306366
 
 */
 static inline uint32_t

From bc5d3176f48b991a00b0130d928258ed00c65263 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Tue, 11 Aug 2026 12:11:57 +0000
Subject: [PATCH 3/3] Fix test, generate hv_approx_fpras.Rd, update doctest
 file

Co-authored-by: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com>
---
 r/man/hv_approx_fpras.Rd                  | 91 +++++++++++++++++++++++
 r/tests/testthat/test-doctest-hv_approx.R | 10 +++
 r/tests/testthat/test-hv_approx.R         |  2 +-
 3 files changed, 102 insertions(+), 1 deletion(-)
 create mode 100644 r/man/hv_approx_fpras.Rd

diff --git a/r/man/hv_approx_fpras.Rd b/r/man/hv_approx_fpras.Rd
new file mode 100644
index 000000000..e1bb09520
--- /dev/null
+++ b/r/man/hv_approx_fpras.Rd
@@ -0,0 +1,91 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/hv_approx.R
+\name{hv_approx_fpras}
+\alias{hv_approx_fpras}
+\title{Approximate the hypervolume indicator via a fully polynomial-time randomized approximation scheme (FPRAS).}
+\usage{
+hv_approx_fpras(
+  x,
+  reference,
+  maximise = FALSE,
+  seed = NULL,
+  epsilon = 0.01,
+  delta = 0.1
+)
+}
+\arguments{
+\item{x}{\code{matrix()}|\code{data.frame()}\cr Matrix or data frame of numerical
+values, where each row gives the coordinates of a point.}
+
+\item{reference}{\code{numeric()}\cr Reference point as a vector of numerical
+values.}
+
+\item{maximise}{\code{logical()}\cr Whether the objectives must be maximised
+instead of minimised. Either a single logical value that applies to all
+objectives or a vector of logical values, with one value per objective.}
+
+\item{seed}{\code{integer(1)}\cr Random seed.}
+
+\item{epsilon}{\code{double(1)}\cr Desired relative error of the approximation,
+\eqn{\epsilon > 0}.}
+
+\item{delta}{\code{double(1)}\cr Desired failure probability
+\eqn{0 < \delta < 1}; \eqn{(1 - \delta)} gives the confidence level.}
+}
+\value{
+A single numerical value.
+}
+\description{
+This function implements the approximation algorithm by
+\citet{BriFri2010approx}.  This algorithm returns, with probability
+\eqn{(1 - \delta)}, an \eqn{\epsilon}-approximation of the hypervolume
+metric with respect to the given reference point, assuming minimization of
+all objectives by default.
+}
+\details{
+This function computes an approximation \eqn{\hat{v}} of the true
+hypervolume \eqn{v = \text{hyp}_r(A)} of the input points in \eqn{A
+\subset \mathbb{R}^m} with respect to the reference point \eqn{r \in
+\mathbb{R}^m}, such that
+
+\deqn{\text{Pr}[(1-\epsilon)v \leq \hat{v} \leq (1+\epsilon)v] \geq (1 - \delta)}
+
+where \eqn{\epsilon > 0} and \eqn{0 < \delta < 1}.
+
+The algorithm requires \eqn{O(\frac{nm}{\epsilon^2}\log\frac{1}{\delta})}.
+That is, it is linear on the number of points and dimensions, but quadratic
+in the approximation error.  In other words, more accurate approximations
+require significantly more time.
+
+In contrast to the (quasi)-Monte-Carlo methods provided by \code{\link[=hv_approx]{hv_approx()}},
+the presence of weakly-dominated points not only increases the runtime, but
+also changes the returned approximation for a fixed random seed.
+
+The implementation uses Walker-Vose's alias method for sampling from a
+discrete probability distribution \citep{Vose1991alias}, which requires
+\eqn{O(1)} per sample.  Using the naive roulette-wheel method would add, at
+least, a factor of \eqn{O(\log n)} to the above runtime.
+}
+\section{Warning}{
+
+Lower values of \code{epsilon} (\eqn{\epsilon}) or \code{delta}
+(\eqn{\delta}) require significantly longer computation time.
+}
+
+\examples{
+x <- matrix(c(5, 5, 4, 6, 2, 7, 7, 4), ncol=2, byrow=TRUE)
+hypervolume(x, ref=10)
+hv_approx(x, ref=10, method="Rphi-FWE+")
+hv_approx_fpras(x, ref=10, epsilon=0.1, delta=0.2, seed=42)
+hv_approx_fpras(x, ref=10, epsilon=0.01, delta=0.2, seed=42)
+}
+\references{
+\insertAllCited{}
+}
+\seealso{
+\code{\link[=hypervolume]{hypervolume()}}, \code{\link[=whv_hype]{whv_hype()}}, \code{\link[=hv_approx]{hv_approx()}}
+}
+\author{
+Manuel \enc{López-Ibáñez}{Lopez-Ibanez}
+}
+\concept{metrics}
diff --git a/r/tests/testthat/test-doctest-hv_approx.R b/r/tests/testthat/test-doctest-hv_approx.R
index bb1a3f101..4e91c2550 100644
--- a/r/tests/testthat/test-doctest-hv_approx.R
+++ b/r/tests/testthat/test-doctest-hv_approx.R
@@ -11,3 +11,13 @@ test_that("Doctest: hv_approx", {
   expect_equal(hv_approx(x, ref = 10, seed = 42, method = "DZ2019-MC"), 38.000806)
 })
 
+test_that("Doctest: hv_approx_fpras", {
+  # Created from @doctest for `hv_approx_fpras`
+  # Source file: R/hv_approx.R
+  x <- matrix(c(5, 5, 4, 6, 2, 7, 7, 4), ncol = 2, byrow = TRUE)
+  expect_equal(hypervolume(x, ref = 10), 38)
+  expect_equal(hv_approx(x, ref = 10, method = "Rphi-FWE+"), 37.999979)
+  expect_equal(hv_approx_fpras(x, ref = 10, epsilon = 0.1, delta = 0.2, seed = 42), 38.1446, tolerance = 1e-4)
+  expect_equal(hv_approx_fpras(x, ref = 10, epsilon = 0.01, delta = 0.2, seed = 42), 37.9541, tolerance = 1e-4)
+})
+
diff --git a/r/tests/testthat/test-hv_approx.R b/r/tests/testthat/test-hv_approx.R
index 9fc744409..197e8115a 100644
--- a/r/tests/testthat/test-hv_approx.R
+++ b/r/tests/testthat/test-hv_approx.R
@@ -13,7 +13,7 @@ test_that("hv_approx_fpras errors", {
     "delta must be strictly within")
   expect_error(hv_approx_fpras(x, reference = c(1, 1), epsilon = 1e-9, delta = 0.001),
     "would require a very long time")
-  expect_equal(hv_approx_fpras(x, reference = c(1, 1)), 0)
+  expect_equal(hv_approx_fpras(matrix(c(2, 2), ncol = 2), reference = c(1, 1)), 0)
 })
 
 for (dim in seq(3L, 10L)) {