exp_repro: Reproducible implementation of exp() by marshallward · Pull Request #1164 · NOAA-GFDL/MOM6

marshallward · 2026-08-18T23:03:53Z

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(10M) 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

The Intel vs exp_repro() accuracy can be visualized below. (Not GCC, despite the title...)

If -fp-model source is enabled then Intel's vector exp() accuracy is restored, although performance drops significantly (see below).

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

With -fp-model source vectorized performance drops significantly:

=== MOM_intrinsic_functions timing ===
npts = 100000, niter = 200
x range: [ -10.0,   10.0]
 
exp() time/elem:                 3.53 ns
  sum (to prevent elision):   1.10142E+08
exp_repro() time/elem:           1.61 ns
  sum (to prevent elision):   1.10142E+08
 
slowdown factor:                 0.46x

Timings and unit tests have been included:

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-accuracy measurement
Floating-point signal correctness
Timing comparisons to exp() for scalar and (vectorized) arrays

This PR does not implement a method for selecting an exp() implementation.
That will (presumably) come in a future commit.

Changes to suppor this PR:

TestSuite was modified to run more easily without MPI
.testing/Makefile added hooks for OPT and COVERAGE builds
A new autoconf macro MOM6_RC_FAST_RINT to replace ieee_rint() with a faster alternative when possible.
./configure is now correctly force-built by .testing/Makefile
makedep was updated so that end procedure did not prematurely end a module parse. Links from submodules were also cleaned up.

adcroft · 2026-08-31T13:47:38Z

I have one question: you showed a 1.55x slow down for exp() itself, but what is the impact on the model going to be, and how will we invoke it in the model?

adcroft

Should we consider making the extreme value handling (e.g. underflow) optional, which I understand would make this faster than the library versions. In the reproducing sums, we take advantage of knowing the range we care about, and we could do the same here.

marshallward · 2026-08-31T20:28:08Z

@adcroft I have an experimental branch which adds an add_scalar() test type.

https://github.com/marshallward/MOM6/blob/fcdf4f8d960729ab283b375531cbc6a1020032fc/src/framework/testing/MOM_intrinsic_functions_tests.F90#L671

marshallward · 2026-09-01T20:57:31Z

In a more realistic case (benchmark), numbers were not quite as promising as in the
timing tests, but still tolerable. For a 1-day benchmark run:

Intel exp (__libm_exp_z0) is 0.09% (8.8M cycles)
exp_repro is 0.18% (17.8M cycles)

so about 2x slower. This was without IPO or other inlining.

marshallward · 2026-09-02T17:52:57Z

Added changes

Consolidated floating point description parameters from MOM_intrinsic_functions.F90 and MOM_exp.F90 into MOM_intrinsic_functions.F90.
- Updated cuberoot to use the renamed parametesr
Fixed a bug in testing which tried to run ULP testing when quad precision is unavaiable. Now it prints a message saying it could not be run.
MOM_exp_data_n128.F90 content now explicitly public
[nondim] added all over

Copilot

🟡 Changes recommended

There are concrete build/test correctness issues (Python IndentationError in ac/makedep, non-deterministic “deterministic” RNG in property tests, and an always-true real128-availability gate) that should be fixed before approval.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Pull request overview

Reminder: MOM6 has a policy on AI-assisted contributions in Consortium-policy-on-AI.md.

This PR adds a new bitwise-reproducible exp_repro() implementation (via a submodule and lookup-table + polynomial approach), along with unit tests, timing drivers, and build-system support to enable a fast rounding path when the compiler honors parentheses in FP expressions.

Changes:

Introduces exp_repro() in MOM_intrinsic_functions with the implementation in a new MOM_exp submodule plus lookup-table data.
Adds a comprehensive unit-test module and new standalone unit/timing driver programs for validating accuracy, IEEE flags, and performance.
Updates testing/build infrastructure (MPI-skipping hooks, autoconf macro, makedep parsing, and .testing Makefile options) to support the new tests and reproducibility requirements.

File summaries

File	Description
src/framework/testing/MOM_intrinsic_functions_tests.F90	New unit test suite covering correctness, properties, ULP checks, and IEEE flags for `exp_repro()`
src/framework/MOM_unit_testing.F90	Allows unit tests to run without MPI sync when MPI-dependent behavior is explicitly skipped
src/framework/MOM_intrinsic_functions.F90	Exposes `exp_repro()` via interface and updates FP-layout constants used by existing intrinsics
src/framework/MOM_error_handler.F90	Adds `query_skip_mpi()` accessor to support conditional MPI-dependent behavior
src/framework/intrinsics/MOM_exp.h	Preprocessor macro wrapper selecting `fast_rint()` vs `ieee_rint()` for range reduction
src/framework/intrinsics/MOM_exp.F90	New submodule implementing reproducible `exp_repro()` and supporting polynomial/rounding routine
src/framework/intrinsics/MOM_exp_data_n128.F90	New lookup-table data for `2**(i/ndiv)` scaling and residual correction (ndiv=128)
config_src/drivers/unit_tests/test_MOM_intrinsic_functions.F90	New unit-test driver program for intrinsic-function tests without MPI
config_src/drivers/timing_tests/time_MOM_intrinsic_functions.F90	New timing driver comparing intrinsic `exp()` vs `exp_repro()`
ac/makedep	Updates Fortran parsing/link logic to handle submodules and avoid premature module termination
ac/m4/mom6_fc_fast_rint.m4	New autoconf macro to detect compiler flags needed to protect FP parentheses (enabling fast rint path)
ac/configure.ac	Hooks in the new fast-rint detection macro and adds discovered flags to `FCFLAGS`
.testing/Makefile	Adds OPT linker flags and ensures `configure` is always rebuilt when invoked from `.testing`

Review details

Suppressed comments (2)

ac/makedep:40

These indented comment lines are at top-level indentation and will cause an IndentationError in Python. Top-level comments must not be indented unless inside a block.

    # NOTE: re_procedure excludes comments and tokens with substrings
    # containing `function` or `subroutine`, but will fail if the keywords
    # appear in other contexts.

src/framework/testing/MOM_intrinsic_functions_tests.F90:672

This 'deterministic' generator uses real arithmetic with a 31-bit modulus and a large multiplier; the intermediate products exceed the exact-integer range of IEEE real64, so the generated values are not guaranteed to be reproducible across compilers/platforms.

  ! Use a simple deterministic sequence for reproducibility
  seed = 0.987654321

  do i = 1, npts
    ! Generate pseudo-random values in a range that avoids overflow/underflow
    ! Keep x in [-300, 300] so both exp(x) and exp(-x) are representable
    seed = mod(seed * 1103515245. + 12345., 2.**31)
    x = (seed / 2.**31) * 600. - 300.

Files reviewed: 13/13 changed files
Comments generated: 5
Review effort level: Lite

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

codecov · 2026-09-03T15:07:15Z

Codecov Report

❌ Patch coverage is 95.59118% with 22 lines in your changes missing coverage. Please review.
✅ Project coverage is 38.57%. Comparing base (7d4202b) to head (7af9fab).

Files with missing lines	Patch %	Lines
...ramework/testing/MOM_intrinsic_functions_tests.F90	95.53%	5 Missing and 15 partials ⚠️
src/framework/MOM_unit_testing.F90	60.00%	0 Missing and 2 partials ⚠️

Additional details and impacted files

@@             Coverage Diff              @@
##           dev/gfdl    #1164      +/-   ##
============================================
+ Coverage     38.27%   38.57%   +0.30%     
============================================
  Files           275      278       +3     
  Lines         93940    94431     +491     
  Branches      18164    18191      +27     
============================================
+ Hits          35953    36425     +472     
- Misses        51218    51222       +4     
- Partials       6769     6784      +15

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

* 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. For safety, fast_rint() is disabled for cross-compile builds, although we may want to provide a way to manually override this in the future.

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

Several integers describing the floating point format were defined in both MOM_intrinsic_functions.F90 and MOM_exp.F90, sometimes with different expressions (but hopefully identical values). This patch moves the MOM_exp.F90 definitions to MOM_intrinsic_functions.F90 and updates the cuberoot solver to use the new values. ULP accuray tests are now correctly disabled if quad precision is undetected by the compiler. (Non)dimensional metadata has also been added to several variable docstrings.

adcroft

Looks good. Testing at https://gitlab.gfdl.noaa.gov/ogrp/mom6ci/MOM6/-/pipelines/31864

marshallward force-pushed the exp-repro branch from e720332 to b63251d Compare August 19, 2026 14:24

adcroft reviewed Aug 31, 2026

View reviewed changes

marshallward mentioned this pull request Sep 1, 2026

Reproducible exp() function MOM6-GPU/MOM6#95

Merged

Hallberg-NOAA reviewed Sep 2, 2026

View reviewed changes

Comment thread src/framework/MOM_intrinsic_functions.F90

Hallberg-NOAA reviewed Sep 2, 2026

View reviewed changes

Comment thread src/framework/intrinsics/MOM_exp_data_n128.F90 Outdated

Hallberg-NOAA reviewed Sep 2, 2026

View reviewed changes

Comment thread src/framework/testing/MOM_intrinsic_functions_tests.F90

Hallberg-NOAA reviewed Sep 2, 2026

View reviewed changes

Comment thread src/framework/intrinsics/MOM_exp.F90

adcroft requested a lite review from Copilot September 3, 2026 14:20

Copilot started reviewing on behalf of adcroft September 3, 2026 14:21 View session

Copilot AI reviewed Sep 3, 2026

View reviewed changes

Comment thread ac/m4/mom6_fc_fast_rint.m4

Comment thread src/framework/testing/MOM_intrinsic_functions_tests.F90

Comment thread src/framework/testing/MOM_intrinsic_functions_tests.F90

Comment thread ac/makedep Outdated

Comment thread src/framework/intrinsics/MOM_exp.F90

marshallward force-pushed the exp-repro branch from a655d4b to 6432d45 Compare September 3, 2026 14:59

marshallward added 2 commits September 3, 2026 11:18

marshallward force-pushed the exp-repro branch 2 times, most recently from 863565d to 5950dc4 Compare September 3, 2026 18:14

marshallward added 5 commits September 3, 2026 14:36

marshallward force-pushed the exp-repro branch from 5950dc4 to 166726b Compare September 3, 2026 18:37

Merge branch 'dev/gfdl' into exp-repro

7af9fab

adcroft approved these changes Sep 14, 2026

View reviewed changes

Navigation Menu

Uh oh!

exp_repro: Reproducible implementation of exp() - #1164

exp_repro: Reproducible implementation of exp()#1164
marshallward wants to merge 8 commits into
NOAA-GFDL:dev/gfdlfrom
marshallward:exp-repro

marshallward commented Aug 18, 2026 •

edited

Loading

Uh oh!

adcroft commented Aug 31, 2026

Uh oh!

adcroft left a comment

Uh oh!

marshallward commented Aug 31, 2026

Uh oh!

marshallward commented Sep 1, 2026 •

edited

Loading

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

marshallward commented Sep 2, 2026 •

edited

Loading

Uh oh!

Copilot AI left a comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

codecov Bot commented Sep 3, 2026 •

edited

Loading

Uh oh!

adcroft left a comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Uh oh!

Conversation

marshallward commented Aug 18, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

adcroft commented Aug 31, 2026

Uh oh!

adcroft left a comment

Choose a reason for hiding this comment

Uh oh!

marshallward commented Aug 31, 2026

Uh oh!

marshallward commented Sep 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

marshallward commented Sep 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

🟡 Changes recommended

Suppressed comments (2)

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

codecov Bot commented Sep 3, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

adcroft left a comment

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

marshallward commented Aug 18, 2026 •

edited

Loading

marshallward commented Sep 1, 2026 •

edited

Loading

marshallward commented Sep 2, 2026 •

edited

Loading

codecov Bot commented Sep 3, 2026 •

edited

Loading