Stable exp_mod_normal

[avehtari]

This PR and PR text were assisted by Sol and Herbie

Summary

• Stabilize exp_mod_normal_lpdf, cdf, lcdf, and lccdf values and gradients in normal tails,
  cancellation regions, and the small-lambda limit.
• Apply equivalent CPU and OpenCL changes (OpenCL needing more lines), including corrected size checks, infinity handling, and
  CDF scale-gradient signs.
• Add independent fixed-reference regressions and CPU/OpenCL comparisons.
• Inlcudes helpers that can be used by other distribution functions

Motivation

The existing CDF family subtracts nearly equal normal-tail terms. This can produce substantially
wrong finite values and near-zero gradients at moderate normalized parameters. For example:

• lcdf: -0.62114 instead of -0.65315, with d_y ≈ -1.5e-14 instead of 0.76490.
• cdf: 0.53733 instead of 0.52041, with d_y ≈ -8e-15 instead of 0.39806.
• lccdf: -22.44429 instead of -22.00452, with d_y ≈ -9e-14 instead of -6.33910.

The LPDF separately evaluates log(erfc(r)) and exp(-r^2) / erfc(r). At r = 27.25, the current
implementation returns -Inf and an infinite gradient, although the finite values are approximately
-709.09642 and 37.56323.

The regression tests have specific failing examples.

Herbie analysis

Herbie was run on scalar value and analytic-gradient expressions copied from the current
implementation, using fixed seeds and targeted central, tail, and cancellation domains. Across three
seeds, it identified repeatable improvements of approximately:

• 8 bits for the LCDF value;
• 9–10 bits for the LCDF gradient;
• 7–8 bits for the LCCDF value;
• 11–12 bits for the LCCDF gradient.

Herbie consistently exposed the cancellation-prone normal terms and suggested erfc- and
FMA-based rewrites. Its output was used for triage only. The final implementation uses manually
derived log-space identities, asymptotic expansions, and independently evaluated regression
references.

Helper rationale

The repair repeatedly needs stable log(Phi(x)), the inverse Mills ratio phi(x) / Phi(x), and
positive-tail erfcx evaluation. These operations are implemented in an internal
normal_tail_utils.hpp helper. Exp-mod-normal-specific log differences, log sums, small-parameter
expansions, and normalized derivatives remain private to exp_mod_normal_utils.hpp.

The same normal-tail pair is required by skew_normal_lpdf and is directly applicable to lognormal
LCDF/LCCDF. Keeping it distribution-independent avoids duplicating sensitive tail approximations
without adding a public API or refactoring the existing normal CDF implementations.

Tests

Regression tests are included

Side Effects

Helper functions will make couple more PRs smaller

Release notes

• Improved stability of exp_mod_normal_lpdf, cdf, lcdf, and lccdf values and gradients

Checklist

• Copyright holder: Aki Vehtari

  The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
  - Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
  - Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

• the basic tests are passing

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • code passes the built in C++ standards checks (make cpplint)

• the code is written in idiomatic C++ and changes are documented in the doxygen

• the new changes are tested