Merge branch 'hotfix/2.0.2'

Author: Ozaq

VERSION

@@ -1 +1 @@
-2.0.1
+2.0.2

src/eckit/geo/util/gaussian_latitudes.cc

@@ -35,8 +35,8 @@ const std::vector<double>& gaussian_latitudes(size_t N, bool increasing) {
         return cache[key];
     }
 
-    std::vector<double> lats(2 * N);
 
+    std::vector<double> lats(2 * N);
 
     // Fourier coefficients of series expansion for the ordinary Legendre polynomials
     std::vector<double> zzfn(N + 1);
@@ -62,25 +62,29 @@ const std::vector<double>& gaussian_latitudes(size_t N, bool increasing) {
     }
 
 
-    // Newton loop (per latitude) to find 0 of Legendre polynomial of degree N (GAWL)
+    /*
+     * Newton loop (per latitude) to find 0 of Legendre polynomial of degree N (GAWL)
+     * The convergence is run more than once, for historical reasons.
+     */
     constexpr size_t Nmax = 20;
     constexpr auto eps    = std::numeric_limits<double>::epsilon() * 1000.;
 
     for (size_t i = 0; i < N; ++i) {
         // First guess for colatitude [rad]
-        double z = static_cast<double>(4 * (i + 1) - 1) * M_PI / static_cast<double>(4 * 2 * N + 2);
-        double x = (z + 1. / (std::tan(z) * static_cast<double>(8 * (2 * N) * (2 * N))));
+        auto z = static_cast<double>(4 * (i + 1) - 1) * M_PI / static_cast<double>(4 * 2 * N + 2);
+        auto x = z + 1. / (std::tan(z) * static_cast<double>(8 * (2 * N) * (2 * N)));
 
         auto converged = false;
 
         for (size_t n = 0; n < Nmax; ++n) {
-            auto f  = 0.5 * zzfn[0];  // normalised ordinary Legendre polynomial == \overbar{P_n}^0
-            auto fp = 0.;             // normalised derivative == d/d\theta(\overbar{P_n}^0)
-
-            for (size_t i = 1; i <= N; ++i) {
-                const auto i2 = static_cast<double>(i * 2);
-                f += zzfn[i] * std::cos(i2 * x);
-                fp -= zzfn[i] * std::sin(i2 * x) * i2;
+            auto f = zzfn[0] * 0.5 +
+                     zzfn[1] * std::cos(2. * x);         // normalised ordinary Legendre polynomial == \overbar{P_n}^0
+            auto fp = zzfn[1] * std::sin(2. * x) * -2.;  // normalised derivative == d/d\theta(\overbar{P_n}^0)
+
+            for (size_t j = 2; j <= N; ++j) {
+                const auto j2 = static_cast<double>(j * 2);
+                f += zzfn[j] * std::cos(j2 * x);
+                fp -= zzfn[j] * std::sin(j2 * x) * j2;
             }
 
             auto dx = -f / fp;

src/eckit/geo/util/gaussian_quadrature_weights.cc

@@ -49,7 +49,7 @@ const std::vector<double>& gaussian_quadrature_weights(size_t N) {
                 zfn[i] *= std::sqrt(1. - 0.25 / (static_cast<double>(j * j)));
             }
 
-            for (size_t j = 2; j <= i - i % 2; j += 2) {
+            for (size_t j = 2; j <= i - (i % 2); j += 2) {
                 zfn[i - j] = zfn[i - j + 2] * static_cast<double>((j - 1) * (2 * i - j + 2)) /
                              static_cast<double>(j * (2 * i - j + 1));
             }
@@ -61,34 +61,37 @@ const std::vector<double>& gaussian_quadrature_weights(size_t N) {
     }
 
 
-    // Newton loop (per latitude) to find 0 of Legendre polynomial of degree N (GAWL)
+    /*
+     * Newton loop (per latitude) to find 0 of Legendre polynomial of degree N (GAWL)
+     * The convergence is run more than once, for historical reasons.
+     */
     constexpr size_t Nmax = 20;
-    constexpr double eps  = std::numeric_limits<double>::epsilon() * 1000.;
+    constexpr auto eps    = std::numeric_limits<double>::epsilon() * 1000.;
 
     for (size_t i = 0; i < N; ++i) {
-        // First guess
+        // First guess for colatitude [rad]
         auto z = static_cast<double>(4 * (i + 1) - 1) * M_PI / static_cast<double>(4 * 2 * N + 2);
         auto x = z + 1. / (std::tan(z) * static_cast<double>(8 * (2 * N) * (2 * N)));
 
         double w       = 0.;
-        bool converged = false;
+        auto converged = false;
 
         for (size_t n = 1; n < Nmax; ++n) {
-            auto f  = 0.5 * zzfn[0];  // normalised ordinary Legendre polynomial == \overbar{P_n}^0
-            auto fp = 0.;             // normalised derivative == d/d\theta(\overbar{P_n}^0)
-
-            for (size_t j = 1; j <= N; ++j) {
-                const auto i2 = static_cast<double>(j * 2);
-                f += zzfn[j] * std::cos(i2 * x);
-                fp -= zzfn[j] * std::sin(i2 * x) * i2;
+            auto f = zzfn[0] * 0.5 +
+                     zzfn[1] * std::cos(2. * x);         // normalised ordinary Legendre polynomial == \overbar{P_n}^0
+            auto fp = zzfn[1] * std::sin(2. * x) * -2.;  // normalised derivative == d/d\theta(\overbar{P_n}^0)
+
+            for (size_t j = 2; j <= N; ++j) {
+                const auto j2 = static_cast<double>(j * 2);
+                f += zzfn[j] * std::cos(j2 * x);
+                fp -= zzfn[j] * std::sin(j2 * x) * j2;
             }
 
             auto dx = -f / fp;
             x += dx;
 
             if (converged) {
-                const auto j = 2 * N - 1 - i;
-                w = weights[i] = static_cast<double>(2 * 2 * N + 1) / (fp * fp);
+                w = static_cast<double>(2 * 2 * N + 1) / (fp * fp);
                 break;
             }
 

tests/geo/CMakeLists.txt

@@ -4,6 +4,7 @@ foreach(_test
         cache
         figure
         figure_earth
+        gaussian
         great_circle
         grid
         grid_healpix

tests/geo/gaussian.cc

@@ -0,0 +1,70 @@
+/*
+ * (C) Copyright 1996- ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ *
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+
+#include "eckit/geo/util.h"
+#include "eckit/maths/FloatingPointExceptions.h"
+#include "eckit/testing/Test.h"
+#include "eckit/types/FloatCompare.h"
+
+
+#define EXPECT_APPROX(a, b) EXPECT(::eckit::types::is_approximately_equal<double>((a), (b), 1e-6))
+
+
+namespace eckit::geo::test {
+
+
+/**
+ * This case verifies that Floating Point Exceptions (FPEs) don't trigger for the functions:
+ * - eckit::geo::util::gaussian_latitudes
+ * - eckit::geo::util::gaussian_quadrature_weights
+ *
+ * FPE trapping is enabled to ensure any FPE-triggering operations are caught.
+ * @see ECKIT-674
+ */
+CASE("gaussian_latitudes") {
+    const auto& lats_1 = geo::util::gaussian_latitudes(1, true);
+
+    EXPECT(lats_1.size() == 2);
+    EXPECT_APPROX(lats_1.front(), -35.264390);
+    EXPECT_APPROX(lats_1.front(), -lats_1.back());
+
+    const auto& lats_2 = geo::util::gaussian_latitudes(80, false);
+
+    EXPECT(lats_2.size() == 2 * 80);
+    EXPECT_APPROX(lats_2.front(), 89.141519);
+    EXPECT_APPROX(lats_2.front(), -lats_2.back());
+}
+
+
+CASE("gaussian_quadrature_weights") {
+    const auto& quad_1 = geo::util::gaussian_quadrature_weights(1);
+
+    EXPECT(quad_1.size() == 2);
+    EXPECT_APPROX(quad_1.front(), 0.5);
+    EXPECT_APPROX(quad_1.front(), quad_1.back());
+
+    const auto& quad_2 = geo::util::gaussian_quadrature_weights(80);
+
+    EXPECT(quad_2.size() == 2 * 80);
+    EXPECT_APPROX(quad_2.front(), 0.000144);
+    EXPECT_APPROX(quad_2.front(), quad_2.back());
+}
+
+
+}  // namespace eckit::geo::test
+
+
+int main(int argc, char** argv) {
+    eckit::maths::FloatingPointExceptions::enable_floating_point_exceptions();
+
+    return eckit::testing::run_tests(argc, argv);
+}