🔢 mathematics

Performance Evaluation of an Adaptive Quadrature and a Double Exponential Formula Using Arbitrary-Precision Floating-Point Arithmetic

This paper evaluates the performance of Ninomiya's adaptive quadrature (AQE11D) and Takahasi-Moris' double exponential (DE) formula using arbitrary-precision arithmetic on 21 test problems, revealing that while AQE11D achieves high accuracy across all cases, the DE formula is significantly more efficient for endpoint singularities but fails on problems with interior singularities.

Original authors: Tomonori Kouya

Published 2026-08-14

📖 7 min read🧠 Deep dive

CC BY 4.0

Original authors: Tomonori Kouya

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

The Quest for Perfect Numbers

Imagine you are trying to measure the area of a strangely shaped pond. You could throw a net over it and count the fish, or you could try to calculate the exact shape using a ruler and a calculator. In the world of science and engineering, this is called numerical integration: the art of finding the exact area under a curve when you can't just use a simple formula. Usually, computers use a standard ruler called "double precision," which is incredibly accurate for most daily tasks, like building bridges or simulating weather. But sometimes, the pond has a bottom that drops off into an infinite abyss, or the water is so turbulent that the standard ruler breaks. These are "ill-conditioned" problems, often involving special functions or extreme precision needed for things like verifying mathematical proofs or calculating the behavior of subatomic particles.

To tackle these impossible shapes, scientists use arbitrary-precision arithmetic. Think of this not as a standard ruler, but as a magical, infinitely expandable tape measure. Instead of being stuck with a fixed number of markings, you can ask the computer to add as many decimal places as you need—50, 100, or even more—to get the answer just right. The paper we are exploring today uses a specific, super-powerful version of this tape measure (called MPFR) to test two different strategies for measuring these tricky areas. One strategy is like a smart, adaptive explorer who keeps zooming in on the tricky parts, while the other is a master of a special transformation that smooths out the bumps. The question is: when the math gets really, really hard, which explorer wins?

The Great Race: The Adaptive Climber vs. The Smooth Slider

In this study, the author, Tomonori Kouya, set up a digital race track with 21 different "Kahaner test problems." These are like 21 different obstacle courses, ranging from smooth, gentle hills to jagged cliffs with sharp spikes and even places where the ground suddenly disappears (singularities). The goal was to see which of two mathematical methods could measure the area of these shapes with extreme accuracy (specifically, getting the answer right to 50 or 100 decimal places) without taking too long or using up too much computing power.

The first contender is AQE11D, an "adaptive quadrature" method. Imagine this as a very careful hiker who carries a map that gets more detailed as they walk. When the hiker sees a flat, easy path, they take big steps. But when they hit a rocky, jagged section or a cliff edge, they stop, split the path in half, and take tiny, cautious steps to measure every nook and cranny. This method is incredibly thorough; it can handle almost any shape, including those with sudden jumps or sharp peaks inside the path.

The second contender is the Double Exponential (DE) formula. This is more like a master surfer who uses a special board to ride a wave. Instead of walking step-by-step, this method transforms the entire landscape. It takes a jagged cliff or a sharp peak and mathematically "smooths" it out into a gentle, rolling hill that decays incredibly fast. Once the shape is smoothed, the method can measure it with huge, sweeping strides, getting a result very quickly.

The Results: Thoroughness vs. Speed

When the race began, the results were a tale of two very different strengths.

The AQE11D hiker was the ultimate survivor. It managed to conquer all 21 obstacle courses, reaching the target accuracy even on the most difficult terrain. However, this thoroughness came at a steep price. When the hiker encountered a "strong endpoint singularity"—a cliff that drops off infinitely at the very start of the path, like the function $1/\sqrt{x}$ —they had to take an enormous number of tiny steps. To get the answer right to 100 decimal places on this specific cliff, the hiker had to take about 54 million steps (function evaluations). It was a marathon that took a long time, but they got there.

