Merge pull request #3 from axionbuster/perf1

Perf1

dyn/newton.h

@@ -102,17 +102,15 @@ template <typename F = float, unsigned short N_MONTE = 30> class Gravity {
       auto q = c1 - r0 * p;
       auto r = std::abs(q);
       if (r > c1.radius) {
-        // Give up a bit of accuracy for speed.
         auto s = F(1) / r;
-        // Cancel out the non-radial component.
-        auto dB = s * s * s * (c1 - r0 * p.real());
+        auto dB = s * s * s * (c1 - r0 * p);
         B += dB;
       }
     }
 
     // Converting constexpr integer N_MONTE to floating point constexpr is
     // helpful because it reduces burden due to type conversion on CPU.
-    return F(0.5) / F(N_MONTE) * m1 * G * B;
+    return F(1) / F(N_MONTE) * m1 * G * B;
   }
 };
 

tests/newton_test.cpp

@@ -42,6 +42,7 @@ TEST_F(NewtonSuite, PassThrough0) {
   auto constexpr S2 = std::numbers::sqrt2_v<float>;
   dyn::Yoshida<> yoshi{c0 + S2 + S2 * 1.0if, 0};
 
+  // 90 steps : 1 second = 1800 steps : 20 seconds
   auto constexpr steps = 1800;
   auto accel = [&](auto xy) { return this->common_accel(xy, c0); };
   for (int i = 0; i < steps; i++) {
@@ -136,5 +137,5 @@ TEST_F(NewtonSuite, Figure8) {
                  std::abs(yoshi2.y0 - yoshis[2].y0)};
 
   for (int i = 0; i < 3; i++)
-    ASSERT_NEAR(1.0f, qs[i], 0.1f) << "(i = " << i << ")";
+    ASSERT_NEAR(0.0f, qs[i], 0.1f) << "(i = " << i << ")";
 }