diff --git a/image_analysis/scale_merge/FrenchWilson.cpp b/image_analysis/scale_merge/FrenchWilson.cpp
index 8b32ff824..7d105b021 100644
--- a/image_analysis/scale_merge/FrenchWilson.cpp
+++ b/image_analysis/scale_merge/FrenchWilson.cpp
@@ -24,8 +24,23 @@ struct Posterior {
 // so the exponentials never overflow/underflow. acentric: p(J) ~ exp(-J/S); centric:
 // p(J) ~ exp(-J/2S)/sqrt(J).
 // `logw` is caller-owned scratch of npts doubles (one per worker), so the integration allocates nothing.
+// The quadrature grid is j = (i + 1/2) * dj, so sqrt(j) and log(j) separate into a per-reflection
+// factor and a term that depends only on i: sqrt(j) = sqrt(dj) * sqrt(i + 1/2) and
+// log(j) = log(dj) + log(i + 1/2). The i-dependent halves are the same for every reflection, so
+// they are tabulated once instead of being recomputed npts times per intensity.
+struct QuadratureGrid {
+    std::vector<double> sqrt_i, log_i;
+    explicit QuadratureGrid(int npts) : sqrt_i(npts), log_i(npts) {
+        for (int i = 0; i < npts; ++i) {
+            const double x = i + 0.5;
+            sqrt_i[i] = std::sqrt(x);
+            log_i[i] = std::log(x);
+        }
+    }
+};
+
 Posterior integrate_posterior(double I, double sigma, double sigma_wilson, bool centric, int npts,
-                              std::vector<double> &logw) {
+                              std::vector<double> &logw, const QuadratureGrid &grid) {
     const double inv_2s2 = 1.0 / (2.0 * sigma * sigma);
     // The posterior is the Gaussian likelihood tilted by the exponential prior, so it peaks at
     // I - sigma^2/S and decays over whichever of sigma and S is TIGHTER. Ranging to I + 8 sigma
@@ -39,25 +54,34 @@ Posterior integrate_posterior(double I, double sigma, double sigma_wilson, bool
     const double j_max = peak + 10.0 * width;
     const double dj = j_max / npts;
 
+    const double log_dj = centric ? std::log(dj) : 0.0;
     double max_logw = -std::numeric_limits<double>::infinity();
     for (int i = 0; i < npts; ++i) {
         const double j = (i + 0.5) * dj;
         const double diff = I - j;
-        const double log_prior = centric ? (-j / (2.0 * sigma_wilson) - 0.5 * std::log(j))
-                                          : (-j / sigma_wilson);
+        const double log_prior = centric
+                                     ? (-j / (2.0 * sigma_wilson) - 0.5 * (log_dj + grid.log_i[i]))
+                                     : (-j / sigma_wilson);
         logw[i] = log_prior - diff * diff * inv_2s2;
         max_logw = std::max(max_logw, logw[i]);
     }
 
+    // exp(-37) is 8e-17: a point that far below the peak cannot change a normalised sum of
+    // doubles, and the posterior is sharply peaked, so most of the grid is skipped outright.
+    constexpr double LOG_NEGLIGIBLE = -37.0;
+    const double sqrt_dj = std::sqrt(dj);
     double sum_w = 0, sum_wI = 0, sum_wF = 0;
     for (int i = 0; i < npts; ++i) {
+        const double shifted = logw[i] - max_logw;
+        if (shifted < LOG_NEGLIGIBLE)
+            continue;
         const double j = (i + 0.5) * dj;
-        const double w = std::exp(logw[i] - max_logw);
+        const double w = std::exp(shifted);
         if (!std::isfinite(w))
             continue;
         sum_w += w;
         sum_wI += w * j;
-        sum_wF += w * std::sqrt(j);
+        sum_wF += w * (sqrt_dj * grid.sqrt_i[i]);
     }
     if (sum_w <= 0.0) {
         const double j = std::max(I, 0.0);
@@ -128,16 +152,17 @@ void ApplyFrenchWilson(std::vector<MergedReflection> &merged, const gemmi::Space
 
     // French-Wilson |F| for one intensity of reflection r (its mean, or one Bijvoet hand); the shell
     // Wilson prior, epsilon and centric flag are the reflection's, shared by all three.
+    // The Wilson prior and the centric flag belong to the reflection, not to the intensity, so they
+    // are looked up once and shared by its mean and both Bijvoet hands - three symmetry lookups
+    // became one.
     auto fw_one = [&](const MergedReflection &r, float I, float sigma, float &F, float &sigF,
-                      std::vector<double> &logw) {
+                      std::vector<double> &logw, const QuadratureGrid &grid, double sigma_wilson,
+                      bool centric) {
         if (!std::isfinite(I) || !std::isfinite(sigma) || sigma <= 0.0f) { naive_one(I, sigma, F, sigF); return; }
         // Strong reflections: the FW correction is negligible, <|F|> = sqrt(I).
         if (I > opts.strong_cutoff * sigma) { naive_one(I, sigma, F, sigF); return; }
-        const auto s = shells.GetShell(r.d);
-        const double sigma_wilson = epsilon(r) * (s ? shell_mean[*s] : global_mean);
-        const bool centric = gops.is_reflection_centric({{r.h, r.k, r.l}});
         const Posterior post = integrate_posterior(I, sigma, sigma_wilson, centric,
-                                                   opts.integration_points, logw);
+                                                   opts.integration_points, logw, grid);
         F = static_cast<float>(post.mean_F);
         sigF = static_cast<float>(std::sqrt(std::max(0.0, post.mean_I - post.mean_F * post.mean_F)));
     };
@@ -146,13 +171,17 @@ void ApplyFrenchWilson(std::vector<MergedReflection> &merged, const gemmi::Space
     const int n = static_cast<int>(merged.size());
     const int nt = std::clamp(opts.num_threads, 1, n);
     const int chunk = (n + nt - 1) / nt;
+    const QuadratureGrid grid(opts.integration_points);
     auto do_chunk = [&](int lo, int hi) {
         std::vector<double> logw(opts.integration_points);
         for (int i = lo; i < hi; ++i) {
             MergedReflection &r = merged[i];
-            fw_one(r, r.I, r.sigma, r.F, r.sigmaF, logw);
-            fw_one(r, r.I_plus, r.sigma_plus, r.F_plus, r.sigmaF_plus, logw);
-            fw_one(r, r.I_minus, r.sigma_minus, r.F_minus, r.sigmaF_minus, logw);
+            const auto s = shells.GetShell(r.d);
+            const double sigma_wilson = epsilon(r) * (s ? shell_mean[*s] : global_mean);
+            const bool centric = gops.is_reflection_centric({{r.h, r.k, r.l}});
+            fw_one(r, r.I, r.sigma, r.F, r.sigmaF, logw, grid, sigma_wilson, centric);
+            fw_one(r, r.I_plus, r.sigma_plus, r.F_plus, r.sigmaF_plus, logw, grid, sigma_wilson, centric);
+            fw_one(r, r.I_minus, r.sigma_minus, r.F_minus, r.sigmaF_minus, logw, grid, sigma_wilson, centric);
         }
     };
     if (nt == 1) {