# \[ANN\] FastTanhSinhQuadrature.jl: High-Performance SIMD-Accelerated Numerical Integration

**URL:** <https://discourse.julialang.org/t/ann-fasttanhsinhquadrature-jl-high-performance-simd-accelerated-numerical-integration/135271>\
**Category:** Package Announcements\
**Created:** [January 26, 2026, 10:41am UTC](https://discourse.julialang.org/t/ann-fasttanhsinhquadrature-jl-high-performance-simd-accelerated-numerical-integration/135271 "2026-01-26T10:41:28Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![svretina](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/svretina/32/206157_2.png) [@svretina](https://discourse.julialang.org/u/svretina)\
**Post date:** [January 26, 2026, 10:41am UTC](https://discourse.julialang.org/t/ann-fasttanhsinhquadrature-jl-high-performance-simd-accelerated-numerical-integration/135271/1 "2026-01-26T10:41:28Z")

</div>

Hi everyone,

I’m excited to announce [**FastTanhSinhQuadrature.jl**](https://github.com/svretina/FastTanhSinhQuadrature.jl), a new package for high-precision numerical integration targeting performance-critical applications and high-precision requirements (Float64, Double64, BigFloat).

### Why another Tanh-Sinh package?

The main differentiator here is **speed and SIMD-friendliness**. Most implementations of Tanh-Sinh (Double Exponential) quadrature require checks for underflow or “breaking” the loop when weights become too small near the endpoints. These branches prevent the compiler from effectively using SIMD instructions.

This package follows the paper **Tanh-sinh quadrature for single and multiple  
integration using floating-point arithmetic** ([Meitner et al., 2020](https://arxiv.org/pdf/2007.15057.pdf)). By utilizing the specific window limits and spacing defined in the paper, we can pre-determine exactly where the weights would underflow. This allows us to run the integration loops without any conditional checks, making the code perfectly suited for `LoopVectorization.jl` and upcoming parallel implementations.

### Key Features

- **Arbitrary Precision Support** : Seamlessly works with `Float32`, `Float64`, `BigFloat`, and extended precision types like `Double64` (from `DoubleFloats.jl`).
- **High Performance** : Specialized `integrate1D_avx`, `integrate2D_avx`, `integrate3D_avx` routines utilize `LoopVectorization.jl` for maximum speed.
- **Multidimensional Support** : Built-in support for **1D** , **2D** , and **3D** integration domains.
- **Memory Efficiency** : Pre-compute quadrature nodes and weights once and reuse them for multiple integrations.
- **Singularity Handling** : Robust handling of functions with singularities at integration boundaries via `quad_split`.
- **Double Exponential Convergence** : Achieve machine precision with few points even for singular integrands.
- **Adaptive Integration** : Adaptive integration with `quad` and `quad_split` functions.
- **Optimal and maximal spacing** : Optimal spacing of nodes and weights for maximum accuracy.
- **Underflow/Overflow Handling** : Robust handling of underflow/overflow when generating nodes and weights.

### Convergence & Performance

The Tanh-Sinh method is famous for its “double exponential” convergence. As shown in the convergence plot, it reaches machine precision extremely rapidly even for challenging integrands:

 ![Convergence Plot](https://global.discourse-cdn.com/julialang/original/3X/c/d/cd6d730702cc08558d414d216d9aa3f4bfb61ff8.svg)

**Benchmarks (Intel Core Ultra 7 155U):**  
By removing the overhead of underflow checks and enabling SIMD, we see substantial improvements in execution time:

| Function | Domain | Points | TS (ns) | TS SIMD (ns) | GQ (ns) | Ratio (TS/GQ) | Ratio (TS SIMD/GQ) |
| --- | --- | --- | --- | --- | --- | --- | --- |
| exp(x) | [-1, 1] | 5 | 33.09 | 14.27 | 33.03 | 1.00 | 0.43 |
| exp(x) | [-1, 1] | 50 | 322.87 | 81.99 | 176.68 | 1.83 | 0.46 |
| exp(x) | [-1, 1] | 500 | 3374.75 | 782.08 | 1665.50 | 2.03 | 0.47 |
| sin(x)^2 | [-1, 1] | 5 | 40.39 | 23.38 | 32.63 | 1.24 | 0.72 |
| sin(x)^2 | [-1, 1] | 50 | 404.32 | 145.48 | 176.75 | 2.29 | 0.82 |
| sin(x)^2 | [-1, 1] | 500 | 4361.71 | 1409.10 | 1768.40 | 2.47 | 0.80 |
| 1/(1+25x^2) | [-1, 1] | 5 | 3.00 | 3.75 | 17.35 | 0.17 | 0.22 |
| 1/(1+25x^2) | [-1, 1] | 50 | 41.82 | 39.17 | 34.77 | 1.20 | 1.13 |
| 1/(1+25x^2) | [-1, 1] | 500 | 477.21 | 213.24 | 320.39 | 1.49 | 0.67 |
| sqrt(1-x^2) | [-1, 1] | 5 | 4.06 | 4.69 | 20.31 | 0.20 | 0.23 |
| sqrt(1-x^2) | [-1, 1] | 50 | 68.01 | 45.31 | 69.68 | 0.98 | 0.65 |
| sqrt(1-x^2) | [-1, 1] | 500 | 636.53 | 313.13 | 743.43 | 0.86 | 0.42 |
| x^2 | [-1, 1] | 5 | 2.01 | 3.04 | 19.86 | 0.10 | 0.15 |
| x^2 | [-1, 1] | 50 | 17.72 | 21.68 | 38.07 | 0.47 | 0.57 |
| x^2 | [-1, 1] | 500 | 213.31 | 54.97 | 273.19 | 0.78 | 0.20 |
| log(1-x) | [-1, 1] | 5 | 38.06 | 30.29 | 34.32 | 1.11 | 0.88 |
| log(1-x) | [-1, 1] | 50 | 371.88 | 187.29 | 214.70 | 1.73 | 0.87 |
| log(1-x) | [-1, 1] | 500 | 4117.38 | 1978.80 | 2117.80 | 1.94 | 0.93 |

### High-Level API: `quad`

The simplest way to integrate is using the `quad` function, which provides adaptive integration:

```julia
using FastTanhSinhQuadrature

# Integrate exp(x) from 0 to 1
val = quad(exp, 0.0, 1.0)
println(val) # ≈ e - 1 ≈ 1.7182818...

# Integrate over default domain [-1, 1]
val = quad(x -> 3x^2)
println(val) # ≈ 2.0

# Handle singularities with quad_split
f(x) = 1 / sqrt(abs(x)) # Singular at x=0
val = quad_split(f, 0.0, -1.0, 1.0) # Split at singularity
println(val) # ≈ 4.0

```

### Pre-computed Nodes

For repeated integrations, pre-compute nodes and weights once:

```julia
using FastTanhSinhQuadrature

# Generate nodes (x), weights (w), and step size (h)
x, w, h = tanhsinh(Float64, 80)

# Integrate multiple functions efficiently
f1(x) = sin(x)^2
f2(x) = cos(x)^2

res1 = integrate1D(f1, 0.0, π, x, w, h)
res2 = integrate1D(f2, 0.0, π, x, w, h)
println("Integrals: $res1, $res2") # Both ≈ π/2

```

### SIMD-Accelerated Integration

For `Float32`/`Float64`, use the `_avx` variants for maximum speed:

```julia
x, w, h = tanhsinh(Float64, 100)

# Standard integration
val1 = integrate1D(exp, x, w, h)

# SIMD-accelerated (2-3x faster)
val2 = integrate1D_avx(exp, x, w, h)

```

### High-Precision Integration

Switch to higher precision types like `BigFloat` or `Double64`:

```julia
using FastTanhSinhQuadrature, DoubleFloats

# Double64 precision (~32 decimal digits)
val = quad(exp, Double64(0), Double64(1); tol=1e-30)

# BigFloat precision (arbitrary)
setprecision(BigFloat, 256)
x, w, h = tanhsinh(BigFloat, 120)
val = integrate1D(exp, x, w, h)

```

### Multidimensional Integration (2D & 3D)

Use `StaticArrays` for defining integration bounds:

```julia
using FastTanhSinhQuadrature, StaticArrays

# 2D: Integrate f(x,y) = x*y over [-1,1] × [-1,1]
x, w, h = tanhsinh(Float64, 40)
low = SVector(-1.0, -1.0)
up = SVector(1.0, 1.0)
val = integrate2D((x, y) -> x * y, low, up, x, w, h)
println(val) # ≈ 0.0

# 3D: Integrate constant 1 over unit cube
val = quad((x, y, z) -> 1.0, [0.0, 0.0, 0.0], [1.0, 1.0, 1.0])
println(val) # ≈ 1.0

```

Feedback and PRs are welcome!

**Repo:** [GitHub - svretina/FastTanhSinhQuadrature.jl: Fast and high-precision numerical integration using Tanh-Sinh (Double Exponential) quadrature in Julia.](https://github.com/svretina/FastTanhSinhQuadrature.jl)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 26, 2026, 11:29am UTC](https://discourse.julialang.org/t/ann-fasttanhsinhquadrature-jl-high-performance-simd-accelerated-numerical-integration/135271/2 "2026-01-26T11:29:38Z")

</div>

Are you registering the package? I setup an Integrals.jl wrapper:

> <https://github.com/SciML/Integrals.jl/pull/310>
>
> \## Summary
> 
> \- Add support for \[FastTanhSinhQuadrature.jl\](https://github.com/svr…etina/FastTanhSinhQuadrature.jl) as a new integration backend via package extension
> \- FastTanhSinhQuadrature.jl provides high-performance SIMD-accelerated Tanh-Sinh (double exponential) quadrature
> \- Supports 1D, 2D, and 3D integration
> 
> Related announcement: https://discourse.julialang.org/t/ann-fasttanhsinhquadrature-jl-high-performance-simd-accelerated-numerical-integration/135271
> 
> \## Changes
> 
> \- Add \`FastTanhSinhQuadratureJL\` algorithm struct with \`tol\` and \`max\_levels\` parameters
> \- Create \`IntegralsFastTanhSinhQuadratureExt\` extension module in \`ext/\`
> \- Add to \`Project.toml\`: weakdeps, extensions, compat, extras, and test targets
> \- Add to interface tests with appropriate capability flags (\`nout=Inf\`, \`min\_dim=1\`, \`max\_dim=3\`, \`allows\_batch=false\`, \`allows\_iip=false\`)
> \- Export \`FastTanhSinhQuadratureJL\` from main module
> 
> \## Usage Example
> 
> \`\`\`julia
> using Integrals, FastTanhSinhQuadrature
> 
> \# 1D integration
> f(x, p) = cos(x)
> prob = IntegralProblem(f, (0.0, 1.0))
> sol = solve(prob, FastTanhSinhQuadratureJL())
> 
> \# 2D integration
> f2d(x, p) = x\[1\] \* x\[2\]
> prob2d = IntegralProblem(f2d, (\[0.0, 0.0\], \[1.0, 1.0\]))
> sol2d = solve(prob2d, FastTanhSinhQuadratureJL(tol=1e-10))
> \`\`\`
> 
> \## Note
> 
> FastTanhSinhQuadrature.jl is not yet registered in the General registry. Users will need to add it from GitHub:
> \`\`\`julia
> Pkg.add(url="https://github.com/svretina/FastTanhSinhQuadrature.jl")
> \`\`\`
> 
> \## Test plan
> 
> \- \[x\] Verified 1D integration works correctly
> \- \[x\] Verified 2D integration works correctly
> \- \[x\] Verified 3D integration works correctly
> \- \[x\] Verified \`reltol\` parameter is passed through correctly
> \- \[\] Interface tests pass (requires FastTanhSinhQuadrature.jl to be registered or manually added to CI)
> 
> 🤖 Generated with \[Claude Code\](https://claude.ai/code)

---

<div class="post-metadata">

**Author:** ![svretina](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/svretina/32/206157_2.png) [@svretina](https://discourse.julialang.org/u/svretina)\
**Post date:** [January 26, 2026, 11:43am UTC](https://discourse.julialang.org/t/ann-fasttanhsinhquadrature-jl-high-performance-simd-accelerated-numerical-integration/135271/3 "2026-01-26T11:43:30Z")

</div>

yes I made the registration to the Registry and it passes the automated checks.  
I was looking at the Integrals.jl to write a wrapper 🙂

---

<div class="post-metadata">

**Author:** ![svretina](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/svretina/32/206157_2.png) [@svretina](https://discourse.julialang.org/u/svretina)\
**Post date:** [January 26, 2026, 12:08pm UTC](https://discourse.julialang.org/t/ann-fasttanhsinhquadrature-jl-high-performance-simd-accelerated-numerical-integration/135271/4 "2026-01-26T12:08:42Z")

</div>

I could implement simd-friendly window limits for the rest of the Double Exponential rules, but the package name then would be missleading. There exists [DoubleExponential](https://github.com/machakann/DoubleExponentialFormulas.jl) already for 1D with checks in the loops breaking simd.  
I am open to suggestions 🙂