ENH: Replace ITKNetlib/SLATEC with std and Eigen special functions

[hjmjohnson]

Replaces the vendored netlib/SLATEC special functions used by the
ITKStatistics chi-square and Student-t distributions with std::tgamma and
the Eigen-backed itk::Math::IncompleteGammaP / RegularizedIncompleteBeta,
and removes the now-unused ITKNetlib third-party module. Closes #6498.

Aligns with the VNL/netlib → Eigen numerics direction (#6403, #6230); follows
the deprecation/benchmark pattern of #6441 and #6454.

What changed

• Vendored Eigen's unsupported SpecialFunctions module (scalar
  numext::igamma, numext::betainc) and registered it in
  UpdateFromUpstream.sh.
• New itkSpecialFunctions.h (itk::Math::IncompleteGammaP,
  RegularizedIncompleteBeta).
• itkChiSquareDistribution/itkTDistribution: Γ → std::tgamma; the
  incomplete gamma/beta calls → the new Eigen-backed API.
• Removed ITKNetlib from ITKStatistics deps and deleted
  Modules/ThirdParty/Netlib (~3,800 lines of Fortran-era C).
• New itkSpecialFunctionsGTest (vs SciPy references); legacy
  itkTDistributionTest scavenged into itkTDistributionGTest.

Phase-5 benchmark: Eigen vs SLATEC (50-digit mpmath truth)

std is preferred where available (Γ → std::tgamma); Eigen backs the two
functions with no std equivalent. Both replacements are far more than accurate
enough for the probability distributions (the distribution GTests pass at
1e-14 CDF / 1e-12 inverse).

Function	Backend	max rel-err	mean rel-err	median ns	speedup
IncompleteGammaP	SLATEC	2.820e-14	3.171e-15	~445	—
IncompleteGammaP	Eigen	1.901e-14	1.832e-15	37.5	~12×
RegIncompleteBeta	SLATEC	7.463e-15	1.224e-15	~640	—
RegIncompleteBeta	Eigen	3.005e-14	2.946e-15	73.8	~8.7×

Eigen wins both axes for the incomplete gamma. For the incomplete beta it is
~8.7× faster with a last-digit accuracy trade (3.0e-14 vs 7.5e-15, both
tens-of-ULP) — accepted because #6498 retires the module entirely and the
difference is irrelevant at the distribution level. Apple M3 Max, clang 21,
-O3, single-thread, median of 25×20000.

[hjmjohnson]

Phase-5 proof: Eigen vs netlib/SLATEC special functions (issue #6498)

Verdict. IncompleteGammaP passes the gate on both axes — Eigen is more
accurate (worst-case 1.9e-14 vs 2.8e-14, mean 1.8e-15 vs 3.2e-15) and
~12× faster. RegularizedIncompleteBeta is ~8.7× faster but marginally
less accurate (worst-case 3.0e-14 vs 7.5e-15, mean 2.9e-15 vs
1.2e-15) — both at the tens-of-ULP level, far below any statistical
tolerance (the TDistribution GTest passes at 1e-14 CDF / 1e-12 inverse). Since
issue #6498 retires ITKNetlib entirely (these were private extern "C" engine
internals, not a kept public API), the last-digit beta difference does not
justify keeping an unmaintained Fortran module alive.

Two-axis comparison (lower is better on both)

Accuracy = max / mean relative error vs a 50-digit mpmath reference over the
sampled grid; time = median ns/call, single-thread.

Function	Backend	max rel-err	mean rel-err	Eigen≤SLATEC	median ns	speedup
IncompleteGammaP (P(a,x))	SLATEC	2.820e-14	3.171e-15	—	~445	—
IncompleteGammaP (P(a,x))	Eigen	1.901e-14 ✅	1.832e-15 ✅	36/48	37.5	~12× ✅
RegIncompleteBeta (I_x(a,b))	SLATEC	7.463e-15	1.224e-15	—	~640	—
RegIncompleteBeta (I_x(a,b))	Eigen	3.005e-14	2.946e-15	22/40	73.8	~8.7× ✅

