[ggorman]

FD weights are rounded with evalf(9) regardless of dtype, so fp64 operators get the same 9-digit literals as fp32 ones (e.g. 2.53968254e-2 for 8/315, 1.196289060 for 1.1962890625). In double precision this sets an accuracy floor:

• an fp64 u.dx2 of x**2 is off by ~5e-6 (should be round-off, ~5e-13);
• an MMS of u.laplace at space order 8 stalls at 1e-7 and gets worse as h shrinks (1.2e-4 at n=1024, vs 7.7e-11 with exact weights);
• a 200-step fp64 acoustic run differs from NumPy with exact weights by 9e-8 instead of 1e-14.

Two changes:

1. numeric_weights turns float literals in x0 (e.g. x + 0.5*h_x) into their exact Rational value before computing the Taylor weights, like the Rational(s) already used in make_shift_x0. This removes the float noise that the 9-digit rounding was hiding, with no tolerance involved.
2. The precision follows the expression dtype: 9 digits as before for float32 (and anything else), 20 for float64/complex128, which keeps the literals within one ulp of the exact weights. It applies to every weight source (taylor, custom weights=, registry rules such as drp).

fp32 generated code is byte-identical to main for the seismic examples (acoustic fwd/adj/grad/born, TTI, elastic, viscoelastic, viscoacoustic; space orders 4-16). The only exception I found is div(v, shift=(.5, .5, .5)), where main emits a spurious - 3.05311332e-16F*v_x[...] term that is now gone. Op counts are otherwise unchanged in fp64 too, and no existing expected strings had to change.

This follows the review on #2977: no float-based tolerance, no extra helper or flag, and custom/drp weights are covered.

Tests: test_shifted_div/test_shifted_grad also run in float64 (they fail without change 1), and test_dx2_exact_double fails on main and passes here.

Fixes #2976

🤖 Generated with Claude Code

• Adds -ffp-contract=off to GCC's safe-math flags: the C++ backend was contracting multiply-adds under
  safe-math while the C backend was not, which is what failed test_gradient_equivalence in the two C++ jobs
  (see the comment below for the numbers).

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 84.15%. Comparing base (078bfe3) to head (bbb5de2).

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #3074   +/-   ##
=======================================
  Coverage   84.15%   84.15%           
=======================================
  Files         258      258           
  Lines       57267    57290   +23     
  Branches     4893     4893           
=======================================
+ Hits        48191    48214   +23     
+ Misses       8249     8248    -1     
- Partials      827      828    +1     

Flag	Coverage Δ
pytest-gpu-aomp-amdgpuX	68.41% <85.71%> (+<0.01%)	⬆️
pytest-gpu-gcc-	78.97% <100.00%> (+0.01%)	⬆️
pytest-gpu-icx-	78.91% <97.50%> (+0.01%)	⬆️
pytest-gpu-nvc-nvidiaX	68.99% <85.71%> (-0.02%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ 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.

[ggorman]

The two C++ CI-core failures on 43c5113 were not caused by the 20-digit literals. The literals and their hex
values are identical in the C and C++ code. What differs is contraction: GCC turns -ffp-contract off by
default only in ISO C mode (-std=c99), so the C++ backend (-std=c++14) has been fusing multiply-adds under
safe-math all along, while the C backend has not.

test_gradient_equivalence compares sum(-u.dt2*v) with sum(-u*v.dt2) per point at atol=rtol=1e-13. Without
FMA, each tap of u.dt2 = -2*r1*u0 + r1*um + r1*up is rounded once and the two additions are exact (same binade),
so the error is the second difference of a consistently rounded sequence and largely cancels under the sum.
g++ fuses that into fma(r1, up, fma(-2*r1, u0, r1*um)) and the cancellation is gone. The test's worst-point
ratio to its tolerance (above 1 fails):

build	main	this branch	this branch + -ffp-contract=off
C, gcc-11	0.33	0.45	0.45
CXX, g++-11	0.98	1.30	0.45 (bitwise = C)
CXXopenmp, g++-11	0.66	1.13	0.45 (bitwise = C)
C with forced -ffp-contract=fast	1.50

So main was passing the C++ jobs with 2% headroom, and any perturbation re-rolls the statistic: across 17
variants (digits 9..30, tn 700..800) 13 fail with FMA and none without; a one-ulp change in the damping mask
alone moves the ratio between 0.85 and 1.30. The same numbers reproduce on a Ryzen 5900X with gcc-11 and
gcc-12.

The last commit adds -ffp-contract=off to the GNU safe-math flags, so C++ safe-math behaves like C
safe-math (no contraction). It only applies when safe-math is on; -ffast-math builds keep FMA. Cost under
safe-math on CPU is within noise (3D acoustic, 384^3, so 8/16, fp32/fp64, 16 threads on a Xeon 8480+ and 12
on a Ryzen 5900X: all within 2%, these kernels are bandwidth-bound).

Two things this does not cover: nvc++ and nvcc also contract under safe-math (the same test fails on OpenACC
with -fast and on CUDA by default, and passes with -Mnofma / --fmad=false; there, fp64 FMA is worth 13-39%
on an RTX 3090), and the test's 1e-13 per-point tolerance rests on this exact-cancellation structure rather than
on an error bound. Raised as #3077 and #3078; #3079 is the GPU CI never running tests/test_gradient.py, and
devitocodespro/devitopro#1028 a DevitoPRO-side problem found on the way.

[mloubout]

needs rebase then gtg

[FabioLuporini]

I'd rather

try:
    precision = _PRECISION[expr.dtype]
except KeyError:
    raise .....
weights = [sympify(scale * w).evalf(precision) for w in weights]

and regardless, no hardcoded magic numbers (eg 9 here)

[mloubout]

I'm sorry @FabioLuporini but no that's not better and that create a mess. THere is a bunch of cases where we have different dtypes than the standard one (with fp16 to start with) and the 9 is not magic it just defaults to fp32 precision,

[ggorman]

Integers do reach this site (13 tests, plus PRO's test_weights_are_real on int32/int64/float16), so a raise is a trap, but the inline 9 was duplicating the fp32 entry. Proposal: precision = _PRECISION.get(expr.dtype, _PRECISION[np.float32]), i.e. the mapper's single-precision digits for anything not listed, no raise, and a comment saying so. The current head has the strict version; I'll swap it for this if you're both happy.

[mloubout]

precision = _PRECISION.get(expr.dtype, _PRECISION[np.float32])

yes

[ggorman]

Swapped the strict lookup for _PRECISION.get(expr.dtype, _PRECISION[np.float32]) as agreed above (bbb5de2): no raise, float16 entry kept, unsupported-dtype test dropped. All 42 checks green; the osx job needed one rerun for the Apple-silicon detection miss that #3048 fixes.