Add erfcx function - #3405
Add erfcx function#3405
Conversation
|
So I'm fine with this version. My one question is whether we want to backport something like https://jugit.fz-juelich.de/mlz/lib/cerf/-/blob/main/lib/erfcx.c?ref_type=heads |
|
I'll make a comparison |
|
@WardBrian just pointed out that the file I pointed to pulls in this file of a huge array of precomputed chebyshev polynomials. Not sure if we want that. I'll explore this a bit, but my default is that the PR as of now seems reasonable and good. |
|
Fetched and benchmarked it against the same
|
|
It seems it would be possible to get maybe 20% speedup for this PR, but to get that 2.3 faster on [0, 4] would require adopting that 464-line generated coefficient table |
|
Thinking about it I think I'm okay with backporting the chebyshev polynomial code as well. If we are going to use this in a lot of places then I think it makes sense to have that lookup table. |
|
I'll create another branch using libcerf code and make more comparisons |
|
I did some studying (with help from Claude). One additional consideration is that that huge array of precomputed chebyshev polynomials is GPU-hostile GPU results — Tesla V100-SXM2-16GB, OpenCL 3.0 CUDA, 80 CUs,
Where cody is the algorithm used in this PR. Do we want to use different algorithm for CPU (libcerf with Chebyshev) and GPU (cody)? |
Jenkins Console Log Machine informationDistributor 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.269 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 |
Oh yeah if we do not leave the table over there that is annoying to transfer. It would be nice to have the more accurate CPU version, but I'm also fine with the cody algorithm this PR uses. Let me have one more look at this over the morning. |
|
I think we can have a Cody version which is
I'll update the PR |
Yes and just better than what we are currently so idt we need to chase something perfect |
|
I just updated the PR and the PR text which explains the current PR algorithm version which beats libcerf for our purposes |
|
This PR is still missing some OpenCL support and OpenCL unit tests. Working on it, and pushing commits tomorrow |
Jenkins Console Log Machine informationDistributor 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.978 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 |
andrjohns
left a comment
There was a problem hiding this comment.
A couple suggestions for reducing duplication and squeezing out a bit more performance
| inline double erfcx_cody_tail(double x) { | ||
| const double u = 1.0 / (x * x); | ||
| double p = 0.0163153871373020978498; | ||
| p = 0.305326634961232344035 + u * p; | ||
| p = 0.360344899949804439429 + u * p; | ||
| p = 0.125781726111229246204 + u * p; | ||
| p = 0.0160837851487422766278 + u * p; | ||
| p = 0.000658749161529837803157 + u * p; | ||
| double q = -1.0; | ||
| q = -2.56852019228982242072 + u * q; | ||
| q = -1.87295284992346047209 + u * q; | ||
| q = -0.527905102951428412248 + u * q; | ||
| q = -0.0605183413124413191178 + u * q; | ||
| q = -0.00233520497626869185443 + u * q; | ||
| return (INV_SQRT_PI + (p / q) * u) / x; | ||
| } |
There was a problem hiding this comment.
If you split the tail correction from the calculation then it can be reused elsewhere (both in this PR and others):
I also find the arrays+looping expression for the coefficients more readable (and compiles to the same code), but that's not a hard requirement!
| inline double erfcx_cody_tail(double x) { | |
| const double u = 1.0 / (x * x); | |
| double p = 0.0163153871373020978498; | |
| p = 0.305326634961232344035 + u * p; | |
| p = 0.360344899949804439429 + u * p; | |
| p = 0.125781726111229246204 + u * p; | |
| p = 0.0160837851487422766278 + u * p; | |
| p = 0.000658749161529837803157 + u * p; | |
| double q = -1.0; | |
| q = -2.56852019228982242072 + u * q; | |
| q = -1.87295284992346047209 + u * q; | |
| q = -0.527905102951428412248 + u * q; | |
| q = -0.0605183413124413191178 + u * q; | |
| q = -0.00233520497626869185443 + u * q; | |
| return (INV_SQRT_PI + (p / q) * u) / x; | |
| } | |
| /** | |
| * Correction factor of the Cody (1969) third-interval rational: | |
| * `erfcx(x) = (INV_SQRT_PI + u * correction(u)) / x` with `u = 1 / x^2`. | |
| * | |
| * @tparam T scalar type | |
| * @param u inverse square of the argument, `0 <= u <= 1/16` | |
| * @return `P(u) / Q(u)` | |
| */ | |
| template <typename T> | |
| inline T erfcx_tail_correction(const T& u) { | |
| static constexpr double p[] | |
| = {0.000658749161529837803157, 0.0160837851487422766278, | |
| 0.125781726111229246204, 0.360344899949804439429, | |
| 0.305326634961232344035, 0.0163153871373020978498}; | |
| static constexpr double q[] | |
| = {-0.00233520497626869185443, -0.0605183413124413191178, | |
| -0.527905102951428412248, -1.87295284992346047209, | |
| -2.56852019228982242072, -1.0}; | |
| T numerator = p[5]; | |
| T denominator = q[5]; | |
| for (int i = 4; i >= 0; --i) { | |
| numerator = p[i] + u * numerator; | |
| denominator = q[i] + u * denominator; | |
| } | |
| return numerator / denominator; | |
| } | |
| inline double erfcx_cody_tail(double x) { | |
| const double u = 1.0 / (x * x); | |
| return (INV_SQRT_PI + u * erfcx_tail_correction(u)) / x; | |
| } |
| q = -0.00233520497626869185443 + u * q; | ||
| return (INV_SQRT_PI + (p / q) * u) / x; | ||
| } | ||
|
|
There was a problem hiding this comment.
With the tail correction split out you can create a reusable helper for the derivative:
/**
* Derivative of `erfcx`, `2 * x * erfcx(x) - 2 / sqrt(pi)`. That difference
* cancels for large `x`; there the tail rational gives it directly.
*
* @tparam T scalar type
* @param x argument
* @param value `erfcx(x)`
* @return derivative of `erfcx` at `x`
*/
template <typename T>
inline T erfcx_derivative(const T& x, const T& value) {
if (x >= 4.0) {
const T u = 1.0 / (x * x);
return 2.0 * u * erfcx_tail_correction(u);
}
return 2.0 * x * value - TWO_OVER_SQRT_PI;
}
| inline double erfcx_small(double x) { | ||
| double p = 3.05977060678449757e-06; | ||
| p = -9.35890030086883823e-06 + x * p; | ||
| p = 2.46655529768908249e-05 + x * p; | ||
| p = -7.08163358203131886e-05 + x * p; | ||
| p = 1.98445679338826757e-04 + x * p; | ||
| p = -5.34506929034156810e-04 + x * p; | ||
| p = 1.38888415444527033e-03 + x * p; | ||
| p = -3.47359067853470795e-03 + x * p; | ||
| p = 8.33333374332981443e-03 + x * p; | ||
| p = -1.91048337772546720e-02 + x * p; | ||
| p = 4.16666666458337179e-02 + x * p; | ||
| p = -8.59717459974174147e-02 + x * p; | ||
| p = 1.66666666667239644e-01 + x * p; | ||
| p = -3.00901111227312890e-01 + x * p; | ||
| p = 4.99999999999992839e-01 + x * p; | ||
| p = -7.52252778063651983e-01 + x * p; | ||
| p = 1.0 + x * p; | ||
| p = -1.12837916709551256 + x * p; | ||
| return 1.0 + x * p; | ||
| } |
There was a problem hiding this comment.
By splitting the polynomial into the odd & even powers you can split the dependency chain and enable a bit more vectorisation in the compiler:
| inline double erfcx_small(double x) { | |
| double p = 3.05977060678449757e-06; | |
| p = -9.35890030086883823e-06 + x * p; | |
| p = 2.46655529768908249e-05 + x * p; | |
| p = -7.08163358203131886e-05 + x * p; | |
| p = 1.98445679338826757e-04 + x * p; | |
| p = -5.34506929034156810e-04 + x * p; | |
| p = 1.38888415444527033e-03 + x * p; | |
| p = -3.47359067853470795e-03 + x * p; | |
| p = 8.33333374332981443e-03 + x * p; | |
| p = -1.91048337772546720e-02 + x * p; | |
| p = 4.16666666458337179e-02 + x * p; | |
| p = -8.59717459974174147e-02 + x * p; | |
| p = 1.66666666667239644e-01 + x * p; | |
| p = -3.00901111227312890e-01 + x * p; | |
| p = 4.99999999999992839e-01 + x * p; | |
| p = -7.52252778063651983e-01 + x * p; | |
| p = 1.0 + x * p; | |
| p = -1.12837916709551256 + x * p; | |
| return 1.0 + x * p; | |
| } | |
| inline double erfcx_small(double x) { | |
| static constexpr double even[] | |
| = {1.0, | |
| 1.0, | |
| 4.99999999999992839e-01, | |
| 1.66666666667239644e-01, | |
| 4.16666666458337179e-02, | |
| 8.33333374332981443e-03, | |
| 1.38888415444527033e-03, | |
| 1.98445679338826757e-04, | |
| 2.46655529768908249e-05, | |
| 3.05977060678449757e-06}; | |
| static constexpr double odd[] | |
| = {-1.12837916709551256, -7.52252778063651983e-01, | |
| -3.00901111227312890e-01, -8.59717459974174147e-02, | |
| -1.91048337772546720e-02, -3.47359067853470795e-03, | |
| -5.34506929034156810e-04, -7.08163358203131886e-05, | |
| -9.35890030086883823e-06}; | |
| const double x2 = x * x; | |
| double e = even[9]; | |
| double o = odd[8]; | |
| for (int i = 8; i >= 1; --i) { | |
| e = even[i] + x2 * e; | |
| o = odd[i - 1] + x2 * o; | |
| } | |
| return even[0] + x2 * e + x * o; | |
| } |
| inline var erfcx(const var& a) { | ||
| double val = erfcx(a.val()); | ||
| return make_callback_var(val, [a, val](auto& vi) mutable { | ||
| a.adj() += vi.adj() * (2.0 * a.val() * val - TWO_OVER_SQRT_PI); | ||
| }); | ||
| } | ||
|
|
||
| /** | ||
| * The scaled complementary error function for matrix variables. | ||
| * | ||
| * @tparam T a matrix type | ||
| * @param a The variable. | ||
| * @return Scaled complementary error function applied elementwise. | ||
| */ | ||
| template <typename T, require_matrix_t<T>* = nullptr> | ||
| inline auto erfcx(const var_value<T>& a) { | ||
| auto val = to_arena(erfcx(a.val())); | ||
| return make_callback_var(val, [a, val](auto& vi) mutable { | ||
| a.adj().array() | ||
| += vi.adj().array() | ||
| * (2.0 * a.val().array() * val.array() - TWO_OVER_SQRT_PI); | ||
| }); | ||
| } |
There was a problem hiding this comment.
You can combine these and then call the pre-defined helper:
template <
typename T, require_var_t<T>* = nullptr,
require_all_not_nonscalar_prim_or_rev_kernel_expression_t<T>* = nullptr>
inline auto erfcx(T&& a) {
auto val = to_arena(erfcx(a.val()));
return make_callback_var(val, [a](auto& vi) mutable {
const auto& deriv = apply_scalar_binary(
[](double x, double v) { return internal::erfcx_derivative(x, v); },
a.val(), val);
as_array_or_scalar(a.adj())
+= as_array_or_scalar(vi.adj()) * as_array_or_scalar(deriv);
});
}| template <typename T> | ||
| inline fvar<T> erfcx(const fvar<T>& x) { | ||
| T v = erfcx(x.val_); | ||
| return fvar<T>(v, x.d_ * (2.0 * x.val_ * v - TWO_OVER_SQRT_PI)); |
There was a problem hiding this comment.
After introducing a reusable erfcx_derivative helper:
| return fvar<T>(v, x.d_ * (2.0 * x.val_ * v - TWO_OVER_SQRT_PI)); | |
| return fvar<T>(v, x.d_ * internal::erfcx_derivative(x.val_, v)); |
|
Thanks, @andrjohns, I included these. The second one found a real defect. The derivative suggestion revealed a defect
The cancellation is exact and removable: with no subtraction left. I missed this because of accidentally not measuring accuracy of Splitting out the tail correctionAdopted, in the array form. I checked it is bit-identical over 2 million points in Splitting the small-branch polynomial into even and odd powersAdopted. It is a real speed win and a small accuracy cost, so here are both numbers. Over
1.80 times faster by latency, 1.67 by throughput, for 0.7 ulp. The whole function is 7.0 ulp, set by the middle Cody interval, so the small branch was never the limit and the trade is clearly worth taking. Effect on the whole function. CPU, Xeon E5-2680 v3, one dedicated core, throughput, ns per call:
The same derivative defect was in the OpenCL path, and is now fixedYour suggestions are all CPU-side, but the second one applies to OpenCL as well. A.adj() += elt_multiply(res.adj(),
elt_multiply(2.0, elt_multiply(A.val(), res.val())) - TWO_OVER_SQRT_PI);so the OpenCL reverse mode had the same 2.55e+11 ulp error. The OpenCL device function now has The existing OpenCL suite never reached the tail at all. Its inputs were Speed on the GPUThe three CPU-side suggestions do not change the OpenCL device function, so the GPU timings are the same before and after. Tesla V100-SXM2-32GB, ns per call:
I then measured whether the even and odd split is worth porting to the OpenCL device function as well. It is not:
Nothing outside noise, against 1.67 to 1.80 on the CPU. The split shortens the dependency chain but does not reduce the instruction count; it raises it slightly, to two chains of 9 plus a combine. A GPU hides chain latency across warps and is limited by the instruction count, so it gains nothing. The OpenCL device function therefore keeps the single chain and stays at 1.5 ulp on that branch while the CPU is at 2.2 ulp. That 0.7 ulp difference is far inside the Merging the reverse-mode overloadsAdopted. Two notes on the snippet. The lambda captures Test statusBoth configurations pass. Plain CPU, no With QuestionTo call the derivative from the kernel generator I registered |
|
Thank you @avehtari for putting this together! |
SteveBronder
left a comment
There was a problem hiding this comment.
One note, generally anywhere we have static constexpr double {NAME}[] I would prefer we use an std::array instead of double {NAME}[] for a raw array
Overall looks good though!
Co-authored-by: Steve Bronder <Stevo15025@gmail.com>
|
@SteveBronder I committed your suggestions, but they are not actually valid. Fixing... |
Two bugs in the applied batch
static_assert(even.size() == odd.size(), "the two chains are stepped by one loop");
for (int i = static_cast<int>(even.size()) - 1; i >= 0; --i) {That is bit-identical to the previous version at all 200 001 points in the branch interval, and clean under ASan. The OpenCL device function had C++ in it.
Speed, measuredXeon E5-2680 v3, one dedicated core, minimum of three repetitions. Tesla V100-SXM2-32GB, 1048576 work items, best of five, baseline kernel subtracted. The CPU changes are neutral. Every range comes out at 1.00, so I kept them: The OpenCL rewrite costs time, so I reverted that file:
ns per call. Rewriting it with the coefficients written directly into the My guess is that the pairing trades a dependency chain for extra live values, which helps a single scalar CPU stream but not a GPU that already hides chain latency across warps. The same thing happened with the even and odd split on the CPU side: 1.67 to 1.80 times faster on the Xeon, and nothing on the V100. |
| const auto& deriv = apply_scalar_binary( | ||
| [](double x, double v) { return internal::erfcx_derivative(x, v); }, |
There was a problem hiding this comment.
| const auto& deriv = apply_scalar_binary( | |
| [](double x, double v) { return internal::erfcx_derivative(x, v); }, | |
| auto deriv = apply_scalar_binary( | |
| [](double x, double v) { return internal::erfcx_derivative(x, v); }, |
| struct { | ||
| double x; | ||
| double d; | ||
| } cases[] = {{-2.0, -436.89199672700738}, {-1.0, -11.14633932862008}, | ||
| {0.0, -1.1283791670955126}, {1.0, -0.27321201478389856}, | ||
| {5.0, -0.021332789764826311}, {10.0, -0.0055593122190608565}}; |
There was a problem hiding this comment.
Replace this with an std::pair
| inline double erfcx_cody_middle(double y) { | ||
| double p = 2.15311535474403846e-8 * y; | ||
| p = (p + 5.64188496988670089e-1) * y; | ||
| p = (p + 8.88314979438837594) * y; | ||
| p = (p + 66.1191906371416295) * y; | ||
| p = (p + 298.635138197400131) * y; | ||
| p = (p + 881.952221241769090) * y; | ||
| p = (p + 1712.04761263407058) * y; | ||
| p = (p + 2051.07837782607147) * y; | ||
| double q = y; | ||
| q = (q + 15.7449261107098347) * y; | ||
| q = (q + 117.693950891312499) * y; | ||
| q = (q + 537.181101862009858) * y; | ||
| q = (q + 1621.38957456669019) * y; | ||
| q = (q + 3290.79923573345963) * y; | ||
| q = (q + 4362.61909014324716) * y; | ||
| q = (q + 3439.36767414372164) * y; | ||
| return (p + 1230.33935479799725) / (q + 1230.33935480374942); | ||
| } |
There was a problem hiding this comment.
Sorry I didn't see in my last review comment I was in the opencl and those were meant for the cpu verions
| inline double erfcx_cody_middle(double y) { | |
| double p = 2.15311535474403846e-8 * y; | |
| p = (p + 5.64188496988670089e-1) * y; | |
| p = (p + 8.88314979438837594) * y; | |
| p = (p + 66.1191906371416295) * y; | |
| p = (p + 298.635138197400131) * y; | |
| p = (p + 881.952221241769090) * y; | |
| p = (p + 1712.04761263407058) * y; | |
| p = (p + 2051.07837782607147) * y; | |
| double q = y; | |
| q = (q + 15.7449261107098347) * y; | |
| q = (q + 117.693950891312499) * y; | |
| q = (q + 537.181101862009858) * y; | |
| q = (q + 1621.38957456669019) * y; | |
| q = (q + 3290.79923573345963) * y; | |
| q = (q + 4362.61909014324716) * y; | |
| q = (q + 3439.36767414372164) * y; | |
| return (p + 1230.33935479799725) / (q + 1230.33935480374942); | |
| } | |
| inline double erfcx_cody_middle(double y) { | |
| constexpr std::array p{ | |
| 1230.33935479799725, | |
| 2051.07837782607147, | |
| 1712.04761263407058, | |
| 881.952221241769090, | |
| 298.635138197400131, | |
| 66.1191906371416295, | |
| 8.88314979438837594, | |
| 5.64188496988670089e-1 | |
| }; | |
| constexpr std::array q{ | |
| 1230.33935480374942, | |
| 3439.36767414372164, | |
| 4362.61909014324716, | |
| 3290.79923573345963, | |
| 1621.38957456669019, | |
| 537.181101862009858, | |
| 117.693950891312499, | |
| 15.7449261107098347, | |
| }; | |
| const double y2 = y * y; | |
| std::array<double, 4> p_vals; | |
| std::array<double, 4> q_vals; | |
| for (int i = 0, j = 0; i < 4; i++, j+=2) { | |
| p_vals[i] = std::fma(p[j + 1], y, p[j]); | |
| q_vals[i] = std::fma(q[j + 1], y, q[j]); | |
| } | |
| double num = std::fma(2.15311535474403846e-8, y2, p[3]); | |
| double den = y2 + q[3]; | |
| for (int i = 2; i >= 0; i--) { | |
| num = std::fma(num, y2, p[i]); | |
| den = std::fma(den, y2, q[i]); | |
| } | |
| return num / den; | |
| } |
andrjohns
left a comment
There was a problem hiding this comment.
Thanks Aki! A few more speed-optimisations and accuracy improvements and then I think this is just about ready
|
|
||
| /** | ||
| * Correction factor of the Cody (1969) third-interval rational: | ||
| * `erfcx(x) = (INV_SQRT_PI + u * correction(u)) / x` with `u = 1 / x^2`. | ||
| * | ||
| * Split out from the value so that the derivative can reuse it. See | ||
| * `erfcx_derivative`. | ||
| * | ||
| * @tparam T scalar type | ||
| * @param u inverse square of the argument, `0 <= u <= 1/16` | ||
| * @return `P(u) / Q(u)` | ||
| */ | ||
| template <typename T> | ||
| inline T erfcx_tail_correction(const T& u) { | ||
| static constexpr std::array p{ | ||
| 0.000658749161529837803157, 0.0160837851487422766278, | ||
| 0.125781726111229246204, 0.360344899949804439429, | ||
| 0.305326634961232344035}; | ||
| static constexpr std::array q{ | ||
| -0.00233520497626869185443, -0.0605183413124413191178, | ||
| -0.527905102951428412248, -1.87295284992346047209, | ||
| -2.56852019228982242072}; | ||
| T numerator = 0.0163153871373020978498; | ||
| T denominator = -1.0; | ||
| for (int i = 4; i >= 0; --i) { | ||
| numerator = numerator * u + p[i]; | ||
| denominator = denominator * u + q[i]; | ||
| } | ||
| return numerator / denominator; | ||
| } | ||
|
|
||
| /** | ||
| * Cody (1969) third-interval rational, valid for `x >= 4`. | ||
| * | ||
| * Gives `erfcx` directly. `x * x` is infinite for `x` large enough, which | ||
| * correctly collapses this to the leading term `INV_SQRT_PI / x` and, at | ||
| * infinity, to zero. | ||
| * | ||
| * @param x argument, `x >= 4` | ||
| * @return scaled complementary error function | ||
| */ | ||
| inline double erfcx_cody_tail(double x) { | ||
| const double u = 1.0 / (x * x); | ||
| return (INV_SQRT_PI + u * erfcx_tail_correction(u)) / x; | ||
| } |
There was a problem hiding this comment.
If you declare the coefficients separately in the internal namespace then they can be reused for multiple corrections below.
Also, by parameterising the tail polynomial calculation to work on iterators (from the coefficient arrays) then you can enable another speedup by walking the coefficients in reverse, allowing one division to be dropped:
| /** | |
| * Correction factor of the Cody (1969) third-interval rational: | |
| * `erfcx(x) = (INV_SQRT_PI + u * correction(u)) / x` with `u = 1 / x^2`. | |
| * | |
| * Split out from the value so that the derivative can reuse it. See | |
| * `erfcx_derivative`. | |
| * | |
| * @tparam T scalar type | |
| * @param u inverse square of the argument, `0 <= u <= 1/16` | |
| * @return `P(u) / Q(u)` | |
| */ | |
| template <typename T> | |
| inline T erfcx_tail_correction(const T& u) { | |
| static constexpr std::array p{ | |
| 0.000658749161529837803157, 0.0160837851487422766278, | |
| 0.125781726111229246204, 0.360344899949804439429, | |
| 0.305326634961232344035}; | |
| static constexpr std::array q{ | |
| -0.00233520497626869185443, -0.0605183413124413191178, | |
| -0.527905102951428412248, -1.87295284992346047209, | |
| -2.56852019228982242072}; | |
| T numerator = 0.0163153871373020978498; | |
| T denominator = -1.0; | |
| for (int i = 4; i >= 0; --i) { | |
| numerator = numerator * u + p[i]; | |
| denominator = denominator * u + q[i]; | |
| } | |
| return numerator / denominator; | |
| } | |
| /** | |
| * Cody (1969) third-interval rational, valid for `x >= 4`. | |
| * | |
| * Gives `erfcx` directly. `x * x` is infinite for `x` large enough, which | |
| * correctly collapses this to the leading term `INV_SQRT_PI / x` and, at | |
| * infinity, to zero. | |
| * | |
| * @param x argument, `x >= 4` | |
| * @return scaled complementary error function | |
| */ | |
| inline double erfcx_cody_tail(double x) { | |
| const double u = 1.0 / (x * x); | |
| return (INV_SQRT_PI + u * erfcx_tail_correction(u)) / x; | |
| } | |
| // Cody (1969) third interval, ascending powers of `u = 1 / x^2` | |
| inline constexpr std::array<double, 6> erfcx_tail_p{ | |
| 0.000658749161529837803157, 0.0160837851487422766278, | |
| 0.125781726111229246204, 0.360344899949804439429, | |
| 0.305326634961232344035, 0.0163153871373020978498}; | |
| inline constexpr std::array<double, 6> erfcx_tail_q{ | |
| -0.00233520497626869185443, -0.0605183413124413191178, | |
| -0.527905102951428412248, -1.87295284992346047209, | |
| -2.56852019228982242072, -1.0}; | |
| // Cody (1969) second interval, ascending powers of `y` | |
| inline constexpr std::array<double, 9> erfcx_middle_p{ | |
| 1230.33935479799725, 2051.07837782607147, 1712.04761263407058, | |
| 881.952221241769090, 298.635138197400131, 66.1191906371416295, | |
| 8.88314979438837594, 5.64188496988670089e-1, 2.15311535474403846e-8}; | |
| inline constexpr std::array<double, 9> erfcx_middle_q{ | |
| 1230.33935480374942, 3439.36767414372164, 4362.61909014324716, | |
| 3290.79923573345963, 1621.38957456669019, 537.181101862009858, | |
| 117.693950891312499, 15.7449261107098347, 1.0}; | |
| // `2 * y * P(y) - 2 / sqrt(pi) * Q(y)`, formed at 512 bits | |
| inline constexpr std::array<double, 10> erfcx_middle_dp{ | |
| -1.38828929641828509756e+03, -1.42023212188952842156e+03, | |
| -8.20531739638677179304e+02, -2.89174075427329356508e+02, | |
| -6.56377752033711740945e+01, -8.87368790370412020304e+00, | |
| -5.65021004636098339946e-01, 5.29779936004038675702e-05, | |
| -2.17311817239589615890e-06, 4.30623070948807692000e-08}; | |
| /** | |
| * The two tail polynomials, taking the coefficients in the order the | |
| * iterators give them as descending powers of `t`. | |
| * | |
| * @tparam T scalar type | |
| * @tparam It coefficient iterator | |
| * @param p first numerator coefficient | |
| * @param q first denominator coefficient | |
| * @param t argument | |
| * @return numerator and denominator | |
| */ | |
| template <typename T, typename It> | |
| inline std::pair<T, T> erfcx_tail_polynomials(It p, It q, const T& t) { | |
| T numerator = *p; | |
| T denominator = *q; | |
| for (std::size_t i = 1; i < erfcx_tail_p.size(); ++i) { | |
| numerator = numerator * t + *++p; | |
| denominator = denominator * t + *++q; | |
| } | |
| return {numerator, denominator}; | |
| } | |
| /** | |
| * Correction factor of the Cody (1969) third-interval rational: | |
| * `erfcx(x) = (INV_SQRT_PI + u * correction(u)) / x` with `u = 1 / x^2`. | |
| * | |
| * @tparam T scalar type | |
| * @param u inverse square of the argument, `0 <= u <= 1/16` | |
| * @return `P(u) / Q(u)` | |
| */ | |
| template <typename T> | |
| inline T erfcx_tail_correction(const T& u) { | |
| const auto pq = erfcx_tail_polynomials(erfcx_tail_p.crbegin(), | |
| erfcx_tail_q.crbegin(), u); | |
| return pq.first / pq.second; | |
| } | |
| /** | |
| * Cody (1969) third-interval rational, valid for `x >= 4`. | |
| * | |
| * `u * correction(u)` is formed from the reversed coefficients in `x^2`, so | |
| * its division does not wait on `1 / x^2`. | |
| * | |
| * @param x argument, `x >= 4` | |
| * @return scaled complementary error function | |
| */ | |
| inline double erfcx_cody_tail(double x) { | |
| // correction is below eps / 8, and x^12 would overflow further out | |
| constexpr double leading_term_only = 0x1p27; | |
| if (x > leading_term_only) { | |
| return INV_SQRT_PI / x; | |
| } | |
| const double s = x * x; | |
| const auto pq | |
| = erfcx_tail_polynomials(erfcx_tail_p.cbegin(), erfcx_tail_q.cbegin(), s); | |
| return (INV_SQRT_PI + pq.first / (s * pq.second)) / x; | |
| } |
| /** | ||
| * Derivative of `erfcx`, `2 * x * erfcx(x) - 2 / sqrt(pi)`. | ||
| * | ||
| * That difference cancels for large `x`: both terms approach | ||
| * `2 / sqrt(pi)` while the result decays like `1 / (sqrt(pi) * x^2)`. | ||
| * Measured against a 50-digit reference, the difference form gives 6.1 ulp | ||
| * at `x = 4` and 2.55e+11 ulp at `x = 1e6`. | ||
| * | ||
| * For `x >= 4` the tail rational gives the derivative with no subtraction, | ||
| * because the constant cancels analytically: | ||
| * | ||
| * `2 * x * (INV_SQRT_PI + u * C(u)) / x - 2 / sqrt(pi) = 2 * u * C(u)` | ||
| * | ||
| * since `2 * INV_SQRT_PI` is `2 / sqrt(pi)`. That form measures 0.1 to | ||
| * 31 ulp over the same range. | ||
| * | ||
| * @tparam T scalar type | ||
| * @param x argument | ||
| * @param value `erfcx(x)` | ||
| * @return derivative of `erfcx` at `x` | ||
| */ | ||
| template <typename T> | ||
| inline T erfcx_derivative(const T& x, const T& value) { | ||
| if (x >= 4.0) { | ||
| const T u = 1.0 / (x * x); | ||
| return 2.0 * u * erfcx_tail_correction(u); | ||
| } | ||
| return 2.0 * x * value - TWO_OVER_SQRT_PI; | ||
| } |
There was a problem hiding this comment.
The gradient calculations are still losing precision in the [0.46875, 4) and >=30 ranges, and can be improved with additional series (at a slight performance cost).
For these series (and the middle-tail calcs) can speed up the evaluation by running over y^2 and evaluating the polynomials in pairs:
y2 = y*y
p = c8
p = p*y2 + (c6 + c7*y)
p = p*y2 + (c4 + c5*y)
p = p*y2 + (c2 + c3*y)
p = p*y2 + (c0 + c1*y)
| /** | |
| * Derivative of `erfcx`, `2 * x * erfcx(x) - 2 / sqrt(pi)`. | |
| * | |
| * That difference cancels for large `x`: both terms approach | |
| * `2 / sqrt(pi)` while the result decays like `1 / (sqrt(pi) * x^2)`. | |
| * Measured against a 50-digit reference, the difference form gives 6.1 ulp | |
| * at `x = 4` and 2.55e+11 ulp at `x = 1e6`. | |
| * | |
| * For `x >= 4` the tail rational gives the derivative with no subtraction, | |
| * because the constant cancels analytically: | |
| * | |
| * `2 * x * (INV_SQRT_PI + u * C(u)) / x - 2 / sqrt(pi) = 2 * u * C(u)` | |
| * | |
| * since `2 * INV_SQRT_PI` is `2 / sqrt(pi)`. That form measures 0.1 to | |
| * 31 ulp over the same range. | |
| * | |
| * @tparam T scalar type | |
| * @param x argument | |
| * @param value `erfcx(x)` | |
| * @return derivative of `erfcx` at `x` | |
| */ | |
| template <typename T> | |
| inline T erfcx_derivative(const T& x, const T& value) { | |
| if (x >= 4.0) { | |
| const T u = 1.0 / (x * x); | |
| return 2.0 * u * erfcx_tail_correction(u); | |
| } | |
| return 2.0 * x * value - TWO_OVER_SQRT_PI; | |
| } | |
| /** | |
| * Horner in `y^2` over adjacent coefficient pairs. | |
| * | |
| * @tparam T scalar type | |
| * @param c coefficients, ascending | |
| * @param y argument | |
| * @param y2 `y * y` | |
| * @return polynomial value | |
| */ | |
| template <typename T, std::size_t N> | |
| inline T erfcx_paired_horner(const std::array<double, N>& c, const T& y, | |
| const T& y2) { | |
| int j = static_cast<int>(N) - (N % 2 ? 3 : 4); | |
| T r = N % 2 ? T(c[N - 1]) : T(c[N - 2] + c[N - 1] * y); | |
| for (; j >= 0; j -= 2) { | |
| r = r * y2 + (c[j] + c[j + 1] * y); | |
| } | |
| return r; | |
| } | |
| /** | |
| * Derivative of `erfcx`, `2 * x * erfcx(x) - 2 / sqrt(pi)`. | |
| * | |
| * That difference cancels for `x > 1`, so there the subtraction is done | |
| * analytically instead: in the coefficients on the middle interval, against | |
| * the leading term of the tail rational above 4, and by the asymptotic series | |
| * above 30, where the tail rational's own fit error dominates. | |
| * | |
| * @tparam T scalar type | |
| * @param x argument | |
| * @param value `erfcx(x)` | |
| * @return derivative of `erfcx` at `x` | |
| */ | |
| template <typename T> | |
| inline T erfcx_derivative(const T& x, const T& value) { | |
| if (x < 0.46875) { | |
| return 2.0 * x * value - TWO_OVER_SQRT_PI; | |
| } | |
| const T x2 = x * x; | |
| if (x < 4.0) { | |
| return erfcx_paired_horner(erfcx_middle_dp, x, x2) | |
| / erfcx_paired_horner(erfcx_middle_q, x, x2); | |
| } | |
| const T u = 1.0 / x2; | |
| if (x < 30.0) { | |
| return 2.0 * u * erfcx_tail_correction(u); | |
| } | |
| // `-sqrt(pi) * x^2 * erfcx'(x)`, asymptotic, ascending powers of `u` | |
| static constexpr std::array<double, 8> series_coefficients{ | |
| 1.0, -1.5, 3.75, -13.125, | |
| 59.0625, -324.84375, 2111.484375, -15836.1328125}; | |
| const T series = erfcx_paired_horner(series_coefficients, u, u * u); | |
| return -INV_SQRT_PI * u * series; | |
| } |
| inline double erfcx_cody_middle(double y) { | ||
| double p = 2.15311535474403846e-8 * y; | ||
| p = (p + 5.64188496988670089e-1) * y; | ||
| p = (p + 8.88314979438837594) * y; | ||
| p = (p + 66.1191906371416295) * y; | ||
| p = (p + 298.635138197400131) * y; | ||
| p = (p + 881.952221241769090) * y; | ||
| p = (p + 1712.04761263407058) * y; | ||
| p = (p + 2051.07837782607147) * y; | ||
| double q = y; | ||
| q = (q + 15.7449261107098347) * y; | ||
| q = (q + 117.693950891312499) * y; | ||
| q = (q + 537.181101862009858) * y; | ||
| q = (q + 1621.38957456669019) * y; | ||
| q = (q + 3290.79923573345963) * y; | ||
| q = (q + 4362.61909014324716) * y; | ||
| q = (q + 3439.36767414372164) * y; | ||
| return (p + 1230.33935479799725) / (q + 1230.33935480374942); | ||
| } |
There was a problem hiding this comment.
After declaring the coefficients separately and specifying the paired-horner implementation function, this can simplify down to:
| inline double erfcx_cody_middle(double y) { | |
| double p = 2.15311535474403846e-8 * y; | |
| p = (p + 5.64188496988670089e-1) * y; | |
| p = (p + 8.88314979438837594) * y; | |
| p = (p + 66.1191906371416295) * y; | |
| p = (p + 298.635138197400131) * y; | |
| p = (p + 881.952221241769090) * y; | |
| p = (p + 1712.04761263407058) * y; | |
| p = (p + 2051.07837782607147) * y; | |
| double q = y; | |
| q = (q + 15.7449261107098347) * y; | |
| q = (q + 117.693950891312499) * y; | |
| q = (q + 537.181101862009858) * y; | |
| q = (q + 1621.38957456669019) * y; | |
| q = (q + 3290.79923573345963) * y; | |
| q = (q + 4362.61909014324716) * y; | |
| q = (q + 3439.36767414372164) * y; | |
| return (p + 1230.33935479799725) / (q + 1230.33935480374942); | |
| } | |
| inline double erfcx_cody_middle(double y) { | |
| const double y2 = y * y; | |
| return erfcx_paired_horner(erfcx_middle_p, y, y2) | |
| / erfcx_paired_horner(erfcx_middle_q, y, y2); | |
| } |
| inline double erfcx_small(double x) { | ||
| // Split into the even and odd powers of x, so the two Horner chains run | ||
| // independently. A single degree-18 chain is 18 dependent operations; two | ||
| // chains of 9 halve that latency. | ||
| static constexpr std::array even = {1.0, | ||
| 4.99999999999992839e-01, | ||
| 1.66666666667239644e-01, | ||
| 4.16666666458337179e-02, | ||
| 8.33333374332981443e-03, | ||
| 1.38888415444527033e-03, | ||
| 1.98445679338826757e-04, | ||
| 2.46655529768908249e-05}; | ||
| static constexpr std::array odd | ||
| = {-1.12837916709551256, -7.52252778063651983e-01, | ||
| -3.00901111227312890e-01, -8.59717459974174147e-02, | ||
| -1.91048337772546720e-02, -3.47359067853470795e-03, | ||
| -5.34506929034156810e-04, -7.08163358203131886e-05}; | ||
| static_assert(even.size() == odd.size(), | ||
| "the two chains are stepped by one loop"); | ||
| // The leading coefficient of each chain seeds the accumulator, and the | ||
| // trailing 1.0 is added at the end, so both arrays hold the interior | ||
| // coefficients only. Derive the bound from the array rather than writing | ||
| // it out, so shortening an array cannot leave the loop reading past it. | ||
| const double x2 = x * x; | ||
| double e = 3.05977060678449757e-06; | ||
| double o = -9.35890030086883823e-06; | ||
| for (int i = static_cast<int>(even.size()) - 1; i >= 0; --i) { | ||
| e = even[i] + x2 * e; | ||
| o = odd[i] + x2 * o; | ||
| } | ||
| return 1.0 + x2 * e + x * o; | ||
| } |
There was a problem hiding this comment.
We can parallelise/split the evaluation even more by using Estrin's method instead of Horner's:
| inline double erfcx_small(double x) { | |
| // Split into the even and odd powers of x, so the two Horner chains run | |
| // independently. A single degree-18 chain is 18 dependent operations; two | |
| // chains of 9 halve that latency. | |
| static constexpr std::array even = {1.0, | |
| 4.99999999999992839e-01, | |
| 1.66666666667239644e-01, | |
| 4.16666666458337179e-02, | |
| 8.33333374332981443e-03, | |
| 1.38888415444527033e-03, | |
| 1.98445679338826757e-04, | |
| 2.46655529768908249e-05}; | |
| static constexpr std::array odd | |
| = {-1.12837916709551256, -7.52252778063651983e-01, | |
| -3.00901111227312890e-01, -8.59717459974174147e-02, | |
| -1.91048337772546720e-02, -3.47359067853470795e-03, | |
| -5.34506929034156810e-04, -7.08163358203131886e-05}; | |
| static_assert(even.size() == odd.size(), | |
| "the two chains are stepped by one loop"); | |
| // The leading coefficient of each chain seeds the accumulator, and the | |
| // trailing 1.0 is added at the end, so both arrays hold the interior | |
| // coefficients only. Derive the bound from the array rather than writing | |
| // it out, so shortening an array cannot leave the loop reading past it. | |
| const double x2 = x * x; | |
| double e = 3.05977060678449757e-06; | |
| double o = -9.35890030086883823e-06; | |
| for (int i = static_cast<int>(even.size()) - 1; i >= 0; --i) { | |
| e = even[i] + x2 * e; | |
| o = odd[i] + x2 * o; | |
| } | |
| return 1.0 + x2 * e + x * o; | |
| } | |
| inline double erfcx_small(double x) { | |
| static constexpr std::array<double, 19> c{1.0, | |
| -1.12837916709551256, | |
| 1.0, | |
| -7.52252778063651983e-01, | |
| 4.99999999999992839e-01, | |
| -3.00901111227312890e-01, | |
| 1.66666666667239644e-01, | |
| -8.59717459974174147e-02, | |
| 4.16666666458337179e-02, | |
| -1.91048337772546720e-02, | |
| 8.33333374332981443e-03, | |
| -3.47359067853470795e-03, | |
| 1.38888415444527033e-03, | |
| -5.34506929034156810e-04, | |
| 1.98445679338826757e-04, | |
| -7.08163358203131886e-05, | |
| 2.46655529768908249e-05, | |
| -9.35890030086883823e-06, | |
| 3.05977060678449757e-06}; | |
| const double x2 = x * x; | |
| const double x4 = x2 * x2; | |
| const double x8 = x4 * x4; | |
| std::array<double, 4> quad; | |
| for (int i = 0; i < 4; ++i) { | |
| const int j = 4 * i + 2; | |
| quad[i] = (c[j] + c[j + 1] * x) + (c[j + 2] + c[j + 3] * x) * x2; | |
| } | |
| const double rest | |
| = (quad[0] + quad[1] * x4) + (quad[2] + quad[3] * x4 + c[18] * x8) * x8; | |
| return (c[0] + c[1] * x) + x2 * rest; | |
| } |
|
@SteveBronder - using CPUs without fma instructions would be pretty old (~2010), how risky do you think it would be to enable by default (adding a |
|
The time zone difference causes some asynchronity issues in this PR discussion. I had gone through @SteveBronder's latest suggestions, but had not yet pushed new commits as it was too late yesterday. I have solution for fma issue. I usually would use -march, but one of the clusters I'm using has heteregenous nodes and -march causes problems. |
I think everything post haswell is going to have an fma instruction. But making sure users have the flag set is very real. We can just use a * b + c since that should still give the fma instruction when march is used |
|
@andrjohns with your suggestions the speed drops at least on CPU Value and derivative together, which is what the reverse pass costs
So the speed does drop, by 18–33 % on the AD path, in exchange for the accuracy:
What is your preference? |
Ahh yeah that's really not worth it, safe to ignore me |
|
All right, took some time to test the effect of each suggestion on speed and accuracy. The decision was to prioritize speed as long as the accuracy is still good. Even 100 ulp is much better than what we had (10^9 to 10^11 ulp). Speed is the ratio against the code before your suggestions and the version with changes listed below, value and derivative together. A ratio above 1 is faster. Xeon E5-2680 v3 and Tesla V100-SXM2, three repetitions, ns per call. CPU.
Value alone: 1.56 on [5, 20], 1.41 on [0.125, 12], 1.38 on [-0.46875, 0.46875], 1.13 on [-4, 4], 1.05 on [-20, 0]. GPU. The device function keeps its own evaluation order, because every candidate was measured on a V100 first and only one of them won there:
So the mathematics is the same on both targets and the evaluation order is not: pairing and Estrin are 1.4–1.8× on the Xeon and 8–13 % slower on the V100. Accuracy. Worst case per range over 4001 points, against mpmath at 40 digits. Here 1 ulp means a relative error of
The value stays within 4.3 ulp, or 4.8e-16 relative, everywhere. The derivative is worse than the value in two places, for two separate reasons. On On Absolute error is not a useful measure on the negative axis: What went in. From @SteveBronder: the paired Cody middle coefficients on CPU, What did not. The middle One change is from neither of you. With the reversed walk the value evaluates the tail rational in Three changes from the suggestions as written. In the 20:03 All suites pass on a V100: |
|
I'm investigating one additonal simple idea to reduce the worst ulp (even it is already good) |
|
I was not able to reduce the worst ulp 192 of
The Mills-ratio route never touches The AD route damps The error is large in relative terms precisely where the derivative is small, so it never reaches the answer. For other distributions calling So the honest conclusion for the PR: the 192 ulp is invisible to the intended consumer. Reducing that 192 ulp is possible, but costs 0.72 on a interval. I think we can stop trying to optimize this further (at least for now) |
Jenkins Console Log Machine informationDistributor 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.000 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 |
Claude and Herbie assisted work. Re-uses existing code from Stan math normal_lcdf.hpp
Summary
Implements #3401: adds
erfcx(x) = exp(x²)·erfc(x)toprim,rev,fwdand the OpenCL kernel generator.#3401 covers the motivation. The final choice of algorithm is a bit different. This text describes the algorithm as it now stands, the measurements behind it, and how it relates to the approximation already in
normal_lcdf.Algorithm
Four branches. Cody (1969) gives rational approximations on three intervals, and on the middle one it is arranged as
erfc(x) = exp(-x²)·R(x), soerfcx(x) = R(x)with the exponential cancelled analytically — noexpand noerfccall.erfcis 66 % of the cost of formingexp(x²)·erfc(x)(27.0 ns of 41.1 on[0, 4], measured), anderfcitself computesexp(-x²)internally, which we would then cancel with our ownexp(x²). Two exponentials to produce a product that analytically has none.x >= 40.46875 <= x < 4|x| < 0.46875x <= -0.468752·exp(x²) - erfcx(-x)expThe whole positive axis is computed without a single library call. On the negative side
expis unavoidable, sinceerfcxgrows like2·exp(x²)asx → -∞; belowx = -6.1,erfcx(-x)is undereps/2of2·exp(x²), so the subtraction is skipped.The small-window expansion. The plain Maclaurin series for
erfcxis only usable to about|x| = 0.125. Four further Chebyshev-economized terms extend it to 0.46875 at no measurable cost on CPU (11.00 vs 11.05 ns), which is what removes the lastexpfrom the positive axis. It also covers both signs with no branch. Its low-order coefficients reproduce the Maclaurin coefficients oferfcxexactly (1,-2/√π,1,-0.752252778063652,0.5, ...), which is a useful independent check that the fit is right.fma, not a Dekker split, on the negative branch. The rounding error ofx²is recovered withfma(x, x, -h)and folded back in, becauseexpamplifies it into roughlyx²·eps— 512 ulp atx = -26if left alone. An earlier revision of this PR used a Veltkamp/Dekker split instead. That is correct on the host but silently wrong on a GPU:t - (t - x)isxover the reals but not in floating point, and NVIDIA's OpenCL compiler simplifies it. Measured on a Tesla V100, the split form returnedx_lo == 0at every one of 4096 test points and scored 512.66 ulp — bit-identical to the uncorrected formula it existed to avoid.fmais a single instruction with nothing to reassociate.Note that no test in this repository could have caught that: the erfcx tests are CPU-only and CI has no GPU, so the device function is never executed against reference values anywhere. It was found by running the shipped kernel string on real hardware. That is a general gap in OpenCL coverage rather than something specific to this PR.
Relation to
normal_lcdfprim/prob/normal_lcdf.hppalready carries the Cody third-interval coefficients for its cdf value, with the same crossover at 4 and for the same reason: it is R'spnormcutoff|y| > M_SQRT_32expressed in theerfcargument,sqrt(32)/sqrt(2) = 4exactly. That part is shared.Everything else is different:
normal_lcdfuses only the tail rational, and covers the interior withexp/erfcplus, for the gradient, five hand-placed Taylor expansions and four residual fits (~350 lines). Those interior cutoffs are not from the literature; they were derived empirically in Bugfix/issue #1284 numerical precision of normal lcdf #1411 to patch regions where the autodiff tester was failing.erfcxuses all three Cody intervals plus a Chebyshev expansion, and needs no Taylor patches anywhere, because it never formsexp(+x²)on the positive axis at all.That difference makes this function better.
normal_lcdf's own documentation records its worst in-range relative gradient error as 7.6e-06 to 6.1e-05, "just inside the 1e-4 relative gradient toleranceexpect_adapplies by default". Measured against mpmath, its gradient is 2.7e+09 to 4.2e+11 ulp, where theerfcxidentity√(2/π)/erfcx(-x/√2)gives 3 to 7 ulp — eight to eleven orders of magnitude better. SciPy'slog_ndtrtakes the same route and needs no rational approximation or Taylor patches belowx = -1for exactly this reason. Thenormal_lcdfimprovement usingerfcxis a separate PR.Accuracy
Worst ulp against an
mpmathreference at 40 decimal digits, 40 001 points per range, CPU and the OpenCL device function scored against the same reference. libcerf 3.8 is included as the reference implementation of this function.libcerf is more accurate on
[0, 4]and this PR is more accurate in both tails. These accuracy differences are negligible. The currentstd_normal_lcdfgradient is 10⁹–10¹¹ ulp. That's the codeerfcxexists to replace, and the erfcx route gives 3–7 ulp — an improvement of 8 to 11 orders of magnitude.Speed
CPU: ns/call on one dedicated core of a Xeon E5-2680 v3 (Haswell, 2.50 GHz,
performancegovernor),-O3, 4096 pseudo-random arguments × 8000 repetitions, three repetitions, spread under 1 %.GPU: ns/call on a Tesla V100-SXM2 (OpenCL 3.0 CUDA, 80 CUs), 1 048 576 work-items × 2000 in-register iterations, harness overhead subtracted. A simple map kernel would be bandwidth-bound and would measure the memory system rather than the function, so each work-item walks the range with an in-register PRNG; that also gives neighbouring lanes genuinely different subintervals, which is the case of interest.
On CPU the two are close: this PR is 1.17× faster on
[0, 4]and 1.07× on[−4, 4], libcerf is 1.04–1.17× faster elsewhere. On GPU this PR is 12× to 68× faster.The GPU gap is the reason the table-based approach was not adopted. libcerf indexes a 16 640-byte
__constanttable at a position derived from the input, so lanes within a warp hit different subdomains and the read serialises. This PR is pure ALU with no memory traffic. The CPU numbers show the table costs nothing there — it is a GPU-specific penalty, and a large one.Absolute cost is strongly machine-dependent: the same CPU benchmark on a 13th-gen laptop part runs about 2.5× faster and gives noticeably different ratios, so the ratios matter more than the absolute numbers.
Edge cases:
erfcx(0) = 1,erfcx(+inf) = 0,erfcx(-inf) = +inf, NaN preserved,erfcx(1e10)collapses to the leading term1/(x√π),erfcx(-26.9) = inf(the true value has already left binary64).The derivative reuses the value:
so
revandfwdneed no secondexporerfccall, the gradient is exact wherever the value is, and the autodiff layers are independent of which value algorithm sits underneath.Tests
Four files, 18 tests.
test/unit/math/prim/fun/erfcx_test.cpp— values against the reference across negative, small, moderate and far-tail arguments;erfcxUpperTailassertsstd::erfc(30.0) == 0alongside a finiteerfcx, i.e. the failure this function exists to avoid;erfcxBranchContinuitychecks all three internal crossovers (0.46875, 4, −0.46875) against the analytic slope, so a future change to any cutoff cannot silently introduce a step;erfcxReflectionSkippedTermasserts that the term skipped belowx = -6.1really is under half an ulp of2·exp(x²);erfcxEdgeCases;erfcxVectorized; anderfcxNormalTailIdentities, which verifies that the two identities in Adderfcx(scaled complementary error function) as a Stan Math primitive #3401 reproducestd_normal_lcdfand the inverse Mills ratio.test/unit/math/rev/fun/erfcx_test.cpp— value and gradient against2x·erfcx(x) − 2/√π, both theEigen::Matrix<var, -1, 1>andvar_value<Eigen::VectorXd>overloads, and an upper-tail gradient atx = 30where the unscaled route cannot produce one.test/unit/math/fwd/fun/erfcx_test.cpp—fvar<double>value/tangent, chain rule,fvar<fvar<double>>second order.test/unit/math/mix/fun/erfcx_test.cpp—expect_adandexpect_ad_vector_matvar, with arguments straddling the crossovers.The mix test uses explicit finite arguments rather than
expect_common_unary_vectorized, becauseerfcx(-inf) = +infand the function grows like2exp(x²)to the left, so the common argument set is not finite-differenceable here.erfcxReflectionSkippedTermdeliberately does not sweep for continuity atx = -6.1:erfcxhas slope about-3.5e17there, so the function's own variation swamps any fixed relative tolerance. It asserts the design claim directly instead.Separately, the Horner chains in both the host and device functions are written flat (
p = c + x * p;repeated) rather than nested, because cpplint's line-length rule rejects deep nesting. Flattening is an easy place to introduce an off-by-one, so the flattened forms were checked bit-identical against the loop forms they came from over 500 000 points per interval.Side Effects
Adds one include to each of the
prim/rev/fwdfun.hppaggregators and a kernel-generator registration; no existing function changes behaviour.Release notes
Added
erfcx(x), the scaled complementary error functionexp(x^2)*erfc(x), with reverse- and forward-mode derivatives and OpenCL support.Checklist
Copyright holder: Aalto University
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
./runTests.py test/unit)make test-headers)make test-math-dependencies)make doxygen)make cpplint)the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested