Detection and Resolution of Periodic Artifacts in OpenDP's Discrete Laplace Sampler

Detection and Resolution of Periodic Artifacts in OpenDP’s Discrete Laplace Sampler Cesare Gerolimetto Fabrello∗† , Valeria Rossi∗† , Alberto Trombetta∗ , and Massimo Caccia∗† ∗ Università degli Studi dell’Insubria

arXiv:2610.01907v1 [cs.CR] 1 Oct 2026

† Random Power Srl

Abstract—Differential privacy implementations rely on precise sampling from noise distributions to provide formal privacy guarantees. We report the discovery of systematic artifacts in OpenDP’s discrete Laplace sampler that manifest as periodic distortions in the output distribution. Through systematic testing, we trace these artifacts to a faulty implementation in the rational arithmetic library used by the bernoulli_exp1 function, a low-level primitive that implements sampling from Bernoulli(exp(−x)) distributions. We present a diagnostic methodology that isolates the faulty component in the nested sampling hierarchy and propose an alternative implementation based on exact rational arithmetic that eliminates the artifacts. Statistical validation with 106 samples confirms that the corrected sampler produces outputs indistinguishable from the theoretical distribution at the tested precision level. Index Terms—Differential Privacy, Discrete Laplace, Random Sampling, Statistical Testing, OpenDP

I. I NTRODUCTION Differential privacy has emerged as the gold standard for privacy-preserving data analysis [1]. Central to any differential privacy implementation is the noise addition mechanism, which relies on sampling from carefully calibrated probability distributions. The correctness of these sampling procedures is critical. Even small deviations can compromise privacy guarantees or degrade utility. While the mathematical foundations are well understood, their realization on finite-precision hardware introduces subtle challenges. For example, Mironov [2] demonstrated that floating-point arithmetic can leak information through the least significant bits. Our work extends this to exact sampling implementations, where numerical precision issues manifest as systematic artifacts. OpenDP [3] is a widely used open-source library implementing the discrete Laplace mechanism following Canonne, Kamath, and Steinke (CKS20) [4]. This construction provides provably correct samplers through nested primitives with formal correctness guarantees. However, implementation details can introduce unexpected behaviors. In this work, we report systematic artifacts observed in OpenDP’s discrete Laplace sampler while empirically investigating the role of randomness in Differential Privacy. These artifacts manifest as periodic distortions in the empirical distribution that persist across different privacy parameters and entropy sources, suggesting a fundamental implementation issue rather than statistical fluctuation or input randomness quality. Our contributions: We present a systematic methodology for diagnosing sampling artifacts in hierarchical implemen-

tations, successfully isolating the root cause to numerical precision issues in a low-level Bernoulli sampler. We propose an alternative implementation based on exact rational arithmetic that eliminates the artifacts while maintaining the theoretical guarantees of the original algorithm. Our approach demonstrates the importance of rigorous statistical testing in validating differential privacy implementations. II. O PEN DP L IBRARY OpenDP [5] is an open-source differential privacy library providing implementations of standard mechanisms. The discrete Laplace sampler follows the construction of CKS20 [4], decomposing the sampling procedure into a hierarchy of primitives with formal correctness proofs. This modular approach enables verification of individual components while providing a clear separation between mathematical guarantees and implementation details. III. A RTIFACT D ISCOVERY A. Experimental Setup All experiments were conducted on a high-performance server with dual Intel Xeon Gold 5416S processors, providing 64 logical cores (32 physical cores with hyper-threading) operating between 800 MHz and 4.0 GHz. The system has 128 GB of RAM organized in a NUMA architecture with two nodes for efficient memory allocation. These resources were required by our experimental demands: millions of iterations for each combination of parameters and randomness sources. Samples were gathered

Fig. 1: Discrete Laplace distribution generated by OpenDP’s sampler showing the characteristic periodic artifacts. The insets on the plot show a zoomed region around the artifacts, with error bars representing 1σ confidence intervals.

(a) Uniformity fit on the last digit of sampled values.

(b) sample_uniform_ubig_below uniformity fit with upper bound 10.

Fig. 2: Uniformity tests on base sampling primitives. Left: last digit distribution of sampled uniform values (χ2 = 9.04, dof = 9, p = 0.433). Right: sample_uniform_ubig_below with upper bound 10 (χ2 = 6.83, dof = 9, p = 0.654). Both tests show χ2red ≈ 1 and p > 0.4, confirming the entropy source operates correctly.

(a) Standard Bernoulli: z = 0.58, p = 0.56.

(b) Rational Bernoulli: z = 0.66, p = 0.51.

Fig. 3: Bernoulli samplers that draw directly from the entropy source. Left: standard Bernoulli with p = 0.5 (p̂ = 0.5003, pth = 0.5000, χ2 = 0.339, p = 0.561). Right: rational Bernoulli with p = 0.1 (p̂ = 0.1002, pth = 0.1000, χ2 = 0.431, p = 0.511). Both pass goodnessof-fit tests with χ2red < 1 and p > 0.5. Insets show 1σ confidence intervals.

using OpenDP v.0.14.0, with all post-processing and statistical analysis performed in MATLAB. Codebase is available on a dedicated Github repository1 .

of which entropy source feeds the library, including highquality quantum random number generators and OpenDP’s native source. This strongly suggests the issue originates in the internal sampling logic rather than in the external randomness.

B. Empirical Distribution Reconstruction We collected N = 106 samples from OpenDP’s discrete Laplace sampler to reconstruct empirical distributions. Our query counts the number of individuals in a dataset, which has sensitivity ∆s = 1. Throughout this analysis, we focus on ε = 0.1, which yields a Laplace scale of b = 1/ε = 10, but the results are reproducible with different parameter choices. Figure 1 shows the reconstructed histogram using a bin size of 1. The samples exhibit evenly spaced artifacts that suggest a systematic implementation issue rather than statistical fluctuation. The choice of unitary bin size is critical: larger bins smooth over these artifacts and obscure the underlying problem. The large sample size also contributes to highlight the discrepancies. These artifacts exhibit clear symmetry with respect to the distribution center and a periodic pattern whose spacing depends on the Laplace scale. Crucially, they appear regardless 1 https://github.com/grlcsr/dp analysis

IV. I MPLEMENTATION S TRUCTURE According to OpenDP’s documentation [5], the implementation follows the construction described by Canonne, Kamath, and Steinke [4]. The sampling logic is decomposed into nested primitives that call each other in a tree-like hierarchy: • sample_discrete_laplace: top-level sampler for Lap(0, scale), which draws a random sign via sam ple_standard_bernoulli and a magnitude via sample_geometric_exp_fast • sample_geometric_exp_fast: geometric sampler with parameter exp(−x), which repeatedly calls sampl e_uniform_ubig_below and sample_bernoull i_exp • sample_geometric_exp_slow: reference geometric sampler defined as a sequence of independent Bernoulli(exp(−x)) trials • sample_standard_bernoulli: Bernoulli(0.5) implemented by extracting a single random bit

(a) bernoulli_exp1: z = 11.3, p ≈ 10−29 .

(b) bernoulli_exp: z = 11.3, p ≈ 10−29 .

Fig. 4: Exponential Bernoulli samplers showing significant deviations from the theoretical model. Left: bernoulli_exp (p̂ = 0.6696, pth = 0.6643, χ2 = 128, p = 0). Right: bernoulli_exp1 (p̂ = 0.6696, pth = 0.6643, χ2 = 128, p = 0). Both fail goodness-of-fit tests with χ2 ≈ 128 and p ≈ 0. Insets zoom on the 1σ confidence intervals.

(a) Geometric slow: good fit to the theoretical model.

(b) Geometric fast: fails with visible periodic artifacts.

Fig. 5: Geometric sampler distributions. Left: slow implementation (mean z = 0.96, p = 0.34; χ2 = 15, p = 0.26) passes with good fit to the theoretical model. Right: fast implementation (mean z = 3.99, p = 6.6 × 10−5 ; χ2 = 1.9 × 103 , p = 0) fails with visible periodic artifacts similar to those observed in the final Laplace distribution. Insets zoom on the 1σ confidence intervals.

sample_bernoulli_exp: Bernoulli(exp(−x)) for general x ≥ 0, implemented by repeatedly invoking sam ple_bernoulli_exp1 on the integer and fractional parts of x • sample_bernoulli_exp1: Bernoulli(exp(−x)) restricted to x ∈ [0, 1], based on an alternating Taylor series • sample_bernoulli_rational: Bernoulli(p) for rational p = a/b • sample_uniform_ubig_below: rejection sampler for uniform big integer on {0, . . . , upper − 1} based on repeated calls to the underlying byte-level RNG.

