[andrjohns]

Fixes #2881. Implements most of #3364: the normal, lognormal and exp_mod_normal lcdf/lccdf through the shared std_normal_lcdf implementation, and the OpenCL part.

Summary

This PR extracts the key value and gradient calculations from the std_normal_lcdf header into reusable functions that other normal-cdf-related functions can delegate to, removing a lot of duplicated/inconsistent implementations.

The value and gradient calculations are also updated to replace the taylor expansions with the Cody rational approximations (previously only used for the tail).

Two other bugs were found & fixed for var_value<Eigen::VectorXd> handling:

• as_array_or_scalar() now returns an array unchanged instead of re-wrapping it
• ref_type_if keeps an rvalue arena_matrix by value rather than referencing a dead temporary

The valid ranges of inputs (before over/underflow) across all functions are expanded:

function	develop value	branch value	develop gradient	branch gradient
normal_lcdf, std_normal_lcdf	[-5e30, 38.5]	[-1e154, 38.5]	(-1e308, 38.6]	(-1e308, 38.6]
normal_lccdf, std_normal_lccdf	[-38.5, 5e30]	[-38.5, 1e154]	[-38.6, 1e308)	[-38.6, 1e308)
normal_cdf, std_normal_cdf	[-37.5, 8.24]	[-38.5, 8.28]	[-37.5, 38.6]	[-38.5, 38.6]
Phi	[-37.5, 8.24]	[-38.5, 8.28]	[-38.6, 38.6]	[-38.6, 38.6]
lognormal_lcdf	[-38.5, 8.28]	[< -500, 38.5]	[-38.5, 37.6]	[< -500, 37.6]
lognormal_lccdf	[-8.28, 38.5]	[-38.5, > 500]	[-38.6, 38.5]	[-38.6, > 500]
lognormal_cdf	[-38.5, 8.28]	[-38.5, 8.28]	[-27.2, 37.6]	[-38.5, 37.6]
exp_mod_normal_lpdf	[-37.5, 1e308)	[-1e154, 1e308)	[-37.5, 1e308)	(-1e308, 1e308)
exp_mod_normal_lcdf	[-7.2, 37.9]	[-2e5, ~745]	[-8.2, ~745]	[-2e5, ~745]
exp_mod_normal_lccdf	[-7.2, ~745]	[-38.4, 1e308)	[-8.2, ~745]	[-38.6, 1e308)
exp_mod_normal_cdf	[-7.2, 37.9]	[-38.4, 37.9]	[-8.2, ~745]	[-38.4, ~745]
skew_normal_lpdf	[-38.5, 1e154]	[-1e154, 1e154]	[-38.5, 1e308)	[-5e307, 1e308)

While performance is either unchanged or significantly improved across inputs:

function	min	slowest case	max	fastest case
normal_lcdf	1.66x	double, N = 100000	3.56x	var_valueEigen::Matrix, N = 1000
normal_lccdf	1.75x	fvar, N = 100000	3.71x	var_value, N = 1000
normal_cdf	0.92x	fvar, N = 100000	2.13x	double, N = 1000
std_normal_lcdf	1.78x	fvar, N = 100000	5.08x	var_valueEigen::Matrix, N = 1000
std_normal_lccdf	2.41x	double, N = 100000	5.02x	var_valueEigen::Matrix, N = 1000
std_normal_cdf	1.35x	double, N = 100000	2.68x	var_valueEigen::Matrix, N = 1000
Phi	1.17x	fvar, N = 100000	1.65x	double, N = 1000
lognormal_lcdf	1.37x	var_value, N = 10000	2.60x	double, N = 1000
lognormal_lccdf	1.32x	fvar, N = 10000	2.43x	double, N = 1000
lognormal_cdf	1.29x	var_value, N = 10000	2.24x	double, N = 1000
exp_mod_normal_lpdf	1.78x	var_value, N = 10000	2.36x	double, N = 1000
exp_mod_normal_lcdf	1.04x	fvar, N = 100000	1.52x	double, N = 1000
exp_mod_normal_lccdf	1.03x	var_valueEigen::Matrix, N = 100000	1.51x	double, N = 1000
exp_mod_normal_cdf	1.00x	var_value, N = 100000	1.48x	double, N = 1000
skew_normal_lpdf	1.58x	var_value, N = 100000	2.49x	double, N = 1000

Full timing comparisons

double (values only)

function	N = 100	N = 1000	N = 10000	N = 100000
normal_lcdf	34.3 / 13.6 (2.52x)	34.5 / 13.4 (2.57x)	39.0 / 18.7 (2.08x)	40.0 / 24.2 (1.66x)
normal_lccdf	35.0 / 13.9 (2.52x)	34.8 / 13.6 (2.57x)	39.0 / 17.5 (2.24x)	40.2 / 19.4 (2.07x)
normal_cdf	28.3 / 13.7 (2.06x)	28.7 / 13.4 (2.13x)	28.4 / 16.9 (1.68x)	28.4 / 18.9 (1.50x)
std_normal_lcdf	36.3 / 14.3 (2.54x)	38.0 / 14.0 (2.71x)	39.8 / 17.8 (2.24x)	40.7 / 19.7 (2.06x)
std_normal_lccdf	36.8 / 13.4 (2.75x)	38.4 / 13.1 (2.93x)	39.9 / 15.4 (2.60x)	41.0 / 17.0 (2.41x)
std_normal_cdf	26.2 / 14.4 (1.82x)	26.8 / 14.1 (1.91x)	26.6 / 17.8 (1.50x)	26.6 / 19.7 (1.35x)
Phi	26.9 / 16.5 (1.64x)	27.2 / 16.5 (1.65x)	27.3 / 19.9 (1.37x)	27.2 / 20.8 (1.31x)
lognormal_lcdf	37.1 / 14.7 (2.53x)	36.8 / 14.1 (2.60x)	41.0 / 18.9 (2.17x)	42.8 / 20.8 (2.05x)
lognormal_lccdf	36.9 / 15.7 (2.36x)	36.7 / 15.1 (2.43x)	40.2 / 19.6 (2.05x)	42.3 / 21.6 (1.95x)
lognormal_cdf	32.0 / 15.3 (2.09x)	33.0 / 14.7 (2.24x)	36.1 / 19.2 (1.88x)	38.5 / 21.5 (1.79x)
exp_mod_normal_lpdf	29.3 / 13.4 (2.20x)	29.8 / 12.6 (2.36x)	31.1 / 14.6 (2.13x)	32.8 / 16.6 (1.97x)
exp_mod_normal_lcdf	57.5 / 38.5 (1.49x)	57.2 / 37.6 (1.52x)	57.7 / 49.9 (1.16x)	57.8 / 51.4 (1.12x)
exp_mod_normal_lccdf	57.6 / 38.8 (1.48x)	57.3 / 37.9 (1.51x)	57.6 / 50.6 (1.14x)	57.9 / 52.0 (1.11x)
exp_mod_normal_cdf	55.8 / 38.2 (1.46x)	55.6 / 37.4 (1.48x)	55.8 / 50.3 (1.11x)	56.0 / 51.8 (1.08x)
skew_normal_lpdf	30.9 / 12.9 (2.39x)	30.7 / 12.4 (2.49x)	33.2 / 17.0 (1.95x)	34.3 / 19.4 (1.77x)

var_value (Eigen::Matrix<var_value> y, var_value parameters, gradient)

function	N = 100	N = 1000	N = 10000	N = 100000
normal_lcdf	62.2 / 18.6 (3.34x)	62.4 / 17.8 (3.50x)	70.1 / 29.5 (2.38x)	70.6 / 34.0 (2.08x)
normal_lccdf	67.2 / 18.7 (3.59x)	67.1 / 18.1 (3.71x)	73.6 / 21.4 (3.44x)	74.5 / 34.4 (2.16x)
normal_cdf	35.9 / 18.7 (1.92x)	36.1 / 17.9 (2.01x)	36.0 / 21.4 (1.69x)	36.1 / 22.9 (1.57x)
std_normal_lcdf	59.9 / 17.3 (3.47x)	61.5 / 17.2 (3.59x)	67.5 / 20.2 (3.34x)	68.2 / 21.6 (3.15x)
std_normal_lccdf	64.7 / 16.6 (3.89x)	66.1 / 16.4 (4.03x)	71.5 / 18.7 (3.83x)	72.0 / 20.2 (3.57x)
std_normal_cdf	32.2 / 18.3 (1.76x)	32.7 / 17.1 (1.91x)	32.7 / 20.1 (1.63x)	32.6 / 21.6 (1.51x)
Phi	37.8 / 26.3 (1.43x)	38.3 / 26.6 (1.44x)	38.5 / 31.4 (1.23x)	38.7 / 32.2 (1.20x)
lognormal_lcdf	32.5 / 20.7 (1.57x)	31.8 / 20.1 (1.58x)	34.6 / 25.2 (1.37x)	36.5 / 26.5 (1.38x)
lognormal_lccdf	32.5 / 21.2 (1.53x)	31.6 / 20.5 (1.54x)	33.9 / 25.6 (1.32x)	35.8 / 26.9 (1.33x)
lognormal_cdf	30.6 / 20.8 (1.47x)	29.9 / 20.2 (1.48x)	32.7 / 25.4 (1.29x)	34.6 / 26.7 (1.29x)
exp_mod_normal_lpdf	32.6 / 18.1 (1.80x)	32.5 / 17.2 (1.89x)	33.8 / 19.0 (1.78x)	35.6 / 20.0 (1.78x)
exp_mod_normal_lcdf	61.9 / 47.9 (1.29x)	61.4 / 46.1 (1.33x)	61.2 / 58.0 (1.06x)	61.8 / 59.5 (1.04x)
exp_mod_normal_lccdf	62.0 / 48.2 (1.29x)	61.3 / 47.3 (1.30x)	61.0 / 58.7 (1.04x)	62.1 / 60.0 (1.04x)
exp_mod_normal_cdf	60.0 / 48.1 (1.25x)	59.7 / 45.9 (1.30x)	59.4 / 57.9 (1.03x)	59.6 / 59.6 (1.00x)
skew_normal_lpdf	34.7 / 18.1 (1.92x)	34.3 / 17.4 (1.97x)	36.9 / 22.8 (1.62x)	38.5 / 24.4 (1.58x)

var_valueEigen::Matrix (var_valueEigen::VectorXd y, var_value parameters, gradient)

function	N = 100	N = 1000	N = 10000	N = 100000
normal_lcdf	60.0 / 17.7 (3.39x)	60.1 / 16.9 (3.56x)	67.1 / 21.7 (3.10x)	67.9 / 33.0 (2.06x)
normal_lccdf	61.9 / 17.8 (3.48x)	61.3 / 17.1 (3.59x)	68.1 / 20.4 (3.35x)	68.9 / 33.3 (2.07x)
normal_cdf	35.0 / 17.8 (1.97x)	35.2 / 16.9 (2.08x)	35.2 / 20.3 (1.73x)	35.2 / 21.8 (1.61x)
std_normal_lcdf	59.0 / 12.1 (4.87x)	60.4 / 11.9 (5.08x)	65.9 / 15.3 (4.30x)	66.6 / 16.9 (3.95x)
std_normal_lccdf	60.2 / 12.5 (4.82x)	60.9 / 12.1 (5.02x)	66.1 / 15.1 (4.38x)	66.9 / 17.0 (3.94x)
std_normal_cdf	31.4 / 12.1 (2.59x)	31.8 / 11.9 (2.68x)	31.6 / 15.3 (2.06x)	31.8 / 16.9 (1.88x)
Phi	29.0 / 18.3 (1.59x)	28.9 / 18.1 (1.59x)	28.9 / 21.4 (1.35x)	29.1 / 22.5 (1.29x)
lognormal_lcdf	32.1 / 19.9 (1.61x)	31.1 / 19.1 (1.63x)	33.6 / 24.1 (1.39x)	35.4 / 25.5 (1.39x)
lognormal_lccdf	31.5 / 20.4 (1.54x)	30.6 / 19.3 (1.58x)	32.7 / 24.6 (1.33x)	34.5 / 25.9 (1.33x)
lognormal_cdf	29.8 / 19.9 (1.50x)	28.9 / 19.1 (1.51x)	31.5 / 24.3 (1.29x)	33.4 / 25.6 (1.31x)
exp_mod_normal_lpdf	31.8 / 17.3 (1.84x)	31.5 / 16.2 (1.95x)	32.7 / 18.0 (1.82x)	34.6 / 19.1 (1.81x)
exp_mod_normal_lcdf	60.8 / 46.9 (1.30x)	60.2 / 45.0 (1.34x)	60.1 / 56.9 (1.06x)	61.1 / 58.7 (1.04x)
exp_mod_normal_lccdf	60.8 / 47.1 (1.29x)	60.0 / 46.3 (1.29x)	60.0 / 57.7 (1.04x)	60.9 / 59.2 (1.03x)
exp_mod_normal_cdf	59.1 / 47.4 (1.25x)	58.4 / 45.1 (1.30x)	58.4 / 56.9 (1.03x)	59.2 / 58.5 (1.01x)
skew_normal_lpdf	34.1 / 17.2 (1.99x)	33.6 / 16.5 (2.03x)	36.0 / 21.8 (1.65x)	37.4 / 23.5 (1.59x)

fvar (tangent)

