fix(special_functions): gamma_func returns Gamma(|x|) for negative non-integer arguments — reflection formula never applied (core Rust)

[spm-lab-repository-audit-bot]

Summary

sparse_ir::special_functions::gamma_func (sparse-ir/src/special_functions.rs:33-121) returns Γ(|x|) for every negative non-integer argument instead of the correct value of the Gamma function at x. The reflection formula (Γ(x) = π / (sin(πx) Γ(1-x))) is computed but never applied, so both the sign and the magnitude are wrong over the whole negative domain.

Affected ownership layer: core Rust crate (sparse-ir). The function is public (pub fn gamma_func inside pub mod special_functions;, sparse-ir/src/lib.rs:39) and is not exposed through the C ABI, Python, or Fortran bindings (no spir_* export and no reference to it in sparse-ir-capi/, python/, or fortran/).

User impact

• Any downstream Rust caller of sparse_ir::special_functions::gamma_func(x) with x < 0, x ∉ ℤ, gets a silently wrong number. For large |x| the error is catastrophic: gamma_func(-12.5) returns Γ(12.5) ≈ 6.84e7 instead of Γ(-12.5) ≈ -1.84e-9.
• Secondary exposure: pub fn cyl_bessel_j(nu, x) (special_functions.rs:127) calls gamma_func(nu + 1.0) (:130). For negative non-integer order nu < -1.0 the argument to gamma_func is negative, so J_ν(x) also comes back with the wrong sign and magnitude. (Internal spherical_bessel_j callers always use non-negative orders and are unaffected.)
• Negative arguments are evidently an intended part of the domain: the code panics specifically for negative integers (s == 0.0, :39-41), which only makes sense if negative non-integers are meant to be handled by reflection.

Evidence

Static analysis (source is sparse-ir/src/special_functions.rs):

• :37-44 — for x < 0, the reflection prefactor is stored in s (s = sinpi(x); ... s *= x;) and then x = -x. From this point on x is always >= 0.
• :70 — the large-x branch guards the reflection with return if x < 0.0 { PI / (res * s) } else { res };. Because x was already negated at :42, x < 0.0 is always false, so the reflection PI / (res * s) is dead code and Γ(|x|) is returned.
• :96-120 — the small-x branch (x <= 11.5) applies no reflection at all; it returns z * p_val / q_val = Γ(|x|) directly and silently discards s. The while x < 0.0 loop at :102-105 is also dead code (it only ever runs for the already-negated x).

Runtime verification (standalone reproducer that copies the function body verbatim from :33-121, compared against Python math.gamma reference values):

gamma_func(-0.5)  = +1.7724538509055159   math.gamma(-0.5)  = -3.544907701811032     (Γ(0.5)  returned)
gamma_func(-1.5)  = +0.8862269254527579   math.gamma(-1.5)  = +2.3632718012073544
gamma_func(-2.5)  = +1.3293403881791368   math.gamma(-2.5)  = -0.9453087204829417    (Γ(2.5)  returned)
gamma_func(-3.5)  = +3.3233509704478426   math.gamma(-3.5)  = +0.27008820585226917   (Γ(3.5)  returned)
gamma_func(-0.1)  = +9.5135076986687306   math.gamma(-0.1)  = -10.686287021193193
gamma_func(-12.5) = +68421682.73278293    math.gamma(-12.5) = -1.8366064838592814e-09
gamma_func(-10.1) = +454760.7514415855    math.gamma(-10.1) = -2.213416583085618e-06
gamma_func( 0.5)  = +1.7724538509055159   math.gamma( 0.5)  = +1.7724538509055159     (sanity check, OK)

Every row shows gamma_func(x) == Γ(|x|) (sign and magnitude wrong). The existing test_gamma_function (special_functions.rs:279-288) only exercises positive inputs (1.0, 2.0, 3.0, 4.0, 0.5), so the defect is untested.

Note: the candidate bug report that led to this issue gave the expected values for -2.5 and -3.5 with the wrong sign (+0.945... and -0.675...); the correct reference values are -0.9453087204829417 and +0.27008820585226917 (math.gamma). This does not affect the conclusion, which is confirmed independently here.

Root cause

The reflection formula is intended (the comment at :42 and the s == 0.0 panic confirm it) but the branch that applies it was never wired up correctly for either code path:

1. Large-x branch: the guard x < 0.0 at :70 tests the already-negated x, so it is always false — the correct PI / (res * s) expression is unreachable.
2. Small-x branch: the stored s is never used; the branch always returns Γ(|x|).

special_functions.rs:32 documents this as "a direct port of the C++ gamma_func implementation". In the style of implementation this mirrors, the reflection decision is based on the original sign (or a saved flag/s != 0.0), not on the post-negation x. The port appears to have replaced that guard with x < 0.0 evaluated after x = -x, which can never be true. (The upstream C++ source was not available offline here, so this is inferred from the Rust code structure and behavior rather than a direct diff against libsparseir.)

Violated rules

Repository-local rules (more specific; cited in preference to shared rules):

• REPOSITORY_RULES.md:49-52 (Numerical Correctness And Precision): "Numerical behavior is part of the public contract. Changes to kernels ... require tests that check values, reconstruction, or residuals rather than only shapes or successful execution." The function's negative-domain behavior is silently wrong and untested.
• REPOSITORY_RULES.md:53-57: the same section requires covering "difficult regimes: small and large ... real and complex data" — the negative real domain is a representative regime not covered by test_gamma_function.
• REPOSITORY_RULES.md:270-283 (Tests And Documentation): core numerical behavior belongs in sparse-ir and a public behavior defect needs value-checking tests.

Shared rules:

• rules/rust/numerical.md:5-6 (Correctness): "Numerical algorithms need tests for representative values, edge cases, ... and error branches."
• rules/common/provenance.md:9-19: a close translation/port must record and preserve the upstream reference behavior; here the port appears to have mistranslated the reflection guard.

Suggested fix (for the eventual PR, not done here)

• Save the original sign before negating, and apply PI / (res * s) in the large-x branch when the input was negative.
• Apply the reflection to the small-x branch as well (same π/(sin(πx)·Γ(1-x)) identity), or restructure so a single post-processing step handles reflection for both branches.
• Extend test_gamma_function with negative half-integer values and a large negative value checked against analytic/math.gamma references.

Verification limitations

• Per the audit protocol, the full workspace test suite and full builds were not run. Verification is by (a) source inspection via CodeGraph and (b) a standalone rustc reproducer that copies the function body verbatim from special_functions.rs:33-121 and compares against math.gamma references. The internal callers (cyl_bessel_j at :130, spherical_bessel_j_small_args at :165) always pass non-negative arguments, so they were not exercised for negative gamma_func input.
• No comparison against the upstream C++ gamma_func (libsparseir) was possible offline; the mistranslation claim is inferred from the Rust code and the Cephes-style reflection idiom.

Proposed audit-rule revision (optional, not required)

• Consider adding a repository audit heuristic for "negate-then-guard-on-sign" patterns in special-function ports (e.g., .cursor/rules/): flag any if x < 0.0 (or x >= 0.0) branch whose condition refers to a variable that was sign-flipped earlier in the same function, and require that sign-sensitive special functions (gamma, reflection/parity branches) have value-checking tests for both signs of their input domain.