•

V. D IAGNOSTIC M ETHODOLOGY To localize the source of the artifacts, we adopted a topdown diagnostic approach. We implemented a traced version of OpenDP’s CKS20 sampler that records representative outputs from every primitive in the call chain for each discrete Laplace sample generated. Specifically, a single call to the traced sample_dis crete_laplace returns both the final discrete Laplace sample and a diagnostic row containing: the first uniform value drawn; the first value from sample_uniform_ubig_b elow; the first outcomes of sample_bernoulli_stan

dard, sample_bernoulli_rational, sample_bern oulli_exp1, and sample_bernoulli_exp; the first magnitudes from sample_geometric_exp_slow and sample_geometric_exp_fast; and the final discrete Laplace sample. By iterating N = 106 times, we saved these cascades of samples, all generated from the same entropy source. We then performed goodness-of-fit tests (χ2 and z-tests) on each column against the corresponding theoretical distribution. Because all columns derive from the same random bytes, any deviation must originate from one of the deterministic transformations. The key insight is that as soon as one primitive distorts the distribution, the discrepancy becomes visible in its column and propagates to all downstream columns, allowing us to pinpoint the earliest faulty step in the sampling pipeline. VI. L AYER - BY-L AYER A NALYSIS A. Uniform Sampling We began by verifying the uniformity of the bit-stream and the rejection sampler function sample_uniform_ubig_b elow. Results are presented in Figure 2. Both distributions yield χ2red ≈ 1 and p-values of p = 0.433 (uniform bit sampling) and p = 0.654 (thresholded uniform sampler) when

tested against the theoretical uniform model. These results confirm that the entropy source and uniform integer sampling layer operate correctly. B. Bernoulli Samplers The next layer consists of the Bernoulli samplers: samp le_standard_bernoulli, sample_bernoulli_ex p, sample_bernoulli_exp1, and sample_bernoull i_rational. The standard Bernoulli is called only by the Laplace sampler to determine the sign and draws directly from the uniform entropy source. The other three call each other sequentially in a nested environment, with only the rational sampler accessing the entropy source directly. To test conformity, we performed both a χ2 goodnessof-fit test and a z-test on the empirical success probability p̂ = k1 /n against the theoretical value pth , using the normal approximation to the binomial: z=p

p̂ − pth pth (1 − pth )/n

b−1

p = E[e

1 X −u/b 1 1 − e−1 ]= e = · , b u=0 b 1 − e−1/b

where U is uniformly sampled from {0, . . . , b − 1}. The sample_standard_bernoulli and sample _bernoulli_rational functions pass both tests. The standard Bernoulli returns z = 0.58 with p = 0.56, and the rational Bernoulli returns z = 0.66 with p = 0.51. In both cases, the z-value sits within 1σ of the expected value, with correspondingly high p-values (Figure 3). The exponential functions, however, fail the tests dramatically. The bernoulli_exp1 function returns z = 11.30 with p ≈ 0 and χ2 = 127.7. The bernoulli_exp function shows identical statistics, which is expected given the calling hierarchy. Since function calls proceed from bernoulli_e xp to bernoulli_exp1 to bernoulli_rational, the failing function is most likely bernoulli_exp1, with the error propagating upward to bernoulli_exp (Figure 4). C. Geometric Samplers Moving to the geometric samplers, we observe a clear divergence. The slow geometric distribution appears correct, while the fast geometric sampler exhibits artifacts similar to those in the final Laplacian output and fails statistical tests. We performed both z-tests and χ2 tests against the geometric distribution P (K = k) = (1 − e−x ) e−xk ,

Draw U ∼ Uniform({0, . . . , 2128 − 1}); u ← U/2128 4 Sprev ← 1; Scurr ← 1 − x; k ← 1; t ← x 5 while True do 6 L ← min(Sprev , Scurr ); H ← max(Sprev , Scurr ) 7 if u < L then 8 return True 3

9 10 11 12 13

The z-value represents the number of standard deviations the observed probability deviates from the expected value; values within |z| < 3 indicate agreement with the theoretical model at 99% confidence. The standard Bernoulli was tested against p = 0.5. The rational Bernoulli was tested against p = 1/b = 0.1 (matching the Laplace scale). Both exponential functions were tested against −U/b

Algorithm 1: Sample Bernoulli(exp(−x)) via Taylor Series Input: x ∈ [0, 1] Output: True with probability exp(−x) 1 if x = 0 then 2 return True

k = 0, 1, 2, . . .

14 15

if u ≥ H then return False k ← k + 1; t ← t · x/k; Sprev ← Scurr if k is even then Scurr ← Scurr + t else Scurr ← Scurr − t

where x = 1 for sample_geometric_exp_slow and x = 1/b for sample_geometric_exp_fast. The slow geometric sampler appears correct despite depending on the faulty bernoulli_exp because it is always called with parameter x = 1, meaning it only uses a single Bernoulli with p = e−1 , implemented via bernoulli_exp 1(1). As discussed in Section IX, the bug in the underlying dashu library produces errors that depend on the specific values involved. For this particular input the resulting bias is negligible, and the function returns a correct result. The fast geometric sampler, in contrast, repeatedly samples from an entire family of slightly biased coins with probabilities e−u/b . This over-selects some samples while under-selecting others, distorting the entire geometric distribution and causing test failures (Figure 5). D. Final Laplace Distribution The generated discrete Laplace distribution (Figure 1) exhibits the artifacts we first observed (mean z = −0.68 with p = 0.49; χ2 = 2.1 × 103 with p = 0). While the mean is well-centered, the distribution fails goodness-of-fit testing due to systematic deviations. The top-down analysis has successfully traced the problem from its visible manifestation in the final output back to its origin in the bernoulli_exp1 function. In the following Section, we propose an alternative implementation of this sampler that eliminates the artifacts while leaving the rest of the CKS20 construction untouched. VII. A LTERNATIVE S AMPLING A LGORITHM Our alternative implementation is based on the proof of Proposition 33 in CKS20 [4] and on an alternating series Bernoulli sampler described in [6]. The algorithm samples from a Bernoulli(exp(−x)) distribution by iteratively refining

original implementation, which produced z ≈ 11 and p ≈ 0. Results were confirmed through independent experiments with different parameter choices. B. Geometric Samplers

Fig. 6: Iterative refinement of the Taylor series bounds. At each step, the interval [L, H] tightens around the true value. The algorithm terminates when u falls outside the current bracket, returning True if u < L or False if u ≥ H.

the partial sums of the Taylor series expansion until a decision can be made. The first step samples a 128-bit uniformly random integer U ∈ {0, . . . , 2128 − 1} and interprets it as the rational u = U/2128 . In our implementation, this fraction is stored as a pair of arbitrary-precision integers (numerator and denominator) using Rust’s dashu::RBig type [7]. The algorithm exploits the Taylor series expansion: e−x =

∞ X (−x)k k=0

k!

=1−x+

x2 x3 − + ··· 2! 3!

Because this is an alternating series with decreasing terms, consecutive partial sums Sk−1 and Sk always bracket the true value e−x : one lies above and the other below. The algorithm initializes Sprev = 1 and Scurr = 1 − x, corresponding to the first two partial sums. At each iteration, it computes bounds L = min(Sprev , Scurr ) and H = max(Sprev , Scurr ), then compares the uniform sample u against these thresholds. If u < L, we know with certainty that u < e−x , so the algorithm returns True. If u ≥ H, we know u ≥ e−x , so it returns False. Otherwise, u falls within the undecided region, and the algorithm computes the next term until convergence. Figure 6 illustrates this process. At each step, the bounds converge toward e−x until u falls outside the shrinking interval, at which point a decision is made. VIII. VALIDATION R ESULTS After replacing the original bernoulli_exp1 function with our Taylor series implementation, we repeated the diagnostic pipeline to verify that the artifacts had been resolved. The results confirm that all samples now correctly fit their theoretical distributions. A. Bernoulli Samplers Figure 7 shows that both bernoulli_exp1 and bernoulli_exp now produce samples consistent with the theoretical model. From n = 106 samples, we observe an empirical success probability of p̂ = 0.6651 against a theoretical value of pth = 0.6643. The z-test yields z = 1.80 with p = 0.073, and the χ2 test gives χ2 = 3.22 with p = 0.073. While the p-values are somewhat low, the z-score sits within 2σ, indicating acceptable statistical fluctuation rather than systematic bias. This stands in stark contrast to the

The improvement propagates through the sampling hierarchy. Figure 8 shows the geometric samplers, both of which now pass goodness-of-fit tests. The slow geometric sampler continues to perform well (χ2 = 9.74 with 13 degrees of freedom, p = 0.715). The fast geometric sampler now adheres to the theoretical distribution (χ2 = 145.5 with 138 degrees of freedom, p = 0.315), with empirical moments closely matching theoretical values: µ̂ = 9.507 against µth = 9.508, 2 and σ̂ 2 = 99.94 against σth = 99.92. The periodic artifacts are no longer visible. C. Discrete Laplace Distribution Figure 8c shows the discrete Laplace distribution produced by the corrected sampler. The characteristic periodic artifacts (compare with Figure 1) have been completely eliminated. The distribution passes the χ2 goodness-of-fit test with χ2 = 275.4 on 276 degrees of freedom (p = 0.499), indicating excellent agreement with the theoretical model. The empirical variance 2 σ̂ 2 = 199.90 matches the theoretical value σth = 199.83, and the sample mean µ̂ = −0.017 is statistically indistinguishable from zero (z = −1.22, p = 0.224). These results confirm that the artifacts originated from the bernoulli_exp1 implementation and that our proposed fix seemingly resolves the issue. The corrected sampler produces discrete Laplace samples that are statistically indistinguishable from the theoretical distribution at the tested precision level. IX. D ISCUSSION Following disclosure with the OpenDP development team, the root cause of the artifacts was identified as an issue in the dashu library’s RBig / UBig operator [8]. Each coin flip requires evaluating Bernoulli(x/k), which internally calls sample_bernoulli_rational with x/k as a rational argument. This division is delegated to the faulty operator. The published release of dashu (v0.4, January 2024) contains a bug in this operator: given x = a/b and divisor k, the implementation computes the denominator as (b/g) · k instead of the correct b · (k/g), where g = gcd(a, k) [8]. The issue has been patched in the dashu master branch [9] but was never included in an official release. As a result, each call to sample_bernoulli_ratio nal within sample_bernoulli_exp1 receives a subtly incorrect probability. These errors depend on the value of k, producing a bias in the coin flip outcomes. When accumulated in the loop, they are then propagated through the geometric and discrete Laplace constructions, resulting in the artifacts. The OpenDP development team bypassed the call to the division by replacing it with a local implementation that computes b · (k/g) directly, preserving the full implementation of CKS20. The fix was integrated through pull request #2639, and pushed in the official release v0.14.2 [10].

(a) bernoulli_exp1: z = 1.80, p = 0.073.

(b) bernoulli_exp: z = 1.80, p = 0.073.

Fig. 7: Exponential Bernoulli samplers with the new implementation. Left: bernoulli_exp (p̂ = 0.6651, pth = 0.6643, χ2 = 3.22, p = 0.0726). Right: bernoulli_exp1 (p̂ = 0.6651, pth = 0.6643, χ2 = 3.22, p = 0.0726). Both now pass goodness-of-fit tests with χ2red ≈ 3 and p > 0.07, demonstrating correct sampling behavior. The insets on the plots show a zoomed region around the expected value, with error bars representing 1σ confidence intervals, showing how the observed values are within 3σ from the expected values.

(a) Geometric slow: χ2 = 9.74, p = 0.715.

(b) Geometric fast: χ2 = 145.5, p = 0.315.

(c) Discrete Laplace: z = −1.22, p = 0.22; χ2 = 2.8 × 102 , p = 0.5.

Fig. 8: Geometric and Discrete Laplace sampler distributions with the corrected implementation. (a) Slow implementation (mean z = −0.38, p = 0.71; χ2 = 9.7, p = 0.72). (b) Fast implementation (mean z = −0.18, p = 0.86; χ2 = 1.5 × 102 , p = 0.31). (c) Discrete Laplace distribution (mean z = −1.22, p = 0.22; χ2 = 2.8 × 102 , p = 0.5). The periodic artifacts visible in Figure 1 have been eliminated. All samplers pass goodness-of-fit tests.