function	N = 100	N = 1000	N = 10000	N = 100000
normal_lcdf	60.4 / 18.1 (3.33x)	60.7 / 17.5 (3.47x)	75.4 / 36.9 (2.04x)	79.0 / 45.4 (1.74x)
normal_lccdf	62.4 / 18.9 (3.30x)	61.7 / 17.9 (3.44x)	68.4 / 21.4 (3.20x)	80.5 / 46.0 (1.75x)
normal_cdf	36.3 / 18.7 (1.94x)	35.8 / 17.7 (2.02x)	36.0 / 21.1 (1.70x)	36.5 / 39.6 (0.92x)
std_normal_lcdf	60.4 / 17.4 (3.47x)	61.4 / 17.1 (3.59x)	67.3 / 20.0 (3.36x)	68.1 / 38.3 (1.78x)
std_normal_lccdf	61.7 / 17.0 (3.64x)	61.7 / 16.3 (3.78x)	67.0 / 18.6 (3.60x)	68.1 / 20.2 (3.37x)
std_normal_cdf	32.5 / 17.6 (1.85x)	32.5 / 17.1 (1.90x)	32.6 / 20.0 (1.62x)	33.0 / 21.6 (1.53x)
Phi	32.6 / 24.1 (1.35x)	32.8 / 24.5 (1.34x)	33.0 / 27.6 (1.20x)	33.1 / 28.3 (1.17x)
lognormal_lcdf	32.8 / 20.8 (1.58x)	31.6 / 19.8 (1.59x)	34.5 / 25.0 (1.38x)	36.3 / 26.2 (1.39x)
lognormal_lccdf	32.4 / 21.3 (1.52x)	31.3 / 20.1 (1.56x)	33.7 / 25.5 (1.32x)	35.8 / 26.7 (1.34x)
lognormal_cdf	30.8 / 20.6 (1.49x)	29.6 / 19.9 (1.48x)	32.5 / 25.2 (1.29x)	34.5 / 26.5 (1.30x)
exp_mod_normal_lpdf	33.2 / 18.0 (1.85x)	32.7 / 16.8 (1.94x)	34.2 / 18.9 (1.80x)	35.9 / 19.9 (1.80x)
exp_mod_normal_lcdf	61.8 / 47.6 (1.30x)	60.4 / 45.8 (1.32x)	61.1 / 57.8 (1.06x)	62.2 / 60.0 (1.04x)
exp_mod_normal_lccdf	61.9 / 47.8 (1.29x)	60.3 / 46.9 (1.28x)	61.2 / 58.7 (1.04x)	61.9 / 59.9 (1.03x)
exp_mod_normal_cdf	60.1 / 47.8 (1.26x)	58.7 / 45.7 (1.29x)	59.4 / 57.8 (1.03x)	60.6 / 60.1 (1.01x)
skew_normal_lpdf	35.1 / 18.0 (1.95x)	34.1 / 17.2 (1.98x)	37.0 / 23.1 (1.60x)	38.6 / 24.0 (1.61x)

Tests

Additional tests exercising the increased numerical stability are added, as well new mix tests for the distributions

Side Effects

N/A

Release notes

Increased numerical stability and performance of the normal, std_normal, lognormal and exp_mod_normal (LC)CDF functions

Checklist

• Copyright holder: (fill in copyright holder information)

  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)
  • docs build, (make doxygen)
  • 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

[andrjohns]

This duplicates a bunch of the implementation over in #3405, so will make some suggestions/comments on that one and then rebase this once that PR lands

[avehtari]

I guess we should coordinate more as I also was preparing PR for std_normal and normal

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.04	0.04	0.98	-2.01% slower
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.16	0.17	1.0	-0.23% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.86	0.88	0.98	-2.07% slower
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	1.01	1.25% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.1	1.05	4.43% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	8.32	8.31	1.0	0.12% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.69	3.71	1.0	-0.34% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	33.64	33.02	1.02	1.83% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.38	0.99	-0.85% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.86	6.86	1.0	0.12% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	133.47	135.53	0.98	-1.54% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.21	3.2	1.0	0.45% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.6	0.59	1.01	1.07% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.27	25.0	1.01	1.07% faster
performance.compilation	393.17	392.02	1.0	0.29% faster
Mean result: 1.0026527691698786

Jenkins Console Log
Jenkins Build Stages
Commit hash: 98c64f32efad1c4821e30c9f76801a0d0c89d086

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           52 bits physical, 57 bits virtual
CPU(s):                                  192
On-line CPU(s) list:                     0-191
Thread(s) per core:                      2
Core(s) per socket:                      48
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              25
Model:                                   17
Model name:                              AMD EPYC 9474F 48-Core Processor
Stepping:                                1
Frequency boost:                         enabled
CPU MHz:                                 1497.515
CPU max MHz:                             4114.4229
CPU min MHz:                             1500.0000
BogoMIPS:                                7189.39
Virtualization:                          AMD-V
L1d cache:                               3 MiB
L1i cache:                               3 MiB
L2 cache:                                96 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-47,96-143
NUMA node1 CPU(s):                       48-95,144-191
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Not affected
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Enhanced / Automatic IBRS; IBPB conditional; STIBP always-on; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Mitigation; Clear CPU buffers
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good amd_lbr_v2 nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 pcid sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba perfmon_v2 ibrs ibpb stibp ibrs_enhanced vmmcall fsgsbase bmi1 avx2 smep bmi2 erms invpcid cqm rdt_a avx512f avx512dq rdseed adx smap avx512ifma clflushopt clwb avx512cd sha_ni avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local user_shstk avx512_bf16 clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin cppc arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif x2avic v_spec_ctrl vnmi avx512vbmi umip pku ospke avx512_vbmi2 gfni vaes vpclmulqdq avx512_vnni avx512_bitalg avx512_vpopcntdq la57 rdpid overflow_recov succor smca fsrm flush_l1d debug_swap

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

[avehtari]

I reviewed this PR (with help from Claude). Looks great

Accuracy against mpmath

Worst ulp over 601 points per range:

range	value	gradient
deep tail [-300, -40]	1.69	0.85
[-40, -10]	1.51	0.87
[-10, -5.66]	1.60	0.85
seam [-5.8, -5.5]	1.54	3.12
[-5.66, -1]	1.71	5.92
seam [-0.7, -0.6]	2.37	2.75
near zero [-0.6, 0.6]	1.86	2.36
[0.6, 4]	7.25	4.81
[4, 10]	30.37	31.53

For comparison, develop's worst in-range relative error is 7.6e-06 to 5.6e-05, which is 3.4e+10 to 2.5e+11 ulp, across every branch except the far Mills tail. So this PR improves the gradient by about ten orders of magnitude.

erfcx looks slightly faster on the lower tail

Beyond removing the duplicated Cody coefficients, #3405 may improve speed. Comparing your internal std_normal_erfcx path against calling erfcx directly, value and gradient computed together:

range (scaled)	CPU	GPU
far tail [-40, -10]	1.16	1.35
Cody [-10, -4]	1.10	2.21
Cody edge [-4, -2.5]	1.06	1.03
interior [-2.5, -0.1]	1.11	1.93

Ratios above 1 mean the direct erfcx call is faster.

[spinkney]

awesome! thanks @andrjohns for putting this together, it's much needed

[avehtari]

Two more functions get the benefit with no extra work

Two functions that benefit from this std_normal_lcdf PR that are missing from the PR text are

function	calls	OpenCL version
ordered_probit_lpmf	std_normal_lcdf x4	none
von_mises_cdf	normal_cdf	none

One function benefits from erfcx directly and is not covered here

wiener4_lcdf_unnorm line 128 is:

return std_normal_lcdf(-x) - std_normal_lpdf(x);

It calls std_normal_lcdf, so it looks like it inherits the repair., but doesn't. That expression is the log Mills ratio, and both terms approach -x^2/2 while their difference stays O(log x), so the subtraction cancels no matter how accurate the inputs are. Assuming both inputs are already correctly rounded, which is the best case after this PR:

x	log Phi(-x)	log phi(x)	result	as a difference	via erfcx
10	-53.2	-50.9	-2.3123	4.6 ulp	0.4 ulp
30	-454.3	-450.9	-3.4023	74.0 ulp	0.0 ulp
100	-5005.5	-5000.9	-4.6053	195.7 ulp	0.3 ulp
300	-45006.6	-45000.9	-5.7038	6132.2 ulp	0.2 ulp

The loss grows without bound. log Phi(-x) - log phi(x) equals log(sqrt(pi/2) * erfcx(x * INV_SQRT_TWO)) exactly, so one erfcx call removes the subtraction. Out of scope for this PR, but worth recording as the one which calls std_normal_lcdf but not benefiting from this PR.

[SteveBronder]

@andrjohns this is very cool! Because of how tied in this and #3405 are do you want to handle the #3405 review and then I'll review this once that is in?

[andrjohns]

@andrjohns this is very cool! Because of how tied in this and #3405 are do you want to handle the #3405 review and then I'll review this once that is in?

Yep sounds like a plan to me!

[SteveBronder]

@andrjohns ping me when this is ready for review

[SteveBronder]

@andrjohns is this ready for me to look at?

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.23	0.23	1.01	1.31% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	0.96	-4.44% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.32	6.31	1.0	0.13% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.98	-2.0% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.45	8.43	1.0	0.23% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.49	4.46	1.01	0.69% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	165.47	165.13	1.0	0.2% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.89	0.89	1.0	-0.08% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.7	0.71	0.99	-0.92% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	42.53	42.72	1.0	-0.45% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.6	0.6	0.99	-0.74% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	21.08	20.78	1.01	1.43% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	-0.4% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.19	3.19	1.0	0.22% faster
performance.compilation	386.86	387.1	1.0	-0.06% slower
Mean result: 0.9969507030598858

Jenkins Console Log
Jenkins Build Stages
Commit hash: 1a419c76226745cd916ab56e54b60394f2ebbdff

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           43 bits physical, 48 bits virtual
CPU(s):                                  256
On-line CPU(s) list:                     0-255
Thread(s) per core:                      2
Core(s) per socket:                      64
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              23
Model:                                   49
Model name:                              AMD EPYC 7742 64-Core Processor
Stepping:                                0
Frequency boost:                         enabled
CPU MHz:                                 1496.992
CPU max MHz:                             3416.0681
CPU min MHz:                             1500.0000
BogoMIPS:                                4491.85
Virtualization:                          AMD-V
L1d cache:                               4 MiB
L1i cache:                               4 MiB
L2 cache:                                64 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-63,128-191
NUMA node1 CPU(s):                       64-127,192-255
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Mitigation; untrained return thunk; SMT enabled with STIBP protection
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Retpolines; IBPB conditional; STIBP always-on; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Not affected
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba ibrs ibpb stibp vmmcall fsgsbase bmi1 avx2 smep bmi2 cqm rdt_a rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif v_spec_ctrl umip rdpid overflow_recov succor smca sev sev_es

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

[avehtari]

This PR fixes #2881

[avehtari]

This PR also implements most of #3364: the normal, lognormal and exp_mod_normal lcdf/lccdf through the shared std_normal_lcdf implementation, and the OpenCL part.

[avehtari]

std_normal_cdf, normal_cdf and lognormal_cdf from erfc

The original PR also computes these three as exp(lcdf), and with double arguments they take 1.4 to 1.9 times as long as on develop.

I pushed commit 0d4ac36, which computes them from erfc in the same way as Phi in my previous comment. Each factor is 0.5 * erfc(-z / sqrt(2)), the cdf is their product, and each partial is the cdf times phi(z) / Phi(z), divided by sigma (and by y for lognormal_cdf).

A new helper, internal::std_normal_cdf_value_ratio in std_normal_lcdf_impl.hpp, returns Phi(z) and phi(z) / Phi(z) for a scalar or an Eigen array, generic in the scalar type so that nested autodiff works. Where Phi underflows it returns a ratio of 0: the product is 0 there anyway, and exp(-z^2/2) / Phi would be 0 / 0. The three functions keep the argument checks and the broadcasting of their lcdf implementations. I did not change the OpenCL versions, which still use exp(lcdf).

Accuracy

Scalar calls with mu = 0.3, sigma = 1.7, all arguments var, z = (y - mu) / sigma from -38 to 8 in steps of 0.25. Worst relative error in ulp against mpmath (100 digits; partials by mpmath.diff), develop / original PR / after 0d4ac36:

function	z	value	d/dy
std_normal_cdf	[-38, -30)	1500 / 950 / 1500	1000 / 950 / 1.7
std_normal_cdf	[-10, -3)	3.5e5 / 51 / 89	65 / 51 / 1.5
std_normal_cdf	[-3, 3)	60 / 6.3 / 9.8	6.3 / 3.3 / 3.2
normal_cdf	[-30, -10)	1700 / 1200 / 1100	1700 / 1200 / 990
normal_cdf	[-10, -3)	3.5e5 / 120 / 160	63 / 120 / 87
normal_cdf	[-3, 3)	59 / 5 / 11	10 / 4.6 / 5.9
lognormal_cdf	[-30, -10)	1700 / 1200 / 1100	returns 0 / 1200 / 990
lognormal_cdf	[-10, -3)	62 / 120 / 150	64 / 110 / 85
lognormal_cdf	[-3, 3)	12 / 4.7 / 10	9.2 / 3.2 / 4.9

The d/dmu and d/dsigma errors have the same size as d/dy, except near z = 0, where the lognormal_cdf d/dsigma crosses zero; its absolute error there is 2e-18 in all three versions. The commit keeps the fixes of the original PR: develop's 3.5e5 ulp values in std_normal_cdf and normal_cdf, and the zero lognormal_cdf gradients below z = -10. No result is NaN or infinite. Compared with exp(lcdf), the values are up to 2 times less accurate, from the rounding of the erfc argument that Phi has too, and the gradients are equal or better.

Speed

Same machine and flags as for Phi, three trees alternated, 3 runs × 3 rounds. K is the number of elements per call; with var, all arguments are var and one reverse pass covers 4096 elements. Time relative to develop, original PR / after 0d4ac36:

function	range	K	double	var
std_normal_cdf	[-3, 3]	1	1.92 / 0.98	1.11 / 1.02
std_normal_cdf	[-3, 3]	64	1.27 / 0.97	0.89 / 0.95
std_normal_cdf	[-8, -3]	1	1.48 / 0.89	1.02 / 1.02
normal_cdf	[-3, 3]	1	1.72 / 1.02	1.06 / 1.01
normal_cdf	[-3, 3]	64	1.10 / 0.91	1.19 / 1.10
normal_cdf	[-8, -3]	1	1.45 / 0.90	0.98 / 0.96
lognormal_cdf	[-3, 3]	1	1.50 / 1.02	1.05 / 1.06
lognormal_cdf	[-3, 3]	64	0.85 / 0.70	1.50 / 1.44
lognormal_cdf	[-8, -3]	1	1.45 / 1.01	1.02 / 1.03

The round-to-round spread reached 14 % in some rows, so differences below about 0.1 are not significant. With double arguments all three are back at develop speed. One slowdown is left: lognormal_cdf of a vector with var arguments takes about 1.5 times as long as on develop, with either form, so the extra time is not in the erfc or exp(lcdf) code. I have not found where it is.

Tests

On this PR's head bf5861a plus 0d4ac36:

• CPU, 19 unit suites, all pass: prim and mix std_normal_cdf_test and normal_cdf_test, mix lognormal_cdf_test (with this PR's tail and y = 0 checks), and the std_normal_lcdf_impl, Phi, inv_Phi, lcdf, lccdf, exp_mod_normal_cdf and skew_normal suites from before.
• The generated test/prob tests of the three cdfs (runTests.py test/prob/normal test/prob/lognormal test/prob/std_normal -f _cdf_0): 17 test programs, 12000 cases, all pass.
• OpenCL on a Tesla V100: opencl/rev std_normal_cdf_test, normal_cdf_test, lognormal_cdf_test and Phi_test pass. They compare the OpenCL results with the changed CPU code.
• A compile check instantiates the three functions with fvar<fvar<var>>, std::vector<var>, Eigen vectors of fvar<double>, var_value column and row vectors, and mixed scalar and vector arguments.
• clang-format 10 and cpplint report nothing on the changed files.

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.22	0.23	0.98	-2.14% slower
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	1.0	-0.29% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.3	6.25	1.01	0.74% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.96	-3.85% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.23	7.04	1.17	14.39% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.45	4.47	0.99	-0.5% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	217.21	218.34	0.99	-0.52% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.89	0.89	1.0	0.16% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.61	0.62	1.0	-0.2% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	42.45	42.64	1.0	-0.43% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.59	0.59	1.01	0.98% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	20.93	20.82	1.01	0.54% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.01	0.65% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.17	3.17	1.0	-0.02% slower
performance.compilation	365.45	365.67	1.0	-0.06% slower
Mean result: 1.0080687934934953

Jenkins Console Log
Jenkins Build Stages
Commit hash: a2236819d562837e97a9bfe97dc9ebd6d8613729

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           43 bits physical, 48 bits virtual
CPU(s):                                  256
On-line CPU(s) list:                     0-255
Thread(s) per core:                      2
Core(s) per socket:                      64
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              23
Model:                                   49
Model name:                              AMD EPYC 7742 64-Core Processor
Stepping:                                0
Frequency boost:                         enabled
CPU MHz:                                 1496.647
CPU max MHz:                             3416.0681
CPU min MHz:                             1500.0000
BogoMIPS:                                4491.85
Virtualization:                          AMD-V
L1d cache:                               4 MiB
L1i cache:                               4 MiB
L2 cache:                                64 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-63,128-191
NUMA node1 CPU(s):                       64-127,192-255
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Mitigation; untrained return thunk; SMT enabled with STIBP protection
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Retpolines; IBPB conditional; STIBP always-on; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Not affected
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba ibrs ibpb stibp vmmcall fsgsbase bmi1 avx2 smep bmi2 cqm rdt_a rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif v_spec_ctrl umip rdpid overflow_recov succor smca sev sev_es

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

[avehtari]

Temporaries that are read more than once are now evaluated once

The vector lognormal_cdf slowdown from my previous comment has a simple cause, and the same cause made three lcdfs slower after this PR was opened. Commits f4cd10e and c30ee6f fix it. They change only when Eigen evaluates an expression, so all results are bitwise identical.

In lognormal_cdf, z = (log(y) - mu) / sigma was a lazy Eigen expression. Eigen evaluated it when the kernel copied it into an array, and again for the sigma partial, so each element paid for log twice. lognormal_lcdf_impl, which lognormal_lccdf also uses, had the same z. f4cd10e wraps both in eval(...).

To find the lcdf slowdown, I timed every commit of this PR with vectors of 1000 elements and var arguments. One commit is responsible: 9a8deff ("Refactor normal log-CDF dispatch and density expressions") replaced evaluated temporaries with lazy expressions. Among them are z and scaled_slope in normal_lcdf_impl and lognormal_lcdf_impl, and log_b, lp, s_b and q in exp_mod_normal_lcdf_impl. Each of these is read two to four times, and Eigen repeats the work at every read: in exp_mod_normal_lcdf, lp = log_diff_exp(log_a, log_b) ran three times per element. c30ee6f evaluates these temporaries once, as the code did before 9a8deff. With double arguments each temporary is read once, so those paths did not change.

ns per element, Eigen::Matrix<var> of 1000 elements with var parameters, value and gradient, z uniform on [-3, 3] (Xeon E5-2680 v3, g++ 14.3.0, Stan default flags, 2 runs × 3 rounds, spread at most 3.7 %):

function	develop	this PR at a28c1c1	before f4cd10e	after c30ee6f
std_normal_lcdf	97.5	42.7	37.5	37.1
normal_lcdf	104.8	48.9	55.9	43.0
lognormal_lcdf	59.3	59.8	65.8	54.7
exp_mod_normal_lcdf	91.2	129.0	218.5	115.3

std_normal_lcdf did not have the defect; its row shows the effect of the shared kernel alone. For lognormal_cdf of 64-element vectors with var arguments, the time relative to develop went from 1.29 to 1.06 with f4cd10e; when only y is var, so that no sigma partial is computed, it was already 1.02.

exp_mod_normal_lcdf is still slower than develop (0.79 of develop's speed with var, 0.72 with double). It evaluates two log cdfs and a log_diff_exp per element. Develop's form is faster but returns NaN or infinite values in large parts of the parameter space.

9a8deff also made exp_mod_normal_lpdf and skew_normal_lpdf 13 to 17 % slower with var, but there the faster kernel from 5f33a10 already makes up the difference, so I did not change them. neg_binomial_2_lpmf has the same pattern with log(mu + phi). Develop has it too, and #3434 already changes the phi partial of that function, so I will fix it there.

Speed of this PR compared with develop and with the table in the description

In short, the PR as opened made the lcdfs faster than develop but Phi and the cdfs slower. Commit 9a8deff then slowed down the var paths. With the commits from my comments, the PR is now faster than develop for every function in the table except exp_mod_normal_lcdf, apart from exp_mod_normal_lpdf and skew_normal_lpdf with var, which are at develop's time. The larger gains in the description most likely come from a slow erf/erfc in the C library of the machine used for it. Commit 9b38ecb removes the dependence of Phi and the cdfs on those two functions.

All tables give the time relative to develop: develop is 1.00, and a value below 1.00 means faster. The setup follows the description: vectors of N = 1000, Eigen::VectorXd with double parameters, and Eigen::Matrix<var> with var parameters and the gradient. The description does not give its inputs; I used z uniform on [-3, 3]. I ran all four versions in one job on an AMD EPYC 9965 (Zen 5, one whole node), with g++ 11.5, glibc 2.34 and the Stan default flags, 2 runs × 2 rounds, alternated; the spread between runs was at most 5.5 %. The columns are:

• description: 1 / the speedup in the description (its N = 1000 column);
• as opened: a28c1c1;
• before my commits: 686fc67, which includes 9a8deff;
• now: 9b38ecb, which includes all my commits (5f33a10 to 9b38ecb).

double:

function	description	as opened	before my commits	now
Phi	0.61	2.27	2.18	0.73
std_normal_cdf	0.53	1.44	1.43	0.82
normal_cdf	0.47	1.35	1.33	0.77
lognormal_cdf	0.45	1.03	1.07	0.77
std_normal_lcdf	0.37	0.93	0.91	0.80
normal_lcdf	0.39	0.89	0.88	0.78
lognormal_lcdf	0.38	0.82	0.85	0.74
exp_mod_normal_lcdf	0.66	1.62	1.67	1.62
exp_mod_normal_lpdf	0.42	0.98	0.93	0.83
skew_normal_lpdf	0.40	1.04	1.00	0.88

var (value and gradient):

function	description	as opened	before my commits	now
Phi	0.69	1.57	1.56	0.80
std_normal_cdf	0.52	1.00	1.02	0.63
normal_cdf	0.50	0.95	1.19	0.73
lognormal_cdf	0.68	1.26	1.72	0.93
std_normal_lcdf	0.28	0.43	0.44	0.39
normal_lcdf	0.29	0.44	0.56	0.40
lognormal_lcdf	0.63	1.06	1.44	0.97
exp_mod_normal_lcdf	0.75	1.62	2.77	1.38
exp_mod_normal_lpdf	0.53	0.92	1.03	0.98
skew_normal_lpdf	0.51	0.91	1.10	1.04

What the columns show:

• As opened, the lcdfs were faster, because develop's std_normal_lcdf is the slow Taylor/A&S code. Phi and the cdfs were slower, because exp(log Phi) adds a log and an exp to every call where develop calls only erf or erfc. exp_mod_normal_lcdf was slower too: it evaluates two log cdfs and a log_diff_exp, but develop returns NaN or infinite values in large parts of the parameter space.
• Before my commits, the double times are the same as at opening, but 9a8deff made the var paths up to 1.7 times slower (exp_mod_normal_lcdf 1.62 → 2.77, lognormal_cdf 1.26 → 1.72). My previous comment explains the cause and the fix (f4cd10e, c30ee6f).
• Now every function is faster than develop with double, except exp_mod_normal_lcdf. With var the same holds, except that exp_mod_normal_lpdf and skew_normal_lpdf are at develop's time (0.98 and 1.04), not at the 0.92 and 0.91 they had when the PR was opened.

For N = 100 and 10000 the "now" values are within about 0.1 of these; at N = 100000 some var values move toward develop's (std_normal_cdf 0.91, normal_lcdf 0.55). With z on [-30, 30], develop's Taylor code is slower, and std_normal_lcdf and normal_lcdf now take 0.33 of develop's time with double and 0.26 to 0.28 with var, as in the description. On a Xeon E5-2680 v3 with g++ 14.3, the PR as opened showed the same pattern (for example Phi 2.2 and std_normal_cdf 1.3 times develop's time with double).

Why the description shows larger gains

The absolute times show where the difference is. ns per element, double, N = 1000:

function	develop, description	develop, EPYC	PR as opened, description	PR as opened, EPYC	now, EPYC
Phi	27.2	8.9	16.5	20.3	6.5
std_normal_cdf	26.8	9.3	14.1	13.4	7.6
normal_cdf	28.7	10.1	13.4	13.7	7.8
lognormal_cdf	33.0	17.2	14.7	17.7	13.2
std_normal_lcdf	38.0	14.4	14.0	13.5	11.5
normal_lcdf	34.5	15.4	13.4	13.7	11.9
lognormal_lcdf	36.8	21.6	14.1	17.7	15.9
exp_mod_normal_lcdf	57.2	24.3	37.6	39.3	39.3
exp_mod_normal_lpdf	29.8	13.2	12.6	12.9	10.9
skew_normal_lpdf	30.7	13.3	12.4	13.9	11.7

The PR as opened ran at about the same speed on both machines (0.8 to 1.1 times the EPYC time), but develop took 1.7 to 3.1 times as long on the machine of the description. The largest factors (2.8 to 3.1) are for Phi, std_normal_cdf and normal_cdf, whose develop code is little more than an erf or erfc call per element. That fits a C library whose erf and erfc are much slower than glibc's. So the description measured a slow develop, not a faster branch. @andrjohns, which OS, compiler and CPU did you measure on?

On such a platform the erfc forms of Phi and the cdfs from bf5861a and 0d4ac36 would also be slow. Commit 9b38ecb therefore computes them without the C library's erf and erfc. With u = x / sqrt(2):

• for |u| < 0.46875, Phi is 0.5 + 0.5 * erf(u) with the Cody (1969) rational for erf that this PR used before 5f33a10 (no exp; Phi is near 1/2, so nothing cancels);
• otherwise it is 0.5 * exp(-(0.5 * x) * x) * erfcx(|u|) for x < 0, and one minus that for x > 0, with the erfcx rationals from erfcx.hpp.

Only exp comes from the C library. internal::std_normal_cdf_value_ratio, which std_normal_cdf, normal_cdf and lognormal_cdf use, computes the same values (bitwise equal to Phi) and the ratio phi / Phi as sqrt(2 / pi) / erfcx(-u) for x < 0 and in the central range; that needs no exp and has no 0 / 0 where Phi underflows. The OpenCL Phi device function uses the same formula, so the Phi and inv_Phi registrations include erfcx_device_function again. On glibc this form is also faster than develop's erf/erfc code, as the "now" column shows; on two Xeons (E5-2680 v3 and Gold 6248) Phi(double) takes 0.66 to 0.86 of develop's time over the input ranges I tried.

The new form is more accurate in the lower tail. Worst relative error of Phi in ulp against mpmath, 47001 points on [-38.5, 8.5]:

x	develop	original PR (exp(lcdf))	erfc (bf5861a)	erfcx (9b38ecb)
[-38.5, -30)	(cut off at -37.5)	1822	1704	515
[-30, -20)	1224	922	1224	258
[-20, -10)	470	440	470	130
[-10, -5)	118	93	118	34
[-5, -1)	766133	19.1	30.5	12.5
[-1, 0)	3.2	4.6	3.0	3.9
[0, 8.5)	1.3	2.6	1.3	1.4

In the tail the only error that grows with |x| is now one rounding of x * x; the erfc forms also round the argument x / sqrt(2), and erfc magnifies that. The ratio of the cdf helper improves more: on [-30, -1) from up to 1420 ulp to 4 to 7 ulp, and below -37.5, where Phi is subnormal and the erfc helper returned a wrong ratio, it is now within 4 ulp. std_normal_cdf follows Phi (worst value error 1500 ulp before, 2.9 now); normal_cdf and lognormal_cdf improve less, because the rounding of (y - mu) / sigma dominates there.

[avehtari]

9a8deff also made skew_normal_lpdf slower with var arguments, for the same reason as the lcdfs in my previous comment: z and score became lazy expressions, and the partials and the sum read z up to five times and score up to three times. New commit 927099e evaluates both once in the autodiff branch. The results are bitwise identical.

A first version also evaluated z in the double path. That made double calls 7 % slower, because there z is read only twice and the allocation costs more than the second read. So the commit keeps z lazy in that path.

Time relative to develop (develop = 1.00, below 1.00 is faster), Xeon Gold 6248, g++ 14.3, vectors of N = 1000, z uniform on [-3, 3], 2 runs × 3 rounds:

skew_normal_lpdf	as opened (a28c1c1)	before (9b38ecb)	after (927099e)
double	0.97	0.87	0.86
var	0.84	0.96	0.85

I also timed the OpenCL device functions on a Tesla V100, with each work item walking through the input range in registers and a baseline kernel subtracted. The device Phi from 9b38ecb is 2 to 4 times faster than the exp(lcdf) form of the original PR, and the device std_normal_lcdf_impl and its derivative are as fast as or faster than then. In [-5, 8] the device Phi takes up to 1.8 times as long as develop's, which uses erf there and loses up to 9e5 ulp on [-5, -1); the built-in erfc was slower than the current form on the V100. The difference is about 6 ps per evaluation, about 1 % of the cost per element of a Stan OpenCL call, so I left the device code as it is.

exp_mod_normal_lpdf has the same lazy temporaries, but evaluating them made it 3 % slower with var on the same machine, and it is already at its time from when the PR was opened (0.87), so I did not change it.

exp_mod_normal_lcdf, _lccdf and _cdf

New commit a17bba2 replaces the computation of these three functions on the CPU. The previous form, log F = log_diff_exp(log Phi(z), b) with b = v^2/2 - v z + log Phi(z - v) and v = lambda * sigma, has two failures at this PR's head (grid of z in [-1e4, 100] and v in [1e-8, 100], against mpmath):

function	region	value	dy	dsigma	dlambda
lcdf	left tail (z <= -30)	2.9e5 ulp	3.4e7	3.4e7	1.5e13
lcdf	small v (<= 1e-4)	non-finite (10 points)
lccdf	small v	1.4e12	1.3e12	1.3e12	6.0e2

In the left tail both terms are of size z^2/2; for small v they are nearly equal. At (z, v) = (-1e4, 0.01) the lcdf value is about -5e7 with an absolute error of 1.6e-3.

With u = -z / sqrt(2), t = u + v / sqrt(2), E = erfcx and g = exp(-z^2/2) / 2, Phi(z) = g E(u) and exp(v^2/2 - v z) Phi(z - v) = g E(t), so F = g (E(u) - E(t)) and the z^2/2 terms cancel exactly. All partials are ratios of E(u), E(t), E'(t) and E(u) - E(t). The difference E(u) - E(t) is computed directly when E(t) < 0.99 E(u), and otherwise from a series in the derivatives of erfcx at the midpoint, or for arguments of at least 4 from the difference of the Cody tail rational; the rounded t - u is never used. For z > 0 the code uses F = 1 - (Phi(-z) + b) where F >= 1/4, and a form with two positive terms where F < 1/4. The three functions share one scalar loop.

Worst error over the value and the four partials after the change, no non-finite result: 510 ulp (lcdf) and 480 ulp (lccdf) on the grid above, and 640 and 900 ulp on a denser grid of 2500 points with |z| in [1e-4, 100] and v in [1e-10, 100]. One exception: dsigma near z = 1/v at large v (3.9e3 ulp at (0.0316, 31.6)), where this partial has a true zero.

The change also makes the functions faster. Time relative to develop, AMD EPYC 9965, vectors of N = 1000, bulk: z on [-3, 3], wide: z on [-30, 30]:

function	inputs	before double	after double	before var	after var
lcdf	bulk	1.43	0.80	1.37	0.94
lcdf	wide	2.06	0.94	1.76	1.09
lccdf	bulk	1.47	0.88	1.38	0.99
lccdf	wide	2.17	0.79	1.91	0.89
cdf	bulk	1.68	0.93	1.53	1.00
cdf	wide	2.52	1.15	2.03	1.15

A Xeon Gold 6230 gives the same pattern. The cdf is exp(lcdf); develop computes it from erf directly, which is not accurate in the tails.

The OpenCL code still uses the previous form. A port of the new algorithm to the device is in testing, and I will add it as a separate commit.

[avehtari]

OpenCL port of the exp_mod_normal cdfs

As announced in my previous comment, new commit 39886e0 moves the OpenCL versions of exp_mod_normal_lcdf, _lccdf and _cdf to the algorithm of a17bba2. The new device function exp_mod_normal_lcdf_part(z, v, upper, which) follows the CPU code branch for branewnch. The kernel generator only has functions with one output, so eight binary wrappers (exp_mod_normal_{lcdf,lccdf}_{value,d_z,d_v,d_s}) select the value or one of the three partials, and opencl/prim/exp_mod_normal_lcdf.hpp uses them. Compiled on the host, the device code agrees with the CPU code to 5e-14 relative on the 2500-point grid (1.2e-12 at the dsigma zero); the difference comes from the device erfcx, which evaluates its polynomials in a different order.

Time relative to develop, Tesla V100, matrix_cl (double, value only) and var_value<matrix_cl> with var parameters (value and gradient), N = 1e6; the run-to-run spread was 10 to 60 %:

function	inputs	before double	after double	before var	after var
lcdf	bulk	0.86	0.79	0.80	0.89
lcdf	wide	0.58	0.52	0.64	0.85
lccdf	bulk	1.13	0.99	1.09	1.21
lccdf	wide	1.23	1.14	1.07	1.39
cdf	bulk	1.19	1.09	1.05	1.19
cdf	wide	1.28	1.20	1.10	1.51

The var path is slower than before because each of the four wrappers runs the whole device function, and the device compiler does not merge the four copies. The extra cost is 1 to 2 ns per element. A dedicated kernel that computes all four outputs in one pass would remove it, but it would not follow the pattern of the other distributions, so I did not write one.

[avehtari]

@SteveBronder I'm now done improving this PR

[SteveBronder]

Thanks! I'm running some benchmarks locally to see what is going on in more detail

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.22	0.22	1.0	0.28% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	0.99	-0.93% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.31	6.26	1.01	0.75% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.96	-4.59% slower
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.29	7.04	1.18	15.06% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.46	4.45	1.0	0.26% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	217.62	218.53	1.0	-0.42% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.9	0.9	1.0	-0.49% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.62	0.61	1.01	0.52% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	42.88	42.7	1.0	0.42% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.58	0.59	1.0	-0.43% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	20.95	20.88	1.0	0.36% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	0.16% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.17	3.19	0.99	-0.7% slower
performance.compilation	366.95	362.77	1.01	1.14% faster
Mean result: 1.0095524481511577

Jenkins Console Log
Jenkins Build Stages
Commit hash: a2236819d562837e97a9bfe97dc9ebd6d8613729

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           43 bits physical, 48 bits virtual
CPU(s):                                  256
On-line CPU(s) list:                     0-255
Thread(s) per core:                      2
Core(s) per socket:                      64
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              23
Model:                                   49
Model name:                              AMD EPYC 7742 64-Core Processor
Stepping:                                0
Frequency boost:                         enabled
CPU MHz:                                 1497.316
CPU max MHz:                             3416.0681
CPU min MHz:                             1500.0000
BogoMIPS:                                4491.85
Virtualization:                          AMD-V
L1d cache:                               4 MiB
L1i cache:                               4 MiB
L2 cache:                                64 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-63,128-191
NUMA node1 CPU(s):                       64-127,192-255
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Mitigation; untrained return thunk; SMT enabled with STIBP protection
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Retpolines; IBPB conditional; STIBP always-on; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Not affected
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba ibrs ibpb stibp vmmcall fsgsbase bmi1 avx2 smep bmi2 cqm rdt_a rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif v_spec_ctrl umip rdpid overflow_recov succor smca sev sev_es

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

[avehtari]

@SteveBronder. I tested your commits on a cluster (AMD EPYC 9965, g++ 11.5, one benchmark binary per function). exp_mod_normal_cdf is 1.2 to 1.3 times faster after e1fea06, and after cf177ef std_normal_lcdf, normal_lcdf, lognormal_lcdf, exp_mod_normal_lpdf and skew_normal_lpdf are 1.2 to 1.4 times faster for z on [-30, 30]. Other functions change by at most 6 %. Accuracy against mpmath stays within about 2.7 x 2^-52 for erfcx and std_normal_lcdf.

1. I found one new edge case. With inv(sigma), a sigma below 1 / DBL_MAX (about 5.6e-309) makes 1 / sigma infinite, and at y == mu the value becomes 0 * inf = NaN in normal_lcdf, normal_lccdf, lognormal_lcdf and lognormal_lccdf. Before cf177ef these returned correctly log(0.5). I pushed d081fff, which keeps your multiplication and recomputes z with a division only when some inv(sigma) is infinite. The benchmark shows no change in speed (0.99 to 1.01 for double and var, default and -march=native), and the commit adds a test for it.

A note for comparing benchmark results: when I first built all 12 functions into one binary, skew_normal_lpdf looked 20 % slower after e1fea06, although that commit does not touch it, and Phi looked 28 % slower with -march=native after cf177ef. With one binary per function both effects were gone, so they came from GCC's inlining decisions for the whole translation unit.

2. I also pushed 1f8b6ba, which changes exp_mod_normal_lcdf_from_delta and its OpenCL port. Far in the lower tail (for example z = -1e5 with lambda * sigma = 1e-300, or z below -1e150) the partials were formed as e_t * (1 / Delta), where 1 / Delta overflows although the product is finite, so the gradient was infinite while the value was correct. At z = -1e300 with v = 1e10, exp_mod_normal_erfcx_h_diff also gave NaN values, because delta * (y + x) overflowed before the multiplication by w0 = 0. On a grid of 105 points from z = -1e300 to 10 and v = 1e-300 to 1e10, non-finite results drop from 38 to 0 for the cdf, from 1 to 0 for the lccdf, and from 38 to 28 for the lcdf, where all 28 are at |z| >= 1e155 and the value -inf and the sigma partial overflow (develop has 89, 89 and 56). The commit adds six mpmath rows for the deep tail to the CPU and OpenCL reference tests. The exp_mod_normal CPU suites and generated tests pass, and the OpenCL suites pass on a V100; with the old kernels the new reference tests fail. The cost is one range check in the hot path: double is unchanged, and var lcdf and cdf are 1 to 3 % slower with default flags and 1 to 6 % slower with -march=native.

[stevieb-bot]

Fact check, part 1 of 2.

Action: Not Ready for Human Review

Commit checked: 150f54df16e5

Summary

Thank you for this PR and for the detailed write-ups in the comments.

The PR moves the normal family of log CDFs onto one shared std_normal_lcdf kernel built on the erfcx helpers, and recomputes Phi, the cdfs and the exp_mod_normal lcdf, lccdf and cdf from erfcx.
The description and comments claim wider valid input ranges, better accuracy against mpmath, large speedups, and a fix for #2881.

I read the code and the tests at 150f54d.
In a container without network access I ran the unit suites the comments list at the PR head, the generated test/prob tests for the cdfs and for exp_mod_normal, and the OpenCL suites on an NVIDIA RTX A6000 (gpu_device-head.out records the device).
I also copied the PR's new tests into the merge base with develop (e4745c8) and ran them there to see which ones fail on develop.

The output of each run is in this gist.
Its README.md gives the command and commit behind each file.

• Claims and Evidence: 7 Confirmed, 12 Partially Confirmed, 25 Not Verified.
• Description vs Code: 18 Confirmed, 4 Partially Confirmed.

Claims and Evidence

Each section below lists facts from the PR summary and discussion that an agent was able to confirm, partially confirm, or not confirm at all, or that a run of the code in the PR contradicted.
Please review the "Partially Confirmed", "Not Verified", and "Contradicted" facts for facts missing evidence found by this agent.

Confirmed Facts

FACT: ordered_probit_lpmf and von_mises_cdf call the changed functions and have no OpenCL version

Quote:

Two functions that benefit from this std_normal_lcdf PR that are missing from the PR text are

(comment by avehtari)

STATUS: Confirmed

ordered_probit_lpmf calls std_normal_lcdf in four places, and von_mises_cdf returns normal_cdf(z, mu, sigma) in its normal approximation branch.
stan/math/opencl/prim/ has neither function.
Both therefore use the CPU kernels that this PR changes.

Evidence:

The claim follows from the code.
I read the two files at the reviewed commit.
Code: ordered_probit_lpmf.hpp#L91-L97, von_mises_cdf.hpp#L63.

What I read: ordered_probit_lpmf.hpp lines 91, 93, 96 and 97 call std_normal_lcdf. von_mises_cdf.hpp line 63 calls normal_cdf.
Neither file exists under stan/math/opencl/prim/.

---

FACT: wiener4_lcdf_unnorm line 128 subtracts std_normal_lpdf from std_normal_lcdf

Quote:

wiener4_lcdf_unnorm line 128 is:

(comment by avehtari)

STATUS: Confirmed

Line 128 of wiener4_lcdf_unnorm.hpp at the reviewed commit is return std_normal_lcdf(-x) - std_normal_lpdf(x);, as the comment quotes.

Evidence:

The claim follows from the code.
I read the line at the reviewed commit.
Code: wiener4_lcdf_unnorm.hpp#L128.

What I read: wiener4_lcdf_unnorm.hpp line 128 matches the quoted line.

---

FACT: 5f33a10 leaves the double code of the erfcx helpers unchanged, so erfcx returns the same bits

Quote:

I did not change their double code, so erfcx returns the same bits as before.

(comment by avehtari)

STATUS: Confirmed

In 5f33a10 the diff of erfcx.hpp only turns erfcx_cody_tail, erfcx_cody_middle and erfcx_small into templates, replacing double with T, and adds @tparam lines.
The arithmetic is unchanged, so the double instantiations are the same code as before.
A later commit, cf177ef, changes erfcx_cody_tail to put the tail rational over one denominator, At the reviewed commit, erfcx for arguments of 4 and above therefore no longer follows the develop formula.
The 'same bits' statement holds for 5f33a10.
It may not hold at the reviewed commit.

Evidence:

The claim follows from the code.
I read the diff of 5f33a10 and the current erfcx.hpp.
Code: 5f33a10d31 diff, erfcx.hpp#L99-L122 at the reviewed commit.

What I read: git show 5f33a10d31 -- stan/math/prim/fun/erfcx.hpp changes only the types and adds template headers.
At the reviewed commit, erfcx.hpp lines 99 to 109 return the numerator INV_SQRT_PI * s_q + pq.first and the denominator x * s_q, and erfcx_cody_tail (lines 118 to 122) divides them.
Commit cf177ef introduced this form.

---

FACT: 5f33a10 touches 5 files, with 59 lines added and 153 deleted

Quote:

The diff touches 5 files: 59 lines added, 153 deleted.

(comment by avehtari)

STATUS: Confirmed

git show --stat 5f33a10d31 lists 5 files with 59 insertions and 153 deletions.

Evidence:

The claim follows from the commit.
Code: 5f33a10d31.

What I read: git show --stat 5f33a10d31 prints 5 files changed, 59 insertions(+), 153 deletions(-).

---

FACT: After 5f33a10 the 20 CPU suites and 11 OpenCL suites pass

Quote:

CPU, 20 suites, all pass: [...]
OpenCL on a Tesla V100, 11 suites, all pass:

(comment by avehtari)

STATUS: Confirmed

The comment lists the suites it ran at 686fc67 plus 5f33a10.
I ran the same suites at the reviewed commit, which is later and includes all of avehtari's and SteveBronder's commits.
The 50 CPU test files all pass.
They cover the listed std_normal_lcdf_impl, erfcx, Phi, and inv_Phi suites in prim, rev, fwd, and mix, and the normal, std_normal, lognormal, exp_mod_normal, and skew_normal suites.
The 18 OpenCL test programs include the 11 named here.
All of them pass on an NVIDIA RTX A6000.
The comment used a Tesla V100.

Evidence:

The suites are in the PR.
The clang-format and cpplint part of the comment was not rerun.
Code: test/unit/math/mix/prob/std_normal_lcdf_impl_test.cpp, test/unit/math/opencl/rev/std_normal_lcdf_test.cpp.

What I ran: python3 runTests.py -j8 on 50 CPU test files at the head exits 0 with 160 tests passing.
python3 runTests.py -j8 on 18 OpenCL test files with STAN_OPENCL=true exits 0 with 160 tests passing.
Output in the gist: cpu_suites-head.out, opencl-head.out.

---

FACT: After bf5861a the 22 CPU suites and the OpenCL suites, including Phi_test and inv_Phi_test, pass

Quote:

CPU, 22 suites, all pass: the 20 listed in my previous comment, plus prim and mix inv_Phi_test.

(comment by avehtari)

STATUS: Confirmed

At the reviewed commit the CPU run includes prim and mix Phi_test and inv_Phi_test with the 20 suites from the earlier comment, and all pass.
MathFunctions.Phi with the 1e-14 checks at -5 and -3, and PhiTails, pass.
The OpenCL run includes opencl/rev Phi_test and inv_Phi_test, and all pass.
The reviewed commit has the erfcx form of Phi from 9b38ecb, which replaced the form this comment tested.

Evidence:

The suites are in the PR.
Code: test/unit/math/prim/fun/Phi_test.cpp, test/unit/math/opencl/rev/Phi_test.cpp.

What I ran: The same two runs as above.
Both exit 0.
Output in the gist: cpu_suites-head.out, opencl-head.out.

---

FACT: After 1f8b6ba the exp_mod_normal CPU suites, generated tests and OpenCL suites pass

Quote:

The exp_mod_normal CPU suites and generated tests pass, and the OpenCL suites pass on a V100

(comment by avehtari)

STATUS: Confirmed

At the reviewed commit the exp_mod_normal unit tests in prim, rev and mix pass in the CPU run.
runTests.py -j8 test/prob/exp_mod_normal builds and runs the 156 generated test programs and exits 0.
The OpenCL run includes exp_mod_normal_lcdf_test, exp_mod_normal_lccdf_test, exp_mod_normal_cdf_test, exp_mod_normal_lpdf_test and exp_mod_normal_lcdf_reference_test, and all pass on an RTX A6000.

Evidence:

The suites are in the PR or the repository.
Code: test/unit/math/rev/prob/exp_mod_normal_lcdf_test.cpp, test/unit/math/opencl/rev/exp_mod_normal_lcdf_reference_test.cpp.

What I ran: python3 runTests.py -j8 test/prob/exp_mod_normal at the head exits 0 after 4178 s.
The CPU and OpenCL runs above also exit 0.
Output in the gist: prob_exp_mod_normal-head.out, cpu_suites-head.out, opencl-head.out.

---

Partially Confirmed Facts

FACT: The original PR's cdfs take 1.4 to 1.9 times as long as develop, and 0d4ac36 restores develop's speed

Quote:

The original PR also computes these three as exp(lcdf), and with double arguments they take 1.4 to 1.9 times as long as on develop.

(comment by avehtari)

STATUS: Partially Confirmed

The first part of the quote follows from the code.
At a28c1c1, the PR as opened, std_normal_cdf, normal_cdf and lognormal_cdf each return exp(...) of the matching lcdf impl.
The comment's speed table for double and var arguments comes from a benchmark that the PR does not include.

Evidence:

The code at a28c1c1.
The benchmark source is missing.
Code: normal_cdf.hpp#L33 at a28c1c1637, std_normal_cdf.hpp#L28 at a28c1c1637, lognormal_cdf.hpp#L16 at a28c1c1637.

What I read: At a28c1c1, normal_cdf.hpp line 33, std_normal_cdf.hpp line 28 and lognormal_cdf.hpp line 16 return exp(internal::..._lcdf_impl<false>(...)).

Request: Could you add the benchmark source for this table?

---

FACT: 9a8deff slowed down the var lcdf paths, and c30ee6f recovers the time

Quote:

One commit is responsible: 9a8deff ("Refactor normal log-CDF dispatch and density expressions") replaced evaluated temporaries with lazy expressions.

(comment by avehtari)

STATUS: Partially Confirmed

The code part of the claim follows from the diffs.
Commit 9a8deff replaces const auto& z = to_ref(...) and to_ref_if<...>(slopes / sigma_val) in normal_lcdf.hpp with plain auto z = ... and auto scaled_slope = ....
It does the same in lognormal_lcdf.hpp and for diff, z, v, log_b, lp, s_b and q in exp_mod_normal_lcdf.hpp.
Commit c30ee6f wraps these in eval(...) again.
The claim that this commit caused the slowdown, and the ns per element in the table, come from a benchmark that the PR does not include.

Evidence:

The diffs of the two commits.
The benchmark that timed each commit is missing.
Code: 9a8deffbef, c30ee6f9ef.

What I read: git show 9a8deffbef removes the to_ref and to_ref_if evaluations of z and scaled_slope in normal_lcdf.hpp and lognormal_lcdf.hpp and of the temporaries in exp_mod_normal_lcdf.hpp.
git show c30ee6f9ef adds eval(...) around them.

Request: Could you add the benchmark used to time each commit, for the timing part of this claim?

---

FACT: OpenCL exp_mod_normal timings before and after 39886e0

Quote:

The var path is slower than before because each of the four wrappers runs the whole device function, and the device compiler does not merge the four copies.

(comment by avehtari)

STATUS: Partially Confirmed

The first part of the quoted sentence follows from the code.
Each of the eight wrappers calls the full exp_mod_normal_lcdf_part(z, v, upper, which) and selects one output.
The PR does not include the benchmark needed to check that the device compiler does not merge the four copies, or the timings in the table (Tesla V100, N = 1e6, run-to-run spread of 10 to 60 %).

Evidence:

The wrapper code.
The benchmark source is missing.
Code: device_functions/exp_mod_normal_lcdf.hpp#L374-L395.

What I read: Lines 374 to 395 define the eight wrappers, and each one returns exp_mod_normal_lcdf_part(z, v, upper, which) with a different which.

Request: Could you add the benchmark source for this table?

---

FACT: The PR fixes #2881 (non-finite gradients of normal_lcdf and normal_lccdf at infinite y)

Quote:

Fixes #2881.

(PR description)

STATUS: Partially Confirmed

#2881 shows normal_lcdf(+inf | mu, sigma) with mu and sigma as parameters, where develop returns a NaN sigma adjoint.
The issue title also names normal_lccdf.
The PR adds mathMixScalFun_normal_lcdf_infinite_endpoints, which checks the value at both infinite endpoints and that the sigma adjoint of normal_lcdf(INFTY, 0.0, sigma) is 0.
On the merge base with develop this test fails with sigma.adj() equal to -nan, and at the head it passes.
mathMixScalFun_lognormal_lcdf_tail_and_endpoints checks the same for lognormal_lcdf and also fails on develop with a NaN sigma adjoint.
The fix is the select(scaled_slope == 0, 0.0, -scaled_slope * z) partial in normal_lcdf_impl, which normal_lccdf shares.
The PR does not test normal_lccdf at y = -inf with a var sigma, or the mu adjoint from the issue's example.

Evidence:

The PR adds a test for the normal_lcdf case of #2881.
A test for normal_lccdf, the other function in the issue title, is missing.
Code: mix/prob/normal_cdf_log_test.cpp#L56-L67, mix/prob/lognormal_cdf_log_test.cpp#L20-L36, normal_lcdf.hpp#L63-L67.

What I read: normal_lcdf.hpp lines 63 to 67 set the sigma partial to 0 where the scaled slope is 0, which covers y = +inf for the lcdf and y = -inf for the lccdf.

What I ran: python3 runTests.py -j8 test/unit/math/mix/prob/normal_cdf_log_test.cpp on the merge base with the PR's tests copied in.
mathMixScalFun_normal_lcdf_infinite_endpoints fails with sigma.adj() = -nan.
The same test passes at the head in the combined CPU run.
The lognormal_lcdf endpoint test also fails on the base with a NaN sigma adjoint, and after that failure the binary aborts on an Eigen size assertion in the same test.
It passes at the head.
Output in the gist: normal_cdf_log_test-base.out, lognormal_cdf_log_test-base.out, cpu_suites-head.out.

Request: Could you add the normal_lccdf case from the issue title, for example normal_lccdf(NEGATIVE_INFTY, mu, sigma) with var mu and sigma, and also check mu.adj() as in the issue's example?

---

FACT: The valid input ranges before over- or underflow are wider than on develop

Quote:

The valid ranges of inputs (before over/underflow) across all functions are expanded:

(PR description, Summary)

STATUS: Partially Confirmed

The table gives develop and branch ranges for the value and the gradient of 15 functions.
The PR does not include the script that found the range ends, but its tests check some of the cells:

• std_normal_lcdf(-1.5e154) is finite (large_finite_log_density), and the kernel matches a reference at -1e150 (scalar_tail_kernel).
  Both pass at the head.
  They are in a file that also calls the new internal kernel, which develop does not have, so the file does not compile on develop.
  Both points lie past the table's branch end of -1e154, so they show that the range reaches at least that far.
• lognormal_lcdf(exp(-50), 0, 1) returns -inf on develop and the reference value -1254.83 at the head, consistent with the develop range ending near -38.5 for log(y).
• exp_mod_normal_lcdf(-10, 0, 1, 1) returns -inf on develop and -55.65 at the head, which fits the develop lower end of -7.2.
  At y = 60, develop returns 0 where the reference is -1.44e-26, which fits the develop upper end of 37.9.
  The deep-tail reference rows from z = -1e5 (y = -199999.5, the table's branch end) down to z = -1e154 pass at the head and return -inf on develop.
• The OpenCL test empty_and_extreme_gradient checks that the std_normal_lcdf gradient at -1.5e308 is 1.5e308, and it passes on the GPU.
  This point also lies past the table's open bound at -1e308.

The PR does not test the other cells, for example the develop and branch ends of normal_cdf, Phi, skew_normal_lpdf and the branch upper ends of exp_mod_normal_lcdf near 745.

Evidence:

Tests at single points for some cells.
The script that located the range ends is missing.
Code: prim/prob/std_normal_lcdf_impl_test.cpp#L6-L43, mix/prob/lognormal_cdf_log_test.cpp#L20-L36, mix/prob/exp_mod_normal_cdf_log_test.cpp#L26-L38, rev/prob/exp_mod_normal_lcdf_reference.hpp#L83-L104, opencl/rev/std_normal_lcdf_test.cpp#L63-L74.

What I ran: The CPU and OpenCL runs at the head pass all of these tests.
On the merge base with the PR's tests copied in, mathMixScalFun_lognormal_lcdf_tail_and_endpoints, mathMixScalFun_exp_mod_normal_lcdf_tails_and_endpoints and ProbDistributionsExpModNormal_deep_tail_reference fail with -inf values.
Output in the gist: cpu_suites-head.out, opencl-head.out, lognormal_cdf_log_test-base.out, exp_mod_normal_cdf_log_test-base.out, exp_mod_normal_lcdf_test-base.out.

Request: Could you add the script that found the range ends in the table?
A test at both ends of each branch range (just inside, finite and accurate) would also let the table be checked at the current head, since several functions changed after the table was written.

---

FACT: Develop's worst gradient error is 7.6e-06 to 5.6e-05, and the PR improves the gradient by about ten orders of magnitude

Quote:

For comparison, develop's worst in-range relative error is 7.6e-06 to 5.6e-05, which is 3.4e+10 to 2.5e+11 ulp, across every branch except the far Mills tail.
So this PR improves the gradient by about ten orders of magnitude.

(comment by avehtari)

STATUS: Partially Confirmed

The PR's expect_branch_accuracy checks the gradient at 15 points in the five former Taylor branches against 60-digit references, now at a relative tolerance of 1e-12.
On the merge base with develop all 15 points fail, with relative errors from 8.4e-12 to 6.09e-05.
The worst error per branch is 3.26e-05, 2.80e-05, 2.30e-05, 6.09e-05 and 7.61e-06, so the largest is somewhat above the 5.6e-05 in the comment.
At the head all 15 points pass at 1e-12.
At the worst point (6.09e-05 on develop) the error therefore drops by at least seven orders of magnitude, and at the other points the bound from the test is smaller.
The test cannot show the full ten orders.
The comment's figure for develop comes from a grid whose script is not in the PR.

Evidence:

A test at 15 points that fails on develop and passes at the head.
The mpmath grid script is missing.
Code: normal_lcdf_tail_test_helpers.hpp#L156-L195.

What I ran: python3 runTests.py -j8 test/unit/math/mix/prob/normal_cdf_log_test.cpp on the merge base with the PR's tests copied in.
mathMixScalFun_normal_lcdf_branch_accuracy fails at all 15 points.
normal_ccdf_log_test and std_normal_ccdf_log_test fail their branch_accuracy tests the same way.
All three pass at the head.
Output in the gist: normal_cdf_log_test-base.out, normal_ccdf_log_test-base.out, std_normal_ccdf_log_test-base.out, cpu_suites-head.out.

Request: Could you add the grid script behind the 7.6e-06 to 5.6e-05 range?
The PR's own test reaches 6.09e-05 at y = 2.1213, so the upper end may need updating.

---

FACT: Develop's 0.5 * (1 + erf(x / sqrt(2))) branch of Phi on [-5, -0.67] loses up to 7.7e5 ulp

Quote:

The change was made for accuracy, and develop does have a defect: its 0.5 * (1 + erf(x / sqrt(2))) branch on [-5, -0.67] loses up to 7.7e5 ulp.

(comment by avehtari)

STATUS: Partially Confirmed

The PR's MathFunctions.Phi test checks Phi(-5) and Phi(-3) against reference values at a relative tolerance of 1e-14.
On the merge base with develop the test fails at -5 with a relative error of 3.9e-11, and passes at -3.
At the head both points pass.
So the test shows a large error in this branch at one point.
It does not measure the worst case of 7.7e5 ulp, which comes from a grid whose script is not in the PR.

Evidence:

A test at two points that fails on develop at one of them.
The grid script is missing.
Code: prim/fun/Phi_test.cpp#L5-L10.

What I ran: python3 runTests.py -j8 test/unit/math/prim/fun/Phi_test.cpp on the merge base with the PR's test copied in.
MathFunctions.Phi fails at Phi(-5.0) with a relative difference of 3.8966718740596207e-11.
It passes at the head.
Output in the gist: Phi_test-base.out, cpu_suites-head.out.

Request: Could you add the script for the 7.7e5 ulp figure, or a test point in [-5, -0.67] where develop's error is largest?

---

FACT: The previous exp_mod_normal form fails in the left tail and for small v

Quote:

The previous form, log F = log_diff_exp(log Phi(z), b) with b = v^2/2 - v z + log Phi(z - v) and v = lambda * sigma, has two failures at this PR's head

(comment by avehtari)

STATUS: Partially Confirmed

The comment describes failures of the log_diff_exp form that the PR used before a17bba2.
The PR's exp_mod_normal_lcdf_test.cpp checks 13 mpmath reference rows, including a left-tail row at z = -30, v = 1e-6, a far-tail row at z = -100 and a small-v row at z = -1, v = 1e-6, at a relative tolerance of 1e-12.
At the head all rows pass for the lcdf, lccdf and cdf.
On the merge base with develop every test in the file fails.
For example, the lcdf returns -inf at z = -30, v = 1e-6 and has a relative error of about 2e-11 at z = -1, v = 1e-6.
That run tests develop's form.
The PR's form before a17bba2, which the comment is about, was not built, and the ulp figures in the table come from a grid whose script is not in the PR.

Evidence:

mpmath reference rows that pass at the head and fail on develop.
The grid script, and a run of the form before a17bba2, are missing.
Code: rev/prob/exp_mod_normal_lcdf_reference.hpp, rev/prob/exp_mod_normal_lcdf_test.cpp#L113-L138.

What I ran: python3 runTests.py -j8 test/unit/math/rev/prob/exp_mod_normal_lcdf_test.cpp on the merge base with the PR's test and reference header copied in.
All seven tests fail.
They pass at the head.
Output in the gist: exp_mod_normal_lcdf_test-base.out, cpu_suites-head.out.

Request: Could you add the grid script behind the failure table, and say which reference rows fail with the form before a17bba2?

---

FACT: With the old kernels the new deep-tail reference tests fail

Quote:

with the old kernels the new reference tests fail

(comment by avehtari)

STATUS: Partially Confirmed

ProbDistributionsExpModNormal_deep_tail_reference checks the six deep-tail rows for the lcdf, lccdf and cdf at a relative tolerance of 1e-12.
On the merge base with develop it fails, with -inf values at z = -1e5.
At the head it passes.
The comment refers to the PR's kernels before 1f8b6ba, which were not built, so the run shows only that develop's kernel fails these rows.

Evidence:

A reference test that fails on develop and passes at the head.
A run of the kernel before 1f8b6ba is missing.
Code: rev/prob/exp_mod_normal_lcdf_test.cpp#L133-L138.

What I ran: python3 runTests.py -j8 test/unit/math/rev/prob/exp_mod_normal_lcdf_test.cpp on the merge base with the PR's test copied in.
ProbDistributionsExpModNormal_deep_tail_reference fails.
It passes at the head.
Output in the gist: exp_mod_normal_lcdf_test-base.out, cpu_suites-head.out.

Request: Could you say which deep-tail rows failed with the kernel before 1f8b6ba, so the claim can be tied to specific rows?

---

FACT: A sigma below 1 / DBL_MAX made normal_lcdf and lognormal_lcdf return NaN at the median after cf177ef

Quote:

With inv(sigma), a sigma below 1 / DBL_MAX (about 5.6e-309) makes 1 / sigma infinite, and at y == mu the value becomes 0 * inf = NaN in normal_lcdf, normal_lccdf, lognormal_lcdf and lognormal_lccdf.

(comment by avehtari)

STATUS: Partially Confirmed

normal_lcdf_small_sigma_test.cpp checks that normal_lcdf, normal_lccdf, lognormal_lcdf and lognormal_lccdf return log(0.5) at the median for sigma down to 4.9e-324, with the expected slopes.
All three tests pass at the head.
The NaN described in the comment came from a PR commit (cf177ef), and develop divides by sigma.
On the merge base with develop, the value checks pass and the scalar normal test fails only on the slope tolerance (for example 0.79788456160 against 0.79788456080 at sigma = 1).
This reflects develop's slope accuracy and is unrelated to the NaN.
Reading cf177ef shows the mechanism.
normal_lcdf and lognormal_lcdf compute inv_sigma = to_ref(inv(sigma_val)) and multiply z by it, which is 0 * inf at the median when 1 / sigma overflows.
The GCC 13 failure that e076f6a reports for d081fff was not run, since that commit was not built.

Evidence:

A test that passes at the head, and the code at cf177ef.
A run of d081fff with GCC 13 is missing.
Code: rev/prob/normal_lcdf_small_sigma_test.cpp, normal_lcdf.hpp#L50-L52 at cf177ef96b, lognormal_lcdf.hpp#L54-L57 at cf177ef96b.

What I read: At cf177ef, normal_lcdf.hpp lines 50 to 52 and lognormal_lcdf.hpp lines 54 to 57 multiply y - mu (or log(y) - mu) by inv_sigma, which gives 0 * inf at the median when inv_sigma is infinite.

What I ran: python3 runTests.py -j8 test/unit/math/rev/prob/normal_lcdf_small_sigma_test.cpp passes at the head.
On the merge base, ProbNormalLcdf_median_with_subnormal_sigma fails on the y.adj() and mu.adj() tolerance, and the lognormal test and the vector normal test pass.
Output in the gist: cpu_suites-head.out, normal_lcdf_small_sigma_test-base.out.

Request: The fix itself does not need a change.
If the NaN matters for the record, could you note in the test which commit it guards against, since develop does not show it?

---

FACT: After 0d4ac36 the unit suites, the generated cdf tests, the OpenCL cdf tests and a compile check all pass

Quote:

The generated test/prob tests of the three cdfs (runTests.py test/prob/normal test/prob/lognormal test/prob/std_normal -f _cdf_0): 17 test programs, 12000 cases, all pass.
[...]
A compile check instantiates the three functions with fvar<fvar<var>>, std::vector<var>, Eigen vectors of fvar<double>, var_value column and row vectors, and mixed scalar and vector arguments.

(comment by avehtari)

STATUS: Partially Confirmed

At the reviewed commit the unit suites the comment lists pass in the CPU run.
The generated tests, run with the comment's command runTests.py test/prob/normal test/prob/lognormal test/prob/std_normal -f _cdf_0, build and pass 17 test programs.
review.py run cut the middle of the output because of its length, so I did not recount the 12000 cases.
The OpenCL std_normal_cdf_test, normal_cdf_test, lognormal_cdf_test and Phi_test pass.
The compile check that instantiates the three cdfs with fvar<fvar<var>>, std::vector<var> and the other types is not in the PR.

Evidence:

The unit, generated and OpenCL suites are in the PR or the repository.
The compile-check source is missing.
Code: test/prob/normal/normal_cdf_test.hpp, test/unit/math/mix/prob/lognormal_cdf_test.cpp.

What I ran: python3 runTests.py -j8 test/prob/normal test/prob/lognormal test/prob/std_normal -f _cdf_0 at the head exits 0 after running 17 test programs.
The CPU and OpenCL runs above also exit 0.
Output in the gist: prob_cdf-head.out, cpu_suites-head.out, opencl-head.out.

Request: Could you add the compile check (for example as a mix test that instantiates the three cdfs with the listed types), so it runs in CI?

---

FACT: 9a8deff passes its targeted tests and 256 propto checks, and the OpenCL tests compile but fail at kernel creation

Quote:

Validation: targeted value and nested-autodiff tests pass, along with 256 proportional-density checks.

(commit 9a8deffb)

STATUS: Partially Confirmed

9a8deff changed the prim and mix std_normal_lcdf_impl_test.cpp files, and at the reviewed commit both pass in the CPU run with the other value and nested-autodiff tests.
The OpenCL tests compile, run and pass on an RTX A6000, so the CL_INVALID_KERNEL failure from the commit message did not occur here.
The commit message does not say which tests make up the 256 proportional-density checks.
That count was not checked.
9a8deff itself was not built.

Evidence:

The tests are in the PR.
The list of the 256 propto checks is missing.
Code: test/unit/math/prim/prob/std_normal_lcdf_impl_test.cpp, test/unit/math/mix/prob/std_normal_lcdf_impl_test.cpp.

What I ran: The CPU and OpenCL runs above, both at the head, exit 0.
Output in the gist: cpu_suites-head.out, opencl-head.out.

Request: Could you say which tests the 256 proportional-density checks are, or drop the count from the record?

---

Not Verified Facts

FACT: The timing tables in the description show speedups of up to 5.08x

Quote:

While performance is either unchanged or significantly improved across inputs:

(PR description, Summary)

STATUS: Not Verified

The summary table gives the smallest and largest speedups, and the four full timing tables give ns per element for develop and the branch at N = 100 to 100000 for double, var, var_value<Eigen::Matrix> and fvar<double>.
The PR does not include the benchmark that produced them, and the description does not name the machine, compiler, flags, inputs or repetitions.
The opening sentence also says performance is "either unchanged or significantly improved", while the fvar table lists 0.92x for normal_cdf at N = 100000.
The comment by avehtari reports different ratios for the same functions on other machines (for example Phi at 2.27 times develop's time as opened).
Those comments do not include their benchmark code either, so this fact check cannot settle which numbers hold.

Evidence:

Tables only.
The benchmark source, the inputs, and the machine and compiler are missing.

Request: Could you add the benchmark source (for example a google-benchmark file or a script) that produced these tables, with the input distribution, the machine, the compiler and flags, and the number of repetitions?
It would also help to re-run it at the current head, since later commits replaced Phi and the cdf code that the tables time.

---

FACT: Worst-case ulp errors of the new kernel against mpmath, by range

Quote:

Worst ulp over 601 points per range:

(comment by avehtari)

STATUS: Not Verified

The table gives the worst ulp error of the value and the gradient over 601 points in each of nine ranges.
The PR does not include the script that evaluated these points against mpmath.
The PR's tests compare the kernel with high-precision references at single points, at a relative tolerance of 1e-12, so they bound the error at those points but do not reproduce the ulp table.

Evidence:

Table only.
The mpmath script and the grid are missing.

Request: Could you add the script that computed this table (the grid, the mpmath precision, and how the ulp error is formed), for example under a test/ or benchmarks/ folder or in a gist linked from the description?

---

FACT: Calling erfcx directly is faster than the internal std_normal_erfcx path

Quote:

Comparing your internal std_normal_erfcx path against calling erfcx directly, value and gradient computed together:

(comment by avehtari)

STATUS: Not Verified

The table gives CPU and GPU time ratios for four ranges.
The PR does not include the benchmark behind it.
The internal std_normal_erfcx helpers it compares against were removed in 5f33a10, so the comparison can no longer be rerun at the head in any case.

Evidence:

Table only.
The benchmark is missing.

Request: If this comparison still matters for the review, could you link the benchmark code that produced it?

---

FACT: The wiener4_lcdf_unnorm difference loses accuracy without bound

Quote:

The loss grows without bound.

(comment by avehtari)

STATUS: Not Verified

The table gives the ulp error of std_normal_lcdf(-x) - std_normal_lpdf(x) and of an erfcx form at x = 10, 30, 100 and 300.
The PR does not include the script that computed it, and the comment marks the change as out of scope for this PR.

Evidence:

Table only.
The script is missing.

Request: Could you link the script behind the table, or move this note to a separate issue for wiener4_lcdf_unnorm with the script attached?

---

FACT: Accuracy of std_normal_lcdf before and after 5f33a10

Quote:

Worst error in ulp (relative error / 2^-53) against mpmath at 50 digits, 2001 points per range, value / gradient:

(comment by avehtari)

STATUS: Not Verified

The table gives worst ulp errors over 2001 points per range against mpmath at 50 digits, before and after the commit, and the text adds 9.6 ulp on a denser grid.
The PR does not include the script.
The commit does add eight mpmath reference rows to std_normal_lcdf_impl_test.cpp, which the head passes at a relative tolerance of 1e-12, but those single points do not measure the ulp table.

Evidence:

Table only.
The mpmath script and the grids are missing.

Request: Could you add the script that produced the before and after table, including the 6.43 and 8.77 ulp central-range values?

---

FACT: 5f33a10 makes std_normal_lcdf_value_grad faster

Quote:

Time relative to before (below 1 is faster):

(comment by avehtari)

STATUS: Not Verified

The table gives time ratios of 0.66 to 1.02 for six ranges on a Xeon E5-2680 v3.
The PR does not include the benchmark that timed internal::std_normal_lcdf_value_grad.

Evidence:

Table and setup description only.
The benchmark source is missing.

Request: Could you add the benchmark source for this table?

---

FACT: Phi as opened (exp of log Phi) is 1.4 to 4.4 times slower than develop

Quote:

it takes 4.4 times as long for |x| < 0.66, 1.9 times on [-3, 3] and on [-37, -5], and 1.4 times on [5, 8].

(comment by avehtari)

STATUS: Not Verified

The comment reports Phi(double) and Phi(var) timings against develop, and a second table with the times after bf5861a.
The PR does not include the benchmark code.
Later commits replaced the exp(log Phi) form, so the head no longer has the code that was timed.

Evidence:

Timings in the text and a table.
The benchmark source is missing.

Request: Could you add the benchmark for Phi?
The current erfcx form could then be timed against develop with the same code.

---

FACT: Worst ulp error of Phi for develop, the original PR and bf5861a

Quote:

In the lower tail every version loses about x^2 / 2 ulp, because the argument x / sqrt(2) is rounded before erfc sees it, and the exp(log Phi) form has the same limit.

(comment by avehtari)

STATUS: Not Verified

The table gives worst ulp errors over 47001 points on [-38.5, 8.5].
The PR does not include the script.
The PR's Phi_test checks two points at 1e-14 relative error, which is covered in the fact check on develop's Phi defect.

Evidence:

Table only.
The mpmath script is missing.

Request: Could you add the script that computed this table?

---

FACT: Accuracy of the three cdfs for develop, the original PR and 0d4ac36

Quote:

Worst relative error in ulp against mpmath (100 digits; partials by mpmath.diff), develop / original PR / after 0d4ac36:

(comment by avehtari)

STATUS: Not Verified

The table gives worst ulp errors of the value and d/dy for z from -38 to 8 against mpmath at 100 digits.
The PR does not include the script.
9b38ecb later replaced the erfc form that the last column measures.

Evidence:

Table only.
The mpmath script is missing.

Request: Could you add the script, ideally run at the current head?

---

FACT: f4cd10e and c30ee6f leave all results bitwise identical

Quote:

They change only when Eigen evaluates an expression, so all results are bitwise identical.

(comment by avehtari)

STATUS: Not Verified

The two diffs only wrap existing expressions in eval(...), so the source arithmetic per element is the same.
Bitwise equality of the compiled results also depends on how the compiler vectorizes the lazy and the evaluated forms.
The PR does not include a test or script that compares the outputs before and after.
The same applies to the bitwise claim for 927099e.

Evidence:

A statement in the comment.
A comparison of outputs before and after is missing.

Request: Could you add the script that compared the outputs, or a sentence saying how bitwise equality was checked (which functions, inputs and compiler)?

---

FACT: The PR is now faster than develop for the functions in the description's table except exp_mod_normal_lcdf

Quote:

With the commits from my comments, the PR is now faster than develop for every function in the table except exp_mod_normal_lcdf, apart from exp_mod_normal_lpdf and skew_normal_lpdf with var, which are at develop's time.

(comment by avehtari)

STATUS: Not Verified

Two tables give the time relative to develop for ten functions at three versions on an AMD EPYC 9965, with double and var, next to the inverse of the description's speedups.
The PR does not include the benchmark.
These ratios are much closer to 1 than the description's.
The PR does not include code for either set of numbers.

Evidence:

Tables only.
The benchmark source and the inputs are missing.

Request: Could you add the benchmark used for these tables?
The description's table could then be replaced with numbers anyone can rerun.

---

FACT: The description's gains come from a slow C-library erf and erfc on its machine

Quote:

So the description measured a slow develop, not a faster branch.

(comment by avehtari)

STATUS: Not Verified

The comment compares absolute times from the description with times on an AMD EPYC 9965 and concludes that develop was slow on the description's machine.
This is an inference from two sets of timings, neither of which has code in the PR.
The description does not say which machine it used.

Evidence:

An inference from timing tables.
The benchmark and the description's platform are missing.

Request: Could the description name its OS, compiler and CPU, as the comment asks, and could both benchmarks be added so the comparison can be rerun?

---

FACT: The erfcx form of Phi (9b38ecb) is more accurate in the lower tail

Quote:

The new form is more accurate in the lower tail.

(comment by avehtari)

STATUS: Not Verified

The table adds a column for 9b38ecb to the earlier Phi table, with worst errors of 515 ulp and below in the tail.
The text also reports improvements of the cdf ratio.
The PR does not include the script.

Evidence:

Table only.
The mpmath script is missing.

Request: Could you add the script that produced this table?

---

FACT: 927099e brings skew_normal_lpdf with var back to 0.85 of develop's time

Quote:

That made double calls 7 % slower, because there z is read only twice and the allocation costs more than the second read.

(comment by avehtari)

STATUS: Not Verified

The comment gives a small table for skew_normal_lpdf on a Xeon Gold 6248 and notes that evaluating z in the double path made it 7 % slower.
The PR does not include the benchmark.

Evidence:

Table only.
The benchmark source is missing.

Request: Could you add the benchmark source for this table?

---

FACT: OpenCL device timings of Phi and std_normal_lcdf on a Tesla V100

Quote:

I also timed the OpenCL device functions on a Tesla V100, with each work item walking through the input range in registers and a baseline kernel subtracted.

(comment by avehtari)

STATUS: Not Verified

The comment reports that the device Phi from 9b38ecb is 2 to 4 times faster than the original form and up to 1.8 times slower than develop on [-5, 8].
The PR does not include the device benchmark kernel.

Evidence:

Numbers in the text only.
The device benchmark is missing.

Request: Could you add the device benchmark kernel and the host code that timed it?

---

FACT: After a17bba2 the worst exp_mod_normal error is 510 to 900 ulp with every result finite

Quote:

Worst error over the value and the four partials after the change, no non-finite result: 510 ulp (lcdf) and 480 ulp (lccdf) on the grid above, and 640 and 900 ulp on a denser grid of 2500 points with |z| in [1e-4, 100] and v in [1e-10, 100].

(comment by avehtari)

STATUS: Not Verified

The comment gives worst errors on two grids against mpmath.
The PR does not include the grid script.
The PR's exp_mod_normal_lcdf_test.cpp checks 13 reference rows and 6 deep-tail rows at a relative tolerance of 1e-12, which is looser than the 510 to 900 ulp figures, so it does not show them.

Evidence:

Numbers in the text only.
The grid and the mpmath script are missing.

Request: Could you add the grid script, or tighten the reference test tolerance to the level the comment reports?

---

FACT: a17bba2 makes the exp_mod_normal cdfs faster

Quote:

The change also makes the functions faster.

(comment by avehtari)

STATUS: Not Verified

The table gives time relative to develop before and after the commit on an AMD EPYC 9965.
The PR does not include the benchmark.

Evidence:

Table only.
The benchmark source is missing.

Request: Could you add the benchmark source for this table?

---

FACT: The OpenCL exp_mod_normal port agrees with the CPU code to 5e-14 relative

Quote:

Compiled on the host, the device code agrees with the CPU code to 5e-14 relative on the 2500-point grid (1.2e-12 at the dsigma zero)

(comment by avehtari)

STATUS: Not Verified

The comment compares the device code, compiled on the host, with the CPU code on a 2500-point grid.
The PR does not include that comparison.
The PR's OpenCL test exp_mod_normal_lcdf_reference_test.cpp compares the device results with mpmath references at a relative tolerance of 1e-12, which is a different check.

Evidence:

A number in the text only.
The host-compiled comparison is missing.

Request: Could you add the host comparison, or a test that compares the device and CPU results on the grid at the stated tolerance?

---

FACT: e1fea06 and cf177ef make several functions 1.2 to 1.4 times faster

Quote:

exp_mod_normal_cdf is 1.2 to 1.3 times faster after e1fea06, and after cf177ef std_normal_lcdf, normal_lcdf, lognormal_lcdf, exp_mod_normal_lpdf and skew_normal_lpdf are 1.2 to 1.4 times faster for z on [-30, 30].

(comment by avehtari)

STATUS: Not Verified

The comment reports speedups from two commits on an AMD EPYC 9965.
The PR does not include the benchmark binaries or their source.

Evidence:

Numbers in the text only.
The benchmark source is missing.

Request: Could you add the benchmark source used for these measurements?

---

FACT: erfcx and std_normal_lcdf stay within about 2.7 x 2^-52 of mpmath

Quote:

Accuracy against mpmath stays within about 2.7 x 2^-52 for erfcx and std_normal_lcdf.

(comment by avehtari)

STATUS: Not Verified

The PR does not include the mpmath comparison behind this bound.

Evidence:

A number in the text only.
The script is missing.

Request: Could you add the script that measured this bound?

---

FACT: 1f8b6ba cuts the non-finite results on a 105-point grid

Quote:

On a grid of 105 points from z = -1e300 to 10 and v = 1e-300 to 1e10, non-finite results drop from 38 to 0 for the cdf, from 1 to 0 for the lccdf, and from 38 to 28 for the lcdf

(comment by avehtari)

STATUS: Not Verified

The comment counts non-finite results on a grid of 105 points for the PR before and after the commit and for develop.
The PR does not include the grid script.
The six deep-tail reference rows that the commit adds are covered in a separate fact check.

Evidence:

Counts in the text only.
The grid script is missing.

Request: Could you add the grid script, or turn the grid into a test that counts non-finite values?

---

FACT: Restoring to_ref cuts the expression-test compile memory from 33 GB to 9 GB

Quote:

Peak g++ -O0 memory compiling the generated expression tests drops from 33 GB to 9 GB for each exp_mod_normal lcdf/lccdf/cdf and by about 35% for the normal and lognormal ones

(commit 201dcab9)

STATUS: Not Verified

The commit message reports peak g++ -O0 memory for the generated expression tests.
The PR does not include the measurement.

Evidence:

Numbers in the commit message only.
The memory measurement is missing.

Request: Could you add the command used to measure peak compile memory (for example /usr/bin/time -v on one generated expression test) and the output?

---

FACT: Reading y_val without to_ref failed the Prim expression tests

Quote:

This failed the Prim expression tests for std_normal_lcdf, std_normal_lccdf and std_normal_cdf.

(commit f37ea7cb)

STATUS: Not Verified

The commit message reports that the double expression argument was evaluated twice and that the Prim expression tests failed.
The PR does not include the failing output or a link to the CI run.

Evidence:

A statement in the commit message.
The failing test output is missing.

Request: Could you link the CI run that failed, or add the expression test output?

---

FACT: Skipping the far upper tail fixes the scalar exp_mod_normal_lpdf benchmark

Quote:

This fixes the scalar exp_mod_normal_lpdf benchmark (z ~ 9.4), where develop's erfc returns early.

(commit 20155578)

STATUS: Not Verified

The commit message says the early return for z > 8.5 fixes a scalar exp_mod_normal_lpdf benchmark at z of about 9.4.
The PR does not include that benchmark.
The code change itself is covered in the description check of 2015557.

Evidence:

A statement in the commit message.
The benchmark is missing.

Request: Could you add the benchmark, with the timings before and after?

---

FACT: Two var_value<Eigen::VectorXd> bugs in as_array_or_scalar and ref_type_if were found and fixed

Quote:

Two other bugs were found & fixed for var_value<Eigen::VectorXd> handling:

(PR description, Summary)

STATUS: Not Verified

The PR adds mathMixScalFun_normal_lccdf_varmat and mathMixScalFun_std_normal_lccdf_varmat, which compare a var_value<Eigen::VectorXd> argument with an Eigen::Matrix<var> argument.
Both pass at the head and on the merge base with develop.
They do not show a defect in develop's helpers.
The code changes themselves are covered in the description check above.
The bugs may only be reachable through the new code paths of this PR, but the PR does not say so or include a test that fails without the two changes.

Evidence:

Two varmat tests that pass on develop and at the head.
A test that fails without the two helper changes is missing.

Request: Could you add a test that fails without the as_array_or_scalar and ref_type_if changes, for example one that calls the helpers on an rvalue arena_matrix or on a var_value holding an array?
If the bugs only appeared with this PR's new code, could the description say so?

---

[stevieb-bot]

Fact check, part 2 of 2.

Description vs Code

Each section below lists facts from the PR summary and discussion that an agent was able to confirm, partially confirm, or not confirm at all, or that a run of the code in the PR contradicted.
Please review the "Partially Confirmed", "Not Verified", and "Contradicted" facts for facts missing evidence found by this agent.

Confirmed Facts

FACT: The value and gradient calculations move into reusable functions that other functions delegate to

Quote:

This PR extracts the key value and gradient calculations from the std_normal_lcdf header into reusable functions that other normal-cdf-related functions can delegate to, removing a lot of duplicated/inconsistent implementations.

(PR description, Summary)

STATUS: Confirmed

stan/math/prim/fun/std_normal_lcdf_impl.hpp defines internal::std_normal_lcdf_value_grad for scalars and arrays (L51-L161).
normal_lcdf, lognormal_lcdf, skew_normal_lpdf, exp_mod_normal_lpdf and std_normal_lcdf call it, and normal_lccdf, lognormal_lccdf and std_normal_lccdf reach it through the shared impl functions.
For example, normal_lcdf calls it at normal_lcdf.hpp#L54.

---

FACT: as_array_or_scalar returns an array unchanged and ref_type_if keeps an rvalue arena_matrix by value

Quote:

as_array_or_scalar() now returns an array unchanged instead of re-wrapping it [...]
ref_type_if keeps an rvalue arena_matrix by value rather than referencing a dead temporary

(PR description, Summary)

STATUS: Confirmed

In rev/fun/as_array_or_scalar.hpp the var_value overload returns the argument as is when its value type is already an Eigen array, and calls .array() otherwise (L19-L26).
In ref_type.hpp the arena_matrix specialization passes T_base&& for an rvalue.
The Eigen specialization (L32-L47) maps that to T_base held by value.
For an lvalue it passes T_base, which the Eigen specialization maps to const T_base& (L49-L56).
On develop both cases became const T_base&.

---

FACT: New stability tests and new mix tests for the distributions are added

Quote:

Additional tests exercising the increased numerical stability are added, as well new mix tests for the distributions

(PR description, Tests)

STATUS: Confirmed

The PR adds new mix test files for exp_mod_normal (lpdf, lcdf, lccdf, cdf), lognormal (lcdf, lccdf, cdf), skew_normal and the shared kernel, and extends the normal and std_normal mix tests.
It also adds tail and reference tests in prim and rev (for example std_normal_lcdf_impl_test.cpp, exp_mod_normal_lcdf_test.cpp and normal_lcdf_small_sigma_test.cpp) and OpenCL tail tests.

---

FACT: 5f33a10 makes the CPU and OpenCL kernels reuse the erfcx helpers and adds central-range reference rows

Quote:

I deleted std_normal_erfcx_tail and std_normal_erfcx_middle; the kernel calls erfcx_cody_tail and erfcx_cody_middle instead.
[...]
OpenCL: the std_normal_lcdf device function calls erfcx_tail_correction, erfcx_cody_tail, erfcx_cody_middle and erfcx_small from the erfcx device function, and I deleted its own coefficient tables.

(comment by avehtari)

STATUS: Confirmed

At 5f33a10 the CPU kernel calls erfcx_cody_tail, erfcx_cody_middle and erfcx_small, and std_normal_erfcx_tail and std_normal_erfcx_middle are gone.
At the reviewed commit the kernel calls the _fraction forms that cf177ef added (std_normal_lcdf_impl.hpp#L111-L113), and std_normal_tail_correction calls erfcx_tail_polynomials with reverse iterators (L33-L39).
The central range computes LOG_HALF - (0.5 * z) * z + log(erfcx_small(-z / sqrt(2))) and the slope SQRT_TWO_OVER_SQRT_PI / erfcx_small(...) (L98-L108).
The OpenCL device function calls erfcx_tail_correction, erfcx_cody_tail, erfcx_cody_middle and erfcx_small and does not have coefficient tables of its own (device_functions/std_normal_lcdf.hpp#L19-L45).
elt_function_cl.hpp registers std_normal_lcdf_impl and std_normal_lcdf_derivative with erfcx_device_function first.
prim/prob/std_normal_lcdf_impl_test.cpp has rows at -3, -0.6, -0.25, 0.25, 0.6, 0.66, 0.67 and 3.
The PR does not include the mpmath script that produced those reference values.

---

FACT: bf5861a computes Phi as 0.5 * erfc(-x / sqrt(2)) and restores the develop OpenCL registrations

Quote:

I pushed commit bf5861a, which computes Phi that way on the CPU and in the OpenCL device function.

(comment by avehtari)

STATUS: Confirmed

In bf5861a, Phi(double) returns 0.5 * erfc(-INV_SQRT_TWO * x), and elt_function_cl.hpp registers Phi with only phi_device_function and inv_Phi without the std_normal_lcdf device function.
9b38ecb later replaced this form with the erfcx form, which is what the reviewed commit has (Phi.hpp#L61-L72).

---

FACT: 0d4ac36 computes the three cdfs from erfc through a new std_normal_cdf_value_ratio helper

Quote:

A new helper, internal::std_normal_cdf_value_ratio in std_normal_lcdf_impl.hpp, returns Phi(z) and phi(z) / Phi(z) for a scalar or an Eigen array, generic in the scalar type so that nested autodiff works.

(comment by avehtari)

STATUS: Confirmed

In 0d4ac36, internal::std_normal_cdf_value_ratio returns 0.5 * erfc(-INV_SQRT_TWO * z) and the ratio, with scalar and array overloads, and std_normal_cdf, normal_cdf and lognormal_cdf call it.
At the reviewed commit the helper uses the erfcx form from 9b38ecb instead (std_normal_lcdf_impl.hpp#L163-L282).
It still returns a ratio of 0 where erfcx(-u) is 0, at z = -infinity.

---

FACT: The OpenCL cdfs still use exp(lcdf)

Quote:

I did not change the OpenCL versions, which still use exp(lcdf).

(comment by avehtari)

STATUS: Confirmed

At the reviewed commit the OpenCL std_normal_cdf, normal_cdf and lognormal_cdf return exp(...) of the matching lcdf impl (opencl/prim/normal_cdf.hpp#L32, std_normal_cdf.hpp#L24, lognormal_cdf.hpp#L19).

---

FACT: f4cd10e and c30ee6f evaluate the reused temporaries once

Quote:

f4cd10e wraps both in eval(...).

(comment by avehtari)

STATUS: Confirmed

f4cd10e wraps z in eval(...) in lognormal_cdf and lognormal_lcdf_impl.
Commit c30ee6f wraps z and scaled_slope in normal_lcdf_impl and scaled_slope in lognormal_lcdf_impl, and diff, z, v, log_b, lp, s_b and q in exp_mod_normal_lcdf_impl.
At the reviewed commit normal_lcdf still evaluates z and scaled_slope (normal_lcdf.hpp#L52-L56).
Commit a17bba2 later replaced the exp_mod_normal code that c30ee6f changed.

---

FACT: 9b38ecb computes Phi and the cdf helper from the Cody erf rational and erfcx, with only exp from the C library

Quote:

Only exp comes from the C library.

(comment by avehtari)

STATUS: Confirmed

Phi(double) returns 0.5 + 0.5 * erf_cody_small(u) for |u| < 0.46875, and otherwise 0.5 * exp(-(0.5 * x) * x) times erfcx_cody_tail or erfcx_cody_middle of |u|, or one minus that for x >= 0 (Phi.hpp#L61-L72).
The scalar std_normal_cdf_value_ratio computes the same value in the same order and returns the ratio SQRT_TWO_OVER_SQRT_PI / c for z < 0 and SQRT_TWO_OVER_SQRT_PI / erfcx_small(-u) in the central range (std_normal_lcdf_impl.hpp#L173-L205).
The array overload for double, added later in 2015557, uses a vectorized exp and documents that its values can differ from Phi(double) in the last bit.
The OpenCL Phi device function uses the same formula, and elt_function_cl.hpp registers Phi and inv_Phi with erfcx_device_function.

---

FACT: 927099e evaluates z and score once in the autodiff path of skew_normal_lpdf and keeps z lazy for double

Quote:

New commit 927099e evaluates both once in the autodiff branch.

(comment by avehtari)

STATUS: Confirmed

In the autodiff branch skew_normal_lpdf evaluates z_eval = eval(z) and score = eval(...) once and uses them in the partials and the sum.
The value-only branch uses the lazy z (skew_normal_lpdf.hpp#L68-L110).

---

FACT: a17bba2 computes the exp_mod_normal cdfs from erfcx differences with a shared scalar loop

Quote:

The difference E(u) - E(t) is computed directly when E(t) < 0.99 E(u), and otherwise from a series in the derivatives of erfcx at the midpoint, or for arguments of at least 4 from the difference of the Cody tail rational; the rounded t - u is never used.

(comment by avehtari)

STATUS: Confirmed

exp_mod_normal_erfcx_diff takes the direct difference when e_hi < 0.99 * e_lo, and otherwise exp_mod_normal_erfcx_diff_close uses the Cody tail form for lo >= 4 and the midpoint series in the odd derivatives below that (exp_mod_normal_lcdf.hpp#L74-L196).
For z > 0 the kernel uses F = 1 - p where F >= 1/4, and where F < 1/4 it uses g (E(u) - E(t)) for t >= 0 or -expm1(log_b) + g (E(-t) - E(-u)) for t < 0 (L330-L427).
The lcdf, lccdf and cdf all call exp_mod_normal_lcdf_impl, which loops over the elements once (L434-L492).

---

FACT: 39886e0 adds the exp_mod_normal_lcdf_part device function and eight binary wrappers

Quote:

The kernel generator only has functions with one output, so eight binary wrappers (exp_mod_normal_{lcdf,lccdf}_{value,d_z,d_v,d_s}) select the value or one of the three partials, and opencl/prim/exp_mod_normal_lcdf.hpp uses them.

(comment by avehtari)

STATUS: Confirmed

The device function exp_mod_normal_lcdf_part(double z, double v, int upper, int which) is at device_functions/exp_mod_normal_lcdf.hpp#L284-L285, followed by the eight wrappers exp_mod_normal_{lcdf,lccdf}_{value,d_z,d_v,d_s} (L374 to L393).
elt_function_cl.hpp registers all eight as binary functions, and opencl/prim/exp_mod_normal_lcdf.hpp calls them (L66-L89).

---

FACT: d081fff keeps the multiplication by 1 / sigma and divides only where it overflows

Quote:

I pushed d081fff, which keeps your multiplication and recomputes z with a division only when some inv(sigma) is infinite.

(comment by avehtari)

STATUS: Confirmed

In d081fff, normal_lcdf and lognormal_lcdf compute z with inv_sigma and recompute it with a division when any(inv_sigma == INFTY), and the commit adds normal_lcdf_small_sigma_test.cpp.
Commit e076f6a later replaced this with a division in all cases.
The reviewed commit always divides z by sigma (normal_lcdf.hpp#L52).

---

FACT: 1f8b6ba adds six mpmath rows for the deep tail to the CPU and OpenCL reference tests

Quote:

The commit adds six mpmath rows for the deep tail to the CPU and OpenCL reference tests.

(comment by avehtari)

STATUS: Confirmed

exp_mod_normal_lcdf_reference.hpp defines exp_mod_normal_lcdf_deep_tail_reference with six rows (L83-L104).
The CPU test ProbDistributionsExpModNormal_deep_tail_reference and the OpenCL test opencl_matches_deep_tail_reference both use it.
The header comment says the rows come from mpmath at 800 digits.
The script is not in the PR.

---

FACT: e1fea06 adds a log_value flag and a one-pass argument check

Quote:

exp_mod_normal_check_args checks y, mu, sigma and lambda in one branch-free pass and reports whether y holds an infinity

(commit e1fea065)

STATUS: Confirmed

exp_mod_normal_lcdf_kernel and exp_mod_normal_lcdf_impl take a log_value template flag, and exp_mod_normal_cdf calls the impl with it unset (exp_mod_normal_cdf.hpp#L15).
exp_mod_normal_check_args checks all four arguments in one loop, returns whether y has an infinity, and reruns the elementwise checks only when that loop fails (exp_mod_normal_check_args.hpp#L36-L69).

---

FACT: 2015557 returns early past z = 8.5 and vectorizes exp in the cdf ratio

Quote:

std_normal_lcdf_value_grad returns -tail and the density slope for z > 8.5 inside its existing z > 40 check

(commit 20155578)

STATUS: Confirmed

The scalar kernel returns -tail and the density slope for z > 8.5, with the z > 40 return nested inside that check (std_normal_lcdf_impl.hpp#L55-L74).
The double array overload of std_normal_cdf_value_ratio computes gs = 0.5 * exp(-0.5 * values.square()) once before its loop unless every |u| < 0.46875, and uses the scalar exp for |u| > 26.5 (L214-L267).

---

FACT: e076f6a divides z by sigma and multiplies only the slopes by 1 / sigma

Quote:

z now divides by sigma as on develop; only the slopes multiply by 1 / sigma, where an infinite slope is the right limit.

(commit e076f6a8)

STATUS: Confirmed

normal_lcdf_impl computes z = eval(sign * (y_val - mu_val) / sigma_val) and scaled_slope = eval(slopes * inv(sigma_val)) (normal_lcdf.hpp#L52-L56).
lognormal_lcdf_impl does the same (lognormal_lcdf.hpp#L56-L59).

---

FACT: 8cff6fd fuses the lpdf reductions and adds vectorized edge-case tests

Quote:

skew_normal_lpdf and exp_mod_normal_lpdf now reduce the log density with a single sum over one elementwise expression instead of separate sums scaled by N / size(arg).

(commit 8cff6fd3)

STATUS: Confirmed

skew_normal_lpdf sums one elementwise expression plus log_normalizer in each propto case (skew_normal_lpdf.hpp#L94-L101), and exp_mod_normal_lpdf does the same (exp_mod_normal_lpdf.hpp#L97-L110).
The tests the message lists are in the PR, namely expect_propto_vectorized.hpp with the rev lpdf tests, the prim closed-form, error and size tests in exp_mod_normal_lpdf_test.cpp and skew_normal_lpdf_test.cpp, and the mix vector-tail tests.
All of them pass at the head in cpu_suites-head.out.

---

Partially Confirmed Facts

FACT: The Taylor expansions are replaced with the Cody rational approximations

Quote:

The value and gradient calculations are also updated to replace the taylor expansions with the Cody rational approximations (previously only used for the tail).

(PR description, Summary)

STATUS: Partially Confirmed

The Taylor expansions of develop's std_normal_lcdf are gone at the reviewed commit.
For |z| / sqrt(2) of 0.46875 and above the kernel uses the Cody (1969) rationals from erfcx.hpp (erfcx_cody_middle_fraction and erfcx_cody_tail_fraction, std_normal_lcdf_impl.hpp#L111-L113).
In the central range it uses erfcx_small (L98-L108), which erfcx.hpp documents as a degree-18 Chebyshev-economized polynomial.
Commit 5f33a10 replaced the Cody erf rational there.

Request: Could the description say that the central range uses the erfcx_small polynomial and the other ranges use the Cody rationals?

---

FACT: normal, lognormal and exp_mod_normal lcdf and lccdf go through the shared std_normal_lcdf implementation, including OpenCL

Quote:

Implements most of #3364: the normal, lognormal and exp_mod_normal lcdf/lccdf through the shared std_normal_lcdf implementation, and the OpenCL part.

(PR description)

STATUS: Partially Confirmed

On the CPU, normal_lcdf and lognormal_lcdf call internal::std_normal_lcdf_value_grad, and normal_lccdf and lognormal_lccdf call the same impl with reflect set.
exp_mod_normal_lcdf, _lccdf and _cdf do not.
Since a17bba2 they compute the cdf from differences of erfcx in exp_mod_normal_lcdf.hpp, which includes erfcx.hpp and not std_normal_lcdf_impl.hpp.
On OpenCL, normal_lcdf and lognormal_lcdf use the std_normal_lcdf_impl device function, while the exp_mod_normal lcdf, lccdf and cdf use the new exp_mod_normal_lcdf_part device function (opencl/prim/exp_mod_normal_lcdf.hpp#L66-L89).

Request: Could the description say that exp_mod_normal lcdf and lccdf share the erfcx helpers and have their own kernel, rather than the std_normal_lcdf implementation?

---

FACT: The deep lower tail keeps log1p(tail_term) and neg_z / (1 + tail_term), and std_normal_extreme_tail_derivatives checks down to z = -1e50

Quote:

The deep lower tail (z <= -4 sqrt(2)) keeps its log1p(tail_term) value and its neg_z / (1 + tail_term) gradient; only tail_term now comes from the shared code.
[...]
std_normal_extreme_tail_derivatives checks third derivatives down to z = -1e50.

(comment by avehtari)

STATUS: Partially Confirmed

The lower-tail branch with log1p(tail_term) and neg_z / (1.0 + tail_term) is still in the kernel (std_normal_lcdf_impl.hpp#L75-L88).
Since cf177ef it runs only for autodiff types, and double inputs take the z < 0 branch.
The test named std_normal_extreme_tail_derivatives checks rows at -1e4, -1e8 and -1e10 (mix/prob/std_normal_lcdf_impl_test.cpp#L9-L36).
The -1e50 row is in the next test, std_normal_extreme_tail_vector_derivatives (L38-L55).
Both pass at the head in cpu_suites-head.out.

Request: Could the comment name std_normal_extreme_tail_vector_derivatives as the test that reaches -1e50, and say that the log1p branch now applies only to autodiff types?

---

FACT: cf177ef makes the array path one pass and shortens the division chains

Quote:

The scalar function no longer branches on the sign of z for doubles.

(commit cf177ef9)

STATUS: Partially Confirmed

The array std_normal_lcdf_value_grad loops once with coeff and coeffRef (std_normal_lcdf_impl.hpp#L139-L161).
The lower-tail log1p branch runs only for autodiff types (L75-L88), the tail correction multiplies by SQRT_PI (L38), erfcx_cody_tail_fraction puts the tail rational over one denominator, and the z < 0 slope divides the denominator by the numerator (L118).
For double the kernel forms |z| with std::fabs (L93-L97), but it still branches on z < 0 at L114 to choose between the log and log1p forms.
The sentence about the sign branch therefore holds only for the step that forms |z|.
The part about normal_lcdf multiplying z by 1 / sigma was later replaced in e076f6a.

Request: Could the commit message (or a follow-up note) say that the removed sign branch is the one that formed |z|, since the kernel still branches on z < 0 for the value?

---

## Not Checked

• Commit headlines that only name a change, such as 'Remove c++20 code', 'Use TEST_F instead of raw TEST' and 'fix neg_binomial_2_lpmf.hpp for log1p', are not listed as separate fact checks.
• Smaller timing remarks in the comments are not listed separately, for example the 1 to 6 % cost of the range check in 1f8b6ba, the 0.99 to 1.01 timing of d081fff, the 3 % exp_mod_normal_lpdf result and the GCC inlining note.
  As with the timing tables above, the PR does not include the benchmark code behind them.
• Statements about an earlier commit are checked against the merge base with develop and the reviewed commit, and description statements also by reading that commit.
  The intermediate commits were not built.
• The generated expression tests were not run, because they need the stanc signatures, which were not fetched.
• The stan-buildbot end-to-end benchmark table, which no author cites as evidence for a claim.
• The clang-format and cpplint results in avehtari's comments, which are repository checks and not evidence for a claim about the code.