The DE surfer, on the other hand, was a speed demon on the smooth and cliff-edge courses. On that same infinite cliff ( $1/\sqrt{x}$ ), the surfer glided over it in just 743 steps, finishing the job in a fraction of the time. The DE formula was roughly 70,000 times faster than the hiker on this specific problem. It successfully conquered 18 out of the 21 courses, proving it is the champion for smooth paths and those with tricky edges.

However, the surfer had a weakness. When the obstacle course had a sudden jump or a sharp spike in the middle of the path (an interior discontinuity), or when the waves were too chaotic (strong oscillation), the surfer's board failed. The DE formula could not find a smooth way to ride these specific waves and simply gave up on three of the problems (numbers 2, 9, and 21). The adaptive hiker, with their ability to stop and zoom in, was the only one who could solve these middle-of-the-road puzzles.

The Parallel Power-Up

The study also tested what happens if you hire a whole team of hikers or surfers to work at the same time (parallel processing).

For the DE surfer, adding more team members worked beautifully. Because the surfer takes many large, uniform steps, the team could easily split the work. With 20 people, the surfer got about twice as fast, and on the hardest problems, the team could handle the massive number of steps efficiently.

For the AQE11D hiker, adding a team didn't help much. The hiker's strategy is to split the path in half, measure the new spots, and then decide what to do next. This decision-making process is a bottleneck. Even with 20 people, the hiker couldn't speed up much because they were stuck waiting for the leader to say, "Okay, now we split this specific tiny rock." In fact, on the hardest cliff problem, the whole team was held back by the single person stuck on that one spot, slowing down the entire group.

The Verdict

So, what is the takeaway for anyone trying to measure the unmeasurable?

If your problem involves smooth curves or tricky edges at the start or end, the Double Exponential formula is the clear winner. It is fast, efficient, and requires far fewer calculations to get a super-precise answer. It's the go-to tool for the majority of difficult integrals.

However, if your problem has a sudden jump or a sharp spike right in the middle of the road, the adaptive quadrature (AQE11D) is the only one that will get the job done, even if it takes a long time and millions of steps.

The paper concludes that while the adaptive method is a reliable "Swiss Army knife" that can solve any problem, the DE formula is a specialized "race car" that is vastly superior for the vast majority of cases, provided the track doesn't have a sudden pothole in the middle. The author suggests that future work might try to combine the best of both worlds or use powerful graphics cards (GPUs) to make the adaptive hiker run faster, but for now, knowing which tool to pick for the job is the most important lesson.

Technical Summary: Performance Evaluation of Adaptive Quadrature and Double Exponential Formula Using Arbitrary-Precision Arithmetic

Problem Statement
Numerical integration often requires precision beyond the standard IEEE 754 binary64 format, particularly for ill-conditioned problems, special function evaluation, and numerical verification. While the GNU MPFR library provides a de facto standard for correctly rounded arbitrary-precision arithmetic, there is a need to evaluate how different quadrature strategies perform when implemented with this library. Specifically, this study compares two distinct approaches: an adaptive quadrature method (AQE11D) and the Double Exponential (DE) formula, focusing on their ability to meet user-specified accuracy requirements (tolerances of $10^{-50}$ and $10^{-100}$ ) across a diverse set of test cases, including those with endpoint singularities, interior discontinuities, and strong oscillations.

Methodology
The author implemented both algorithms using the GNU MPFR library on a CPU with OpenMP parallelism.

AQE11D: This is an extension of Ninomiya's adaptive 9-point Newton–Cotes rule (AQNN9). It constructs an error estimation functional $\bar{e}$ $\overset{e}{ˉ}$ from 11 function values. To address the limitations of fixed-order rules, AQE11D employs a sequence of higher-order rules (13-, 19-, 27-, and 41-point) for unconverged subintervals. Crucially, for high-precision arithmetic, the algorithm was modified to:
1. Generalize the rounding-error noise floor from a double-precision constant to a quantity proportional to $2^{5-PB}$ (where $PB$ is the mantissa precision).
2. Reverse the order of anomaly detection tests to correctly distinguish between logarithmic singularities (arithmetic sequences) and discontinuities (constant sequences) at arbitrary precision.
3. Utilize semi-analytic processing for endpoint singularities (e.g., $1/\sqrt{x}$ ) to avoid evaluating the divergent endpoint directly.
Double Exponential (DE) Formula: This method applies a variable transformation $x(t) = c + r \tanh(\frac{\pi}{2} \sinh t)$ to transform the integrand into a function that decays double-exponentially. The integral is then approximated using the trapezoidal rule. To prevent cancellation errors or evaluation of divergent values near endpoints in high precision, the implementation computes the relative distance from the endpoint directly rather than the coordinate itself.
Experimental Setup: The algorithms were tested on Kahaner's 21 benchmark problems. The mantissa precision ($PB$) was set to 256 bits for $\epsilon_a = 10^{-50}$ and 512 bits for $\epsilon_a = 10^{-100}$ . Performance was measured by attained accuracy, number of function evaluations, serial execution time, and OpenMP parallel speedup (both intra-integral and task-level).

Key Results

Accuracy and Convergence:
- AQE11D: Successfully attained the target accuracy on all 21 test problems at both tolerance levels. However, for problems with strong endpoint singularities (specifically Problem 7, $1/\sqrt{x}$ ), the method required a massive number of function evaluations ( $\approx 5.4 \times 10^7$ ) at $10^{-100}$ tolerance. This is roughly $7 \times 10^4$ times more evaluations than the DE formula for the same problem.
- DE Formula: Converged on 18 of the 21 problems. It failed to converge on problems featuring interior discontinuities (Problem 2), strong oscillations (Problem 9), or sharp peaks (Problem 21). For the problems where it did converge, it required significantly fewer function evaluations and less time than AQE11D.
Execution Time:
- On problems with endpoint singularities (e.g., Problems 3, 7, 19), the DE formula was orders of magnitude faster. For Problem 7 at $10^{-50}$ , DE took 1.2 ms compared to 125.5 ms for AQE11D.
- On problems where DE failed to converge, AQE11D was the only viable option.
Parallel Performance:
- Intra-integral Parallelism: Parallelizing function evaluations within a single integral yielded little to no speedup for AQE11D due to the small number of points (six) evaluated per bisection and high synchronization overhead. Conversely, the DE formula showed moderate speedups (2–4x) on singular/oscillatory problems and up to 6x on unconverged problems where sample points reached $\approx 10^5$ .
- Task Parallelism: Parallelizing across the 21 problems showed AQE11D achieving a speedup of $\approx 5.0$ at $10^{-50}$ , dropping to 2.6 at $10^{-100}$ due to the serial bottleneck of Problem 7. The DE formula achieved a consistent speedup of $\approx 2.0$ across both precisions, as the unconverged problems could not be subdivided.

Significance and Claims
The paper claims to provide design guidelines for numerical integrators that must meet strict, user-specified accuracy requirements using arbitrary-precision arithmetic. The primary contribution is the empirical demonstration that:

AQE11D is a robust, universal solver capable of handling a wide variety of singularities and discontinuities, but it incurs a prohibitive computational cost for strong endpoint singularities at very high precision.
The DE Formula is highly efficient for smooth functions and endpoint singularities but lacks the adaptability to handle interior discontinuities or complex oscillatory behavior without modification.
Arbitrary-Precision Implementation: The study highlights specific algorithmic adjustments (such as noise floor generalization and test reordering) required to make adaptive quadrature stable at arbitrary precision, noting that standard double-precision heuristics do not directly translate.

The author concludes that while the DE formula is superior in speed and efficiency for its applicable domain, a hybrid or adaptive approach is necessary to cover the full spectrum of integration challenges. Future work is proposed to address load balancing, subinterval-level parallelization, and GPU acceleration using CUDA ports of MPFR/MPC.