We independently repeated the diagnostic pipeline of Section VI on the patched implementation, and all samplers in the chain pass the goodness-of-fit tests, with no visible artifacts in the final discrete Laplace distribution. We therefore continue this section with a brief performance benchmark of the three implementations and a remark on the practical privacy implications of the original faulty implementation. A. Performance Benchmarking While statistical correctness makes the OpenDP fix and the alternative Taylor-series sampler equivalent, their computational complexity differs. The original and the patched implementations follow the same CKS20 structure and differ only in where the rational division x/k is carried out (delegated to dashu in the former, performed locally in the latter) so we expect their computational cost to be essentially the same. The Taylor sampler, on the other hand, replaces the rejection scheme with a looping alternating-series construction. We measured the mean time per call to sample_discre te_laplace over a run of 106 samples each. The original sampler took 14.68±0.01µs, the Fixed sampler 15.08±0.01µs, and the Taylor sampler 17.45 ± 0.01µs.

The original and fixed costs sit within 2.7% of each other, indicating that the local implementation of the division does not change the complexity of the algorithm, as expected. The Taylor sampler is 15.7% slower than the fixed implementation, but remains informative: two structurally different fixes output the same empirical behavior, which suggests that the problem was correctly localized in sample_bernoulli_exp1.

B. Implications on the Privacy Guarantee A natural question is whether the artifacts described translate into a measurable degradation of the theoretical privacy guarantee. The privacy guarantees of CKS20 are established for noise drawn from the exact discrete Laplace distribution [4]. Once the underlying Bernoulli(exp(−γ)) primitive is biased, those bounds are no longer guaranteed to apply. To examine the consequences in practice, we employed the operational definition of Differential Privacy [1], and reconstructed the empirical Privacy Loss Random Variable (PLRV), a precise tool to measure the privacy guarantees of a differentially private mechanism [11]. For a randomized mechanism M and a pair of neighboring

datasets D1 , D2 , the privacy loss at output y is defined as   P (M(D1 ) = y) D1 ,D2 LM (y) = ln , P (M(D2 ) = y)

(a) Original

(b) Taylor series

1 ,D2 and a mechanism is ε-DP if |LD (y)| ≤ ε. The PLRV M thus provides a direct measurement of the privacy guarantee for each output: any bin whose value exceeds the ±ε band corresponds to an output for which the privacy guarantee is not held. We executed N = 106 counting queries on two neighboring datasets D1 and D2 of sizes 10,000 and 10,001 respectively. The counting query returns the number of entries in the database and has unit sensitivity (∆s = 1), since adding or removing a single entry changes the result by at most one. We set ε = 0.1, and estimated the PLRV as the bin-wise ratio of the resulting output histograms. We work with the ratio rather than its logarithm for computational simplicity, with the corresponding bound becoming e±ε instead of ±ε. Under correct ε-DP the empirical PLRV must lie, up to statistical fluctuations, within e±ε across all bins. Figure 9 compares the empirical PLRV obtained with the three implementations. For the original (faulty) sampler (Figure 9a), the fraction of bins whose PLRV deviates from the theoretical bound by more than 3σ is well above the ∼ 0.3% rate expected under Gaussian fluctuations, suggesting that the violations are systematic rather than random. Moreover, the bins that visibly deviate from the model can be traced back to the periodic artifacts identified in Figure 1: the same outputs that were over and under-sampled by the faulty Laplace distribution are those for which the ε-DP bound is violated. In other words, the artifacts visible in the noise distribution propagate to the PLRV and translate into visible violations of the ε-DP bound. Both the Taylor-series sampler (Figure 9b) and the official fix (Figure 9c) produce PLRV values that remain within the theoretical bounds at 3σ confidence across the full histogram, with exceedances compatible with the expected statistical fluctuations. We present this as a conjecture based on the empirical evidence above, not as a rigorous PLRV analysis. Turning the qualitative observation that the bound is violated into a quantitative statement about the effective privacy violation would require a dedicated privacy-loss accounting analysis, which we do not attempt here.

X. C ONCLUSION (c) OpenDP v0.14.2

Fig. 9: Empirical PLRV from the output distributions of two neighboring datasets (ε = 0.1, N = 106 , bin size = 1). Dashed orange lines highlight the theoretical bounds e±ε . Error bars are 1σ confidence intervals. The faulty implementation (a) shows points scattered beyond e±ε , while both corrections (b, c) remain within bounds at 3σ confidence level.

We have identified and resolved systematic artifacts in OpenDP’s discrete Laplace sampler through a systematic diagnostic methodology that isolated the root cause to numerical precision issues in the bernoulli_exp1 primitive. Our alternative implementation eliminates these artifacts at the tested precision level, in agreement with the official OpenDP patch released after our disclosure. Empirical analysis of the privacy loss suggests the artifacts also violate the e±ε bound under the faulty sampler, while both corrections restore it.

More broadly, this work demonstrates the importance of rigorous statistical validation for differential privacy implementations. Even mathematically correct algorithms can produce unexpected behaviors when translated to finite-precision systems, and subtle implementation details can compromise privacy guarantees or degrade utility. Our diagnostic approach provides a template for isolating similar issues in hierarchical sampling implementations and emphasizes the need for comprehensive testing beyond unit tests of individual components. This analysis was shared with the OpenDP development team, leading to the identification of the bug in the dashu rational arithmetic library [8], and for subsequent integration of extended statistical testing throughout the OpenDP library in response to this disclosure [12]. ACKNOWLEDGEMENTS We would like to thank Michael Shoemate2 of the OpenDP development team for his valuable feedback and contributions following this disclosure. This project was partially funded by the European Union’s “Next Generation EU” Program under the Piano Nazionale di Ripresa e Resilienza (PNRR) through NQSTI Spoke8. R EFERENCES [1] Cynthia Dwork and Aaron Roth. The algorithmic foundations of differential privacy. Found. Trends Theor. Comput. Sci., 9(3–4):211–407, August 2014. [2] Ilya Mironov. On significance of the least significant bits for differential privacy. In Proceedings of the 2012 ACM Conference on Computer and Communications Security, CCS ’12, pages 650–661, New York, NY, USA, 2012. Association for Computing Machinery. [3] OpenDP Project. OpenDP: An Open-Source Project for Differential Privacy. https://opendp.org/, 2025. [4] Clement Canonne, Gautam Kamath, and Thomas Steinke. Discrete Gaussian for differential privacy. Journal of Privacy and Confidentiality, 12(1), July 2022. [5] Michael Hay, Marco Gaboardi, and Salil Vadhan. A programming framework for OpenDP. In 6th Workshop on the Theory and Practice of Differential Privacy (TPDP 2020), 2020. https://projects.iq.harvard. edu/files/opendp/files/opendp programming framework 11may2020 1 01.pdf. [6] Peter Occil. Bernoulli factory algorithms. https://peteroupc.github.io/be rnoulli.html, 2025. Version dated 2025-11-06. [7] Jacob Zhong. dashu: A library set of arbitrary precision numbers. https: //lib.rs/crates/dashu, 2024. Rust crate. [8] Michael Shoemate. Incorrect rational division - Issue #57. https://gith ub.com/cmpute/dashu/issues/57, February 2026. GitHub issue, opened February 19, 2026. [9] Zhuo Zhang. Implement Div<RBig> for UBig and IBig. https: //github.com/cmpute/dashu/commit/62378fc84750bb6390ca7424ffca8 c46fc882566, February 2024. Commit to master branch, February 13, 2024. Fix not included in any official release. [10] OpenDP Project. Additive noise sampler bias, statistical testing (PR #2639). https://github.com/opendp/opendp/pull/2639, 2026. Merged in OpenDP v0.14.2, 2026-03-10. Accessed: 2026-04-28. [11] Damien Desfontaines. The privacy loss random variable. https://desf ontain.es/blog/privacy-loss-random-variable.html, March 2020. Ted is writing things (personal blog). [12] Michael Shoemate. Bias in additive noise samplers - Issue #2638. https: //github.com/opendp/opendp/issues/2638, February 2026. GitHub issue, opened February 25, 2026.

2 https://github.com/Shoeboxam

Related documents

• Detection and Resolution of Periodic Artifacts in OpenDP's Discrete Laplace Sampler#1157178
• Verification of OpenDP Transformations#347590
• 5 CFR Part 1630 — Privacy Act Regulations#504139
• TTPrint: Evidence-Grounded TTP Extraction via Diverge-then-Converge Verification#224399
• 36 CFR Part 1008 — Requests Under the Privacy Act#509072