Gamma grid: a ∈ {0.5,1,2,3.5,5,11,25,50} × x ∈ {0.1,0.5,0.9,1.0,1.5,3.0}·a
(48 points spanning the series/continued-fraction regimes). Beta grid:
(a,b) ∈ {(0.5,0.5),(2,3),(5,2),(11,0.5),(3,7),(25,25),(0.5,5),(8,4)} × x ∈
{0.1,0.3,0.5,0.7,0.9} (40 points across the symmetry crossover).

Environment

Apple M3 Max, macOS 26.5.1, Apple clang 21.0.0, -O3 -DNDEBUG,
single-thread (ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1), median of 25 reps ×
20000 inner iterations, load ≈ 5. SLATEC built from the pre-deletion sources
(HEAD~3) linked against libitkv3p_netlib; Eigen path = Eigen::numext::igamma
/ betainc from the vendored unsupported/Eigen/SpecialFunctions. Reference =
mpmath 1.4.1 at 50 decimal digits.

Reproducer

Documentation/.../phase5_bench/ (scratch): restore SLATEC .c from HEAD~3,
compile with -I.../v3p/netlib, link libitkv3p_netlib; driver scores both
backends against the mpmath grid. Run: ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1 ./bench.

[hjmjohnson]

@greptile review

[greptile-apps]

Greptile Summary

Replaces the vendored SLATEC/Netlib Fortran-era C special functions (dgami_, dbetai_, dgamma_) used by ITKStatistics with std::tgamma, std::lgamma, and a new Eigen-backed itkSpecialFunctions.h (itk::Math::IncompleteGammaP, RegularizedIncompleteBeta), then removes the now-empty ITKNetlib third-party module (~3,800 lines of code deleted).

• Correctness: All three mathematical mappings are exact equivalences — the chi-square CDF/PDF and the Student-t CDF/PDF produce the same mathematical expressions as before, with the gamma-ratio in TDistribution::PDF switched to exp(lgamma-lgamma) to avoid overflow at large degrees of freedom.
• Testing: Legacy itkTDistributionTest.cxx is replaced by itkTDistributionGTest.cxx covering the same CDF/InverseCDF reference values, the parameter-vector and separate-DOF overloads, wrong-parameter throws, and variance semantics; itkSpecialFunctionsGTest.cxx adds SciPy reference values, boundary identities, and symmetry checks.
• Upstream tracking: UpdateFromUpstream.sh registers the full unsupported/Eigen/SpecialFunctions tree so future Eigen bumps carry it automatically.

Confidence Score: 5/5

Safe to merge — all three mathematical translations are provably equivalent, the Netlib dependency is cleanly severed, and the new GTests verify both reference accuracy and API contracts.

The argument-order and regularisation of every replaced call site has been verified against the SLATEC signatures and Eigen's igamma/betainc conventions. The lgamma-difference in TDistribution::PDF is a strict improvement over the direct ratio at large dof. Tests pass at the tolerances stated in the benchmark table and the GTest assertions confirm expected boundary behaviour. No regressions were identified.

No files require special attention. The only minor open point is the static_assert in itkSpecialFunctions.h accepting long double, which Eigen does not specialise for — harmless for current callers but could produce confusing errors for future ones.

Important Files Changed

