Reproducible exp() function - #95
Merged
Merged
Conversation
* re_end was looking for "end module|procedure", but this did not
account for submodules or other blocks.
This was modified to a general "block" re_end:
"^ *end *(module|submodule|program)
Not sure if it's more robust but it fixed a problem with incorrect
handling of module procedure in submodules.
* link_obj was updated to follow `use` statements from submodules.
This includes new macros to support LDFLAGS extensions for OPT and COVERAGE builds (which also happen to be used by the timing and unit test builds). It also forces rebuild of the ./configure script on rebuilds. Since it is now a more semi-permanent file that is a product of multiple builds, it is more important to keep it in sync across different Makefiles and builds.
This adds a macro to test if parentheses are protected, as needed by the fast_rint() function to appear in a subsequent commit for exp_repro(). The naming suggests a much more narrow scope at the moment, but we could generalize it in the future if appropriate.
The current TestSuite would do aggressive PE syncs, even if there is only one core, and even if MPI was never initialized. This led to problems. For now, this patch wraps those syncs in an initialization flag check, but there is probably a better way to approach this in the future.
This function offers a reproducible alternative to the exp() intrinsic,
whose implementation is inherently ambiguous.
As an outline of the implementation:
* A range-reduction is applied from x to r = x - K ln2 - i/N ln2 where
K = nint(x/ln2) and i is a subdivision within the interval [-ln2/2, ln2/2].
This reduces the problem to exp(x) = 2**K 2**(i/N) exp(r).
* exp(r) is estimated using a Remez minimax polynomial, optimized to the
subinterval range [-ln2/2N, ln2/2N].
2**(i/N) is obtained from a hard-coded lookup table.
* The final result is computed by applying the exact 2**K scaling, along
with additional steps to account for subnormal values.
Results are identical across Intel, GCC, NVIDIA CPU and NVIDIA GPU over
several tested ranges: O(1M) points between -1:1, -10:10, and -700:700.
Results are nearly correctly-rounded, with almost all below 0.5 ULP.
Typical estimates are shown below.
=== scalar exp_repro() accuracy
Tested 10000000 points in [-10, 10]
max abs err: 1.83371E-12 at x = 9.9062
max rel err: 1.11419E-16 at x = 5.5476
max ULP err (vs quad): 0.5093838179 at x = 1.9410
mean abs err: 4.49703E-14
mean rel err: 4.00488E-17
RMS err: 1.70446E-13
correct (<0.5 ULP): 9989058 ( 99.89%)
above 0.5 ULP: 10942 ( 0.11%)
above 1 ULP: 0 ( 0.00%)
FMAs produce different answers, but comparable accuracy.
=== scalar exp_repro() accuracy
Tested 10000000 points in [-10, 10]
max abs err: 1.82849E-12 at x = 9.7173
max rel err: 1.11419E-16 at x = 5.5476
max ULP err (vs quad): 0.5079601302 at x = 2.6019
mean abs err: 4.49838E-14
mean rel err: 4.00434E-17
RMS err: 1.70520E-13
correct (<0.5 ULP): 9991443 ( 99.91%)
above 0.5 ULP: 8557 ( 0.09%)
above 1 ULP: 0 ( 0.00%)
=== vector exp_repro() matches scalar
Nonfinite numbers (Inf, NaN) are respected, and IEEE signals match
expected results.
IEEE flags summary: exp() exp_repro
IOUXZ IOUXZ
exact (0): ..... .....
normal: ...X. ...X.
overflow: .O.X. .O.X.
underflow: ..UX. ..UX.
near overflow: ...X. ...X.
near underflow: ...X. ...X.
largest float: ...X. ...X.
smallest normal: ...X. ...X.
+Inf: ..... .....
-Inf: ..... .....
NaN: ..... .....
sNaN: I.... I....
Performance is slower than peak vectorized intrinsics, although accuracy
is far greater in these cases. ifx 2025.2 results with -O3 -xHost and
-ipo are shown below.
=== MOM_intrinsic_functions timing ===
npts = 100000, niter = 200
x range: [ -10.0, 10.0]
exp() time/elem: 0.79 ns
exp_repro() time/elem: 1.22 ns
slowdown factor: 1.55x
=== scalar loop-carried timing ===
baseline scalar time/call: 8.29 ns
exp() scalar time/call: 16.13 ns
exp_repro() scalar time/call: 20.87 ns
scalar slowdown factor: 1.29x
exp() minus baseline: 7.83 ns
exp_repro() minus baseline: 12.58 ns
adjusted slowdown factor: 1.61x
Timing difference will be much lower in compilers which cannot inline.
Having said that, most exp() calls will not be vectorized, and
This does not implement a method for selecting an exp() implementation.
That will (presumably) come in a future commit.
Several classes of unit and timing tests for exp_repro(). * exp(0.) = 1. exact test * Several tolerance tests * Property tests: * exp(a+b) = exp(a)*exp(b) * exp(-x) = 1./exp(x) * Nonfinites: * exp(-Inf) = +0. * exp(x)=x for +Inf,+/-NaN * Subnormal evaluation * ULP-accurate measurement * Floating-point signal correctness * Timing comparisons to exp() for scalar and (vectorized) arrays
GPU Port CoverageOverall: 2176 / 5577 portable executed lines ported (39.0%) Files touched by this PR: 0 / 6 portable executed lines ported (0.0%) Full per-file / per-routine breakdown: see the "gpu-port-report" job summary and artifact. |
edoyango
approved these changes
Sep 2, 2026
edoyango
left a comment
edoyango
left a comment
Collaborator
There was a problem hiding this comment.
Tested with set_viscous_ML and it works well! To avoid team size of 1, need to add bind(teams,parallel) to openmp loop directives. Performance isn't noticeably impacted, but register usage is a bit higher.
Member
Author
|
I've been using this for testing: we might want to preprocess this statement somehow to swap between the two, if it impacts register consumption. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Still under review at GFDL (NOAA-GFDL#1164) but seems unlikely to change at this point.
Any changes will probably be new PRs so this should be OK to take in now.