Filename	Overview
Modules/Core/Common/include/itkSpecialFunctions.h	New public header exposing IncompleteGammaP and RegularizedIncompleteBeta via Eigen numext; template wrappers with static_assert guard are correct
Modules/Numerics/Statistics/src/itkChiSquareDistribution.cxx	Replaces dgami_/dgamma_ SLATEC calls with std::tgamma and Math::IncompleteGammaP; argument mapping (dofon2, xon2) is mathematically equivalent
Modules/Numerics/Statistics/src/itkTDistribution.cxx	Replaces dbetai_/dgamma_ with Math::RegularizedIncompleteBeta and std::exp(lgamma-lgamma); argument order (pin,qin,bx) matches old SLATEC signature; lgamma-difference avoids overflow at large dof
Modules/Numerics/Statistics/test/itkTDistributionGTest.cxx	Full GTest replacement of itkTDistributionTest covering CDF/InverseCDF for dof=1 and dof=11, parameter-vector API, wrong-parameter throws, and variance semantics; tolerances match Eigen accuracy
Modules/Core/Common/test/itkSpecialFunctionsGTest.cxx	New GTest with SciPy reference values, boundary identities, beta symmetry, and float overload coverage for both new wrappers
Modules/Numerics/Statistics/itk-module.cmake	Removes ITKNetlib from ITKStatistics DEPENDS; ITKCommon (which now owns itkSpecialFunctions.h) is already a dependency
Modules/ThirdParty/Eigen3/UpdateFromUpstream.sh	Registers unsupported/Eigen/SpecialFunctions and its src/ tree in the upstream-sync script so future Eigen bumps carry the module

Flowchart

%%{init: {'theme': 'neutral'}}%%
flowchart TD
    subgraph ITKCommon ["ITKCommon (new)"]
        SF["itkSpecialFunctions.h\nIncompleteGammaP\nRegularizedIncompleteBeta"]
        E["Eigen::numext::igamma\nEigen::numext::betainc"]
        SF --> E
    end

    subgraph ITKStatistics
        CHI["itkChiSquareDistribution.cxx\nPDF: std::tgamma\nCDF: IncompleteGammaP"]
        TDIST["itkTDistribution.cxx\nPDF: exp(lgamma-lgamma)\nCDF: RegularizedIncompleteBeta"]
    end

    subgraph OLD ["Removed: ITKNetlib / SLATEC"]
        NET["dgamma_\ndgami_\ndbetai_"]
    end

    SF --> CHI
    SF --> TDIST
    OLD -. replaced .-> CHI
    OLD -. replaced .-> TDIST

    subgraph TESTS
        SFTEST["itkSpecialFunctionsGTest\n(SciPy reference, boundaries,\nbeta symmetry, float overload)"]
        TDTEST["itkTDistributionGTest\n(CDF/InvCDF dof=1,11,\nparam APIs, throws, variance)"]
    end

    SF --> SFTEST
    TDIST --> TDTEST
    CHI --> TDTEST

Loading

%%{init: {'theme': 'base', 'themeVariables': {"darkMode": true, "background": "#0d1117", "primaryColor": "#21262d", "primaryTextColor": "#e6edf3", "primaryBorderColor": "#8b949e", "lineColor": "#8b949e", "textColor": "#e6edf3", "edgeLabelBackground": "#161b22", "actorBkg": "#21262d", "actorBorder": "#8b949e", "actorTextColor": "#e6edf3", "actorLineColor": "#8b949e", "signalColor": "#8b949e", "signalTextColor": "#e6edf3", "noteBkgColor": "#373320", "noteBorderColor": "#d4a72c", "noteTextColor": "#f0e6c0", "labelBoxBkgColor": "#21262d", "labelBoxBorderColor": "#8b949e", "labelTextColor": "#e6edf3", "loopTextColor": "#e6edf3", "activationBkgColor": "#30363d", "activationBorderColor": "#8b949e"}}}%%
flowchart TD
    subgraph ITKCommon ["ITKCommon (new)"]
        SF["itkSpecialFunctions.h\nIncompleteGammaP\nRegularizedIncompleteBeta"]
        E["Eigen::numext::igamma\nEigen::numext::betainc"]
        SF --> E
    end

    subgraph ITKStatistics
        CHI["itkChiSquareDistribution.cxx\nPDF: std::tgamma\nCDF: IncompleteGammaP"]
        TDIST["itkTDistribution.cxx\nPDF: exp(lgamma-lgamma)\nCDF: RegularizedIncompleteBeta"]
    end

    subgraph OLD ["Removed: ITKNetlib / SLATEC"]
        NET["dgamma_\ndgami_\ndbetai_"]
    end

    SF --> CHI
    SF --> TDIST
    OLD -. replaced .-> CHI
    OLD -. replaced .-> TDIST

    subgraph TESTS
        SFTEST["itkSpecialFunctionsGTest\n(SciPy reference, boundaries,\nbeta symmetry, float overload)"]
        TDTEST["itkTDistributionGTest\n(CDF/InvCDF dof=1,11,\nparam APIs, throws, variance)"]
    end

    SF --> SFTEST
    TDIST --> TDTEST
    CHI --> TDTEST

Loading

_{Reviews (3): Last reviewed commit: "ENH: Add GoogleTests for itk::Math speci..." | Re-trigger Greptile}

[blowekamp]

Awesome! Great work!

Why are the stdlib gamma functions used directly but the Eigen gamma and beta functions are? Was there an issue with direct use of the Eigen functions?

[hjmjohnson]

Why are the stdlib gamma functions used directly but the Eigen gamma and beta functions are? Was there an issue with direct use of the Eigen functions?

@blowekamp There is no functional issue with the Eigen calls — Eigen::numext::igamma/betainc work fine called directly. The split is a deliberate API-boundary choice, not a workaround:

RULE: Use std:: whenever possible to minimize maintenance burden. (avoiding transitive Eigen exposure in the distribution TUs (cf. #6230). When Eigen simply aliases c++17 features (i.e. in some Eigen API predate pre c++11, or 14, or 17 was a local algorithm, now is simple std:: reference).

• std::tgamma is standard, stable, and self-documenting, so there's nothing to insulate — calling it directly is the simplest correct thing.

• igamma not in std::

• std::beta(x,y) (C++17) is the complete beta B(x,y) only; there's no incomplete/regularized form

• igamma/betainc live in Eigen's unsupported SpecialFunctions module, so I put them behind itk::Math::IncompleteGammaP / RegularizedIncompleteBeta to:
  a. Insulate the call sites from an unsupported, tersely-named third-party API — if Eigen changes or we later swap backends, only the wrapper moves.
  b. Name the regularized semantics explicitly (P(a,x) = γ(a,x)/Γ(a); regularized Iₓ(a,b)), which igamma/betainc don't convey at the call site.
  c. Confine the heavy unsupported/Eigen/SpecialFunctions include (Bessel + AVX/NEON/GPU packet headers) to one wrapper,
  d. Add a static_assert(std::is_floating_point_v<…>) so an integer-literal call gives a clear diagnostic instead of Eigen's deep-template THIS_TYPE_IS_NOT_SUPPORTED.

[hjmjohnson]

⛔ Blocked on #6514 (CI-ordering only — no code dependency).

The red Pixi-Cxx checks here are an unrelated IODCMTK test include-path failure introduced upstream by #6496 (itkDCMTKFileReader.h: fatal error: dcmtk/dcmdata/dcdict.h: No such file or directory), which currently breaks main itself. It is not from this PR — the Statistics/Common/Eigen build (incl. itkSpecialFunctionsGTest) compiles clean; only the DCMTK test driver, a module this PR doesn't touch, fails.

#6514 ("COMP,BUG: Fix IODCMTK test include paths", approved) fixes it. Once it lands I'll rebase this branch onto the fixed main so the builds are no longer confounded.

[hjmjohnson]

Rebased onto main (now includes #6514, which fixed the IODCMTK test include paths from #6496). The previously-red Pixi-Cxx builds were confounded by that unrelated DCMTK breakage; with #6514 in the base, CI should now reflect only this PR's changes. New tip: 8b27bee.

[hjmjohnson]

The mac failures are fixed in already